Lévy Diffusion Under Power-Law Stochastic Resetting

1. Introduction

In nature, random motion manifests ubiquitously across microscopic to macroscopic scales, with Lévy walk emerging as a fundamental paradigm alongside Brownian motion. Characterized by heavy-tailed jump-length distributions, Lévy walk facilitates efficient long-range exploration through its distinctive combination of localized search and sporadic abrupt relocation [1]. This dynamical strategy has been observed in diverse biological systems [2,3,4], including albatross foraging patterns [2], bacterial chemotaxis [3], and cell migration in vivo [4], and theoretically validated as an optimal search mechanism in sparse environments [5,6]. Unlike Brownian motion, whose diffusive behavior is governed by short-range fluctuations, Lévy walk’s power-law statistics enables it to overcome the “curse of dimensionality” in target search, establishing it as a cornerstone for modeling transport in non-equilibrium systems [7].

When searching for a target in a crowd, if one cannot find it within an affordable time interval, then a more efficient way to complete the task is to go back to the beginning and start the process again. Stochastic resetting (SR) - defined as the stochastic re-initiation of a system to a reference state - represents another ubiquitous mechanism in natural search processes [8]. For instance, animal foraging often involves periodic returns to shelter or nesting sites [1], effectively “resetting” the search trajectory to balance exploration and resource exploitation [9]. Similarly, microbial cells navigating chemical gradients may periodically abandon fruitless paths via spontaneous resetting, a strategy that optimizes nutrient acquisition in heterogeneous environments [10,11]. Current research had recourse to SR mainly to optimize search strategies [9,12,13,14,15,16]. For example, Nagar and Gupta explored the optimal search of Brownian motion under SR with waiting times following a power-law distribution [12]; Kusmierz et al. investigated search strategy optimization for heavy-tailed Lévy flights combined with SR [13]. Such resetting events, whether driven by energy constraints, environmental cues, or internal behavioral rhythms, impose a hierarchical organization on random walks, transforming diffusive dynamics into controlled transport processes.

Recently, attention has shifted to the impact of SR on diffusion dynamics. Liu et al. investigated the diffusive behavior of Lévy walks under SR with rests [17]; Kuśmierz et al. [18] and Zhou et al. [19] studied the subdiffusive dynamics of continuous-time random walks (CTRW) under exponentially distributed stochastic resetting times, respectively. Bodrova and Sokolov explored diffusion under power-law resetting [20], but focused on CTRWs with finite step lengths, thereby neglecting the interplay between long-range jumps and power-law resetting. Targeting the long-range jump properties of Lévy walks, our study examines their behavior under power-law resetting [12,20,21], aiming to reveal novel phenomena beyond the scope of CTRW models [18,19,20]. Indeed, in natural systems, processes such as animal foraging [2] and cellular transport [4] exhibit both Lévy-type long-distance migration and power-law resetting (e.g., periodic returns to nesting sites). Our approach provides quantitative explanations for such phenomena, adding to a more realistic representation of biological and physical system dynamics compared to CTRW models.

2. Lévy Walks Under Stochastic Resetting

One dimensional (1D) unbiased Lévy walks constitute a class of random walks characterized by step lengths distributed according to a Lévy distribution [7]. Analytically, symmetric Lévy distributions are defined through their characteristic function,

$$F\left(k\right)=exp\left[-{c\left|k\right|}^{\alpha }\right],0<\alpha \le 2,$$

(1)

where the scale parameter c may be conveniently set to unity ($c\equiv 1$) and the (Lévy) index $\alpha$ spans the interval $(0,2]$. This characteristic function corresponds to analytical expressions in the spatial domain only for three specific $\alpha$ values: Lévy distributions for $\alpha =1/2$, Cauchy distributions for $\alpha =1$, and Gaussian distributions for $\alpha =2$.

Lévy-distributed random number generation necessitates numerical methods such as the Chambers–Mallows–Stuck algorithm [22] or the Mantegna algorithm [23]. The latter generates Lévy-distributed step lengths through two sets of normally distributed random numbers, providing simplified formulation and computational efficiency, though limited to symmetric distributions, such as those in Eq. (1). Consequently, we had recourse to the Mantegna algorithm for all numerical simulations presented herein.

We conducted systematic MSD simulations for Lévy walks and compared numerical results with theoretical predictions. For each generated step length $s_i$, we calculated the corresponding step duration ${\Delta }_i=\left|s_i\right|/v$ for a particle with constant speed v. In our simulations, we conventionally set $v\equiv 1$. Iterative application of this procedure generates 1D Lévy walk trajectories (Figure 1). Specifically, we fixed the Lévy index $\alpha$ and generated $N={10}^6$ independent trajectories originating from the initial condition $x\left(0\right)=0$. For each trajectory, we recorded its position-time series and computed the ensemble-averaged MSD

$$\langle x^2\left(t\right)\rangle ={\displaystyle \frac{1}{N}}\sum _{k=1}^N{|x_k\left(t\right)-x_k\left(0\right)|}^2.$$

(2)

For an untruncated Lévy walk, the asymptotic MSD [i.e., the $t\to \infty$ limit of Eq. (2)] is known to scale like [24,25]

$$M_0\left(t\right)=\underset{t\to \infty }{lim}\langle x^2\left(t\right)\rangle \sim \left\{\begin{array}{cc}(1-\alpha )t^2,\hfill & 0<\alpha <1,\hfill \\ {\displaystyle \frac{t^2}{lnt}},\hfill & \alpha =1,\hfill \\ 2F_{\alpha }{t}^{3-\alpha },\hfill & 1<\alpha <2,\hfill \\ 2t,\hfill & \alpha =2,\hfill \end{array}\right.$$

(3)

with $F_{\alpha }={|\Gamma (1-\alpha )|(\alpha -1)}^2/[\alpha \Gamma (4-\alpha )]$ and $\Gamma (\dots )$ denoting the Gamma function. We note that by setting $c=v=1$, we implicitly established dimensionless units for both x and t.

Except for $\alpha =1$, the 1D asymptotic Lévy walk diffusion can be formulated as

$$M_0\left(t\right)=K{t}^{{\gamma }_0},$$

(4)

where the free diffusion prefactor K and exponent ${\gamma }_0$ are functions of $\alpha$ that can be directly obtained from Eq. (3). However, in view of Eq. (1), at short times, $ct<1$ (i.e., $t<1$), the particle propagates at constant speed, $v=1$, yielding $\langle x^2\left(t\right)\rangle \simeq t^2$ independently of $\alpha$, as demonstrated by the short-time dat fits in Fig. Figure 2. This effect is an artifact of our Lévy step generator and is not further considered in the following.

We employed Eq. (4) to fit simulated MSD curves and extract numerical estimates for K and ${\gamma }_0$, as summarized in Table 1. Notwithstanding increased systematic errors in the Mantegna algorithm near $\alpha =1$ and 2, all relative errors remain below $20\%$, confirming numerical reliability. In particular, simulated exponents consistently match theoretical predictions across repeated trials, validating our numerical code for subsequent Lévy walk simulations under SR conditions.

Under SR, the Lévy walker is reset to the initial position $X_r$ at random time intervals $\tau$ (as illustrated in Figure 1) drawn from a distribution function $\psi \left(\tau \right)$, to be specified. In 1D, its coordinate can be expressed as

$$\begin{array}{c}\hfill x\left(t\right)=X_r+\tilde{x}(t-t_{n-1})\mathrm{for}t\in [t_{n-1},t_n),\end{array}$$

(5)

where $t_{n-1}$ and $t_n$ denote the times of two consecutive random resetting events and $\tilde{x}\left(t\right)$ represents the same Lévy walk reinitialized with $\tilde{x}\left(0\right)=0$ (Method I in notation of Ref. [18]) and truncated after the resetting interval ${\tau }_n=t_n-t_{n-1}$. In our simulations, we set $X_r=0$ and, contrary to the authors of Ref. [26], we allowed resetting to occur at any time, i.e., also during Lévy steps (continuous SR).

As an Ansatz for our analysis, we postulate that the asymptotic MSD under resetting obeys the same scaling law, Eq. (4), as for a free Lévy walker, that is

$$M\left(t\right)=\underset{t\to \infty }{lim}\langle x^2\left(t\right)\rangle =K^{\prime }{t}^{\gamma },$$

(6)

with appropriate prefactor $K^{\prime }$ and diffusion exponent $\gamma$ to be determined based on the specific choice of $\psi \left(\tau \right)$. Logarithmic corrections to the aforementioned scaling law, occurring for $\alpha =1$, according to Eq. (3), and for certain power laws of the $\tau$ distributions, will be treated separately.

3. Lévy Diffusion Under Exponential Stochastic Resetting

In this section, we revisit Lévy walks subject to SR with exponentially distributed resetting times [17,18,19]. Our goal is twofold: first, to develop and validate the renewal theory framework, which will be extended in Section 4 to the more complex case of power-law resetting; and second, to use exponential resetting as a benchmark to illustrate how the interplay between the power-law statistics of Lévy walk steps and resetting times gives rise to distinct diffusive regimes.

We begin by assuming that the resetting times, $\tau$, follow an exponential distribution,

$$\psi \left(\tau \right)={\displaystyle \frac{1}{{\tau }_e}}{e}^{-\tau /{\tau }_e},$$

(7)

with finite first and second moments, respectively,

$$\langle \tau \rangle ={\tau }_e\mathrm{and}\langle {\tau }^2\rangle =2{\tau }_e^2.$$

(8)

Figure 2 illustrates how exponential resetting influences the temporal evolution of Lévy walk MSD for various $\alpha$ values. Colored symbols represent MSD data for different resetting parameters ${\tau }_e$, while black symbols denote corresponding MSD in the absence of resetting. For $t\ll {\tau }_e$, low resetting probability results in MSD closely approximating free diffusion. More significantly, for $t\sim {\tau }_e$, frequent resetting induces MSD curve deviations from free growth, ultimately leading to horizontal saturation plateaus for $t\gg {\tau }_e$. Plateau heights depend on both the Lévy index $\alpha$ and resetting parameter ${\tau }_e$. This phenomenon demonstrates that regardless of underlying dynamics (ballistic, superdiffusive, etc.), exponential stochastic resetting establishes an upper bound for MSD at long times, preventing unbounded diffusion [9,15].

Exact expressions for the dependence of Lévy walk MSD on the tunable parameters $\alpha$ and ${\tau }_e$ are derived using a renewal theory method and numerically validated in the following subsection.

3.1. Scaling Analysis and Diffusion Transition Times

By handling the resetting process of Eq. (5) as a renewal process [27], the contribution of each resetting event to the MSD can be accounted for individually. The renewal theory provides a concise and straightforward analytical framework. The general renewal equation for the MSD in the asymptotic regime is

$$M\left(t\right)=S\left(t\right)M_0\left(t\right)+{\int }_0^t\psi \left(\tau \right)M(t-\tau )d\tau ,$$

(9)

where $M_0\left(t\right)$ is the known free MSD of Eq. (3) and $\psi \left(\tau \right)$ is an arbitrary resetting time distribution function. $S\left(t\right)$ is the survival function associated with $\psi \left(\tau \right)$, i.e., the probability that resetting happens after a waiting time larger than t,

$$S\left(t\right)={\int }_t^{\infty }\psi \left(\tau \right)d\tau .$$

(10)

When the mean $\langle \tau \rangle$ of the resetting time distribution $\psi \left(\tau \right)$ exists (i.e., is finite), one can make use of the closed form of the asymptotic “age distribution”,

$$f_{\mathrm{age}}\left(t\right)={\displaystyle \frac{S\left(t\right)}{\langle \tau \rangle }},$$

(11)

to formally solve Eq. (9) for $t\gg \langle \tau \rangle$ [28], that is

$$M\left(t\right)={\int }_0^t{f}_{\mathrm{age}}\left(t^{\prime }\right)M_0\left(t^{\prime }\right)d{t}^{\prime }.$$

(12)

In this way, the asymptotic MSD is formulated in terms of the average Lévy diffusion within single resetting events.

Equation (12) provides the most direct way to calculate the plateau ${\mathrm{M}}_{\mathrm{st}}\equiv \underset{t\to \infty }{lim}M\left(t\right)$ for the asymptotic Lévy MSD under exponential resetting. From the exponential distribution in Eq. (7), we obtain the survival function $S\left(t\right)=e^{-t/{\tau }_e}$ and the corresponding age distribution $f_{\mathrm{age}}\left(t\right)=e^{-t/{\tau }_e}/{\tau }_e$. Accordingly, for $\alpha \ne 1$, Eq. (12) can be rewritten as

$${\mathrm{M}}_{\mathrm{st}}\equiv \underset{t\to \infty }{lim}M\left(t\right)={\displaystyle \frac{K}{{\tau }_e}}{\int }_0^{\infty }{t}^{{\gamma }_0}{e}^{-t/{\tau }_e}dt,$$

(13)

where we used the formal expression for $M_0\left(t\right)$ introduced in Eq. (4). This integral can be easily performed in terms of Gamma functions, to obtain

$${\mathrm{M}}_{\mathrm{st}}=K{\tau }_e^{{\gamma }_0}\Gamma (1+{\gamma }_0),\alpha \ne 1.$$

(14)

For $\alpha =1$, the Lévy walk MSD in Eq. (3) does not follow a power law, so we must handle it separately. Inserting the free MSD expression, $M_0\left(t\right)=t^2/lnt$ for $\alpha =1$, into Eq. (12), we end up with an integral expression for ${\mathrm{M}}_{\mathrm{st}}$, namely

$${\mathrm{M}}_{\mathrm{st}}={\int }_0^{\infty }{\displaystyle \frac{1}{{\tau }_e}}{e}^{-t/{\tau }_e}{\displaystyle \frac{t^2}{lnt}}dt.$$

(15)

[We recall that x and t are expressed in the dimensionless units corresponding to $c=v=1$.] This integral can be separated into two distinct integrals by introducing the auxiliary variable $y=t/{\tau }_e$, that is,

$${\mathrm{M}}_{\mathrm{st}}={\displaystyle \frac{{\tau }_e^2}{ln{\tau }_e}}{\int }_0^{\infty }\left(1-{\displaystyle \frac{lny}{ln{\tau }_e}}\right)y^2{e}^{-y}dy,$$

(16)

where ${\int }_0^{\infty }{y}^2{e}^{-y}dy=2$ and ${\int }_0^{\infty }{y}^{2}lny{e}^{-y}dy=2\psi \left(3\right)$, with $\psi \left(z\right)={\displaystyle \frac{d}{dz}}ln\Gamma \left(z\right)$. For an integer n, $\psi \left(n\right)=-{\gamma }_{\mathrm{E}}+{\sum }_{z=1}^{n-1}1/z,$ where ${\gamma }_{\mathrm{E}}\approx 0.5772$ is the Euler-Mascheroni constant; hence, for $n=3$, $\psi \left(3\right)\approx 1$. In conclusion, Eq. (16), can be approximated to

$$\begin{array}{cc}\hfill {\mathrm{M}}_{\mathrm{st}}& \simeq {\displaystyle \frac{2{\tau }_e^2}{ln{\tau }_e}}\left(1-{\displaystyle \frac{1}{ln{\tau }_e}}\right),\alpha =1.\hfill \end{array}$$

(17)

Deriving the superdiffusive regime of $M\left(t\right)$ is generally a more complex task, requiring the solution of the full renewal equation (9). In the case of exponential SR, we need to extract the time dependence of $M\left(t\right)$ for $t\ll {\tau }_e$. Upon Laplace transformation, Eq.(9) can be equivalently reformulated as

$$\tilde{M}\left(s\right)={\tilde{S}}_0\left(s\right)+\tilde{\psi }\left(s\right)\tilde{M}\left(s\right)⟹\tilde{M}\left(s\right)={\displaystyle \frac{{\tilde{S}}_0\left(s\right)}{1-\tilde{\psi }\left(s\right)}},$$

(18)

where ${\tilde{S}}_0\left(s\right)=\mathcal{L}\left\{S\left(t\right)M_0\left(t\right)\right\}\equiv \mathcal{L}\left\{S_0\left(t\right)\right\}$. The large-t (short-t) behavior of $M\left(t\right)$ can be readily determined by analyzing Eq. (18) in the limit of small (large) s.

Using Eq. (4) for $M_0\left(t\right)$ with $\alpha \ne 1$ and the explicit expression for $S\left(t\right)$ reported above, calculating $\tilde{M}\left(s\right)$ reduces to evaluating the Laplace transform of elementary functions (powers and exponentials); hence, $\tilde{M}\left(s\right)=K\Gamma (1+{\gamma }_0){\tau }_e^{{\gamma }_0}/\left[s{(1+s{\tau }_e)}^{{\gamma }_0}\right]$. Inverse transforming $\tilde{M}\left(s\right)$ yields

$$M\left(t\right)=M_0\left(t\right),t\ll {\tau }_e,$$

(19)

as anticipated in Fig. Figure 2. Of course, in the opposite limit, $t\gg {\tau }_e$ one recovers the plateau value of Eq. (14). The same conclusion is readily extended to the case of logarithmic corrections to the scaling equation (4) for $M_0\left(t\right)$, which corresponds to the marginal case $\alpha =1$ in Eq. (3).

To validate our analytical predictions, in Figure 3 we display the datasets from Fig. Figure 2 with t and MSD rescaled respectively by ${\tau }_e$ and ${\mathrm{M}}_{\mathrm{st}}$. For ${\mathrm{M}}_{\mathrm{st}}$, we utilized the appropriate analytical expressions from Eqs. (14) and (17), employing numerical values for K obtained from simulation data fitting [rather than theoretical values from Eq. (3)]. This approach ensures greater consistency, accounting for numerical inaccuracies inherent in Mantegna’s method for Lévy step generation (Table 1). For comparative purposes, we include in each panel (i.e., for all $\alpha$ values) the corresponding superdiffusive curve (black) in the absence of SR. For graphical convenience, t and the no-resetting MSD data have been rescaled by ${\tau }_e={10}^4$ and the corresponding ${\mathrm{M}}_{\mathrm{st}}$, respectively. All colored curves for finite ${\tau }_e$ converge asymptotically to $\mathrm{MSD}/{\mathrm{M}}_{\mathrm{st}}\approx 1$ for $t\gg {\tau }_e$, in close agreement with theory. This occurs because, when the observation time exceeds the mean SR time, ${\tau }_e$, accumulated reset events drive the system from free superdiffusion into a regime of reset-governed localization. In contrast, Lévy diffusion remains largely insensitive to SR for $t<{\tau }_e$. This establishes a robust reference framework for subsequent comparisons under varying reset distributions.

Finally, Figure 4 illustrates, for various Lévy indices $\alpha$, the dependence of ${\mathrm{M}}_{\mathrm{st}}$ on the mean resetting time ${\tau }_e$. Colored points represent simulation-measured asymptotic ${\mathrm{MSD}}_{\mathrm{st}}$; solid lines of the same color denote the corresponding theoretical predictions for ${\mathrm{M}}_{\mathrm{st}}$ from Eqs. (14) and (17). Our Lévy walk simulations result to closely overlap the theoretical curves for $M_0\left(t\right)$ in Eq. (3) only for $t\gtrsim 100$ (see Figure 2). Therefore, not surprisingly, substantial discrepancies between numerical data and analytical predictions are observed in Figure 4 for ${\tau }_e<100$, while excellent agreement is achieved for ${\tau }_e\gtrsim 100$ across the entire $\alpha$ range.

4. Lévy Diffusion Under Power-Law Stochastic Resetting

Resetting times in numerous complex systems are characterized by heavy-tailed distributions. For instance, animal foraging-pause durations often exhibit approximate power-law or multiscale characteristics [29], and inter-spike intervals in neuronal network activity may also demonstrate multiscale dynamics [30]. Therefore, we next examine the case of resetting times distributed according to a power-law.

Let the Lévy walker be now instantaneously reset to $X_r=0$ after a random time $\tau$ drawn from a Pareto (Type I) distribution,

$$\psi \left(\tau \right)=\left\{\begin{array}{cc}{\displaystyle \frac{\beta }{{\tau }_0}}{\left(\frac{{\tau }_0}{\tau }\right)}^{1+\beta },\hfill & \tau \ge {\tau }_0,\hfill \\ 0,\hfill & \tau <{\tau }_0,\hfill \end{array}\right.$$

(20)

where ${\tau }_0$ denotes a scale parameter and exponent $\beta$ characterizes the distribution tail decay. The $\tau$ first and second moments,

$$\begin{array}{cc}\hfill \langle \tau \rangle & ={\displaystyle \frac{\beta {\tau }_0}{\beta -1}},\beta >1\hfill \\ \hfill \langle {\tau }^2\rangle & ={\displaystyle \frac{\beta {\tau }_0^2}{\beta -2}},\beta >2,\hfill \end{array}$$

(21)

diverge for $\beta <1$ and $\beta <2$, respectively. In our simulations, unless otherwise specified, we assumed a small resetting timescale, namely ${\tau }_0=0.1$. The principal findings of this section can be readily generalized to Lomax distributions with identical ${\tau }_0$ and $\beta$ parameters, and extended in principle to other heavy-tailed $\tau$ distributions.

Figure 5 illustrates the temporal evolution of Lévy walk MSD under Pareto-distributed resetting times for various Lévy indices $\alpha$. Colored curves for different $\beta$ values are plotted in each panel; black curves for $\beta =0$ represent free Lévy diffusion. In each panel, MSD curves for $\beta <1$ run parallel to a black curve for $\beta =0$, exhibiting approximate free Lévy superdiffusive scaling. As $\beta$ increases, the diffusion exponent gradually decreases until, above a critical $\beta$ threshold, MSD reaches an asymptotic plateau, indicating complete diffusion suppression. This critical $\beta$ threshold appears to decrease with $\alpha$. All curves of all panels coincide at short times, as for $t<{\tau }_0=0.1$ resetting has not yet commenced, while for $t<1$ the walk is purely ballistic and $\alpha$ insensitive. At long times, the curves in each panel diverge demonstrating specific SR diffusion behavior.

Figure 5 graphically demonstrates the strong correlation between exponents $\alpha$ and $\beta$, which will be analyzed in detail in the forthcoming subsections. As anticipated in Section 2, we aim to determine how SR modifies the prefactor and exponent of the free MSD, $M_0\left(t\right)$ from Eq. (4). For this purpose, we must exclude from our analysis Lévy walks with $\alpha =1$, where Eq. (3) indicates that additional logarithmic time dependence must be incorporated. Cases with $\alpha =1$ will be treated separately in Section 4.4.

4.1. $0<\beta \le 1$: SR Independent Diffusion Exponent

As the mean resetting time in Eq. (8) diverges for $\beta <1$, we start our analysis from the Laplace transform of the full renewal equation in Eq. (18). Here, the survival function of the Pareto distribution of Eq. (20), reads

$$S\left(t\right)=\left\{\begin{array}{cc}{\int }_t^{\infty }\psi \left(\tau \right)d\tau ={\left({\displaystyle \frac{{\tau }_0}{t}}\right)}^{\beta },\hfill & t\ge {\tau }_0,\hfill \\ 1,\hfill & t<{\tau }_0.\hfill \end{array}\right.$$

(22)

Using Eq. (4) for $M_0\left(t\right)$ and neglecting contributions to $S\left(t\right)$ for $t<{\tau }_0$, we can approximate $S_0\left(t\right)=S\left(t\right)M_0\left(t\right)$ to $S_0\left(t\right)=K{\tau }_0^{\beta }{t}^{{\gamma }_0-\beta }$ for $\alpha \ne 1$ and calculate its Laplace transform,

$${\tilde{S}}_0\left(s\right)={\tau }_0^{\beta }\Gamma (1+{\gamma }_0-\beta )s^{\beta -(1+{\gamma }_0)}.$$

(23)

The Laplace transform of the Pareto distribution is special at $\beta =1$ (critical point) and needs to be calculated separately. Let’s calculate first the Laplace transform of the Pareto distribution for $0<\beta <1$ and $s\to 0$:

$$\tilde{\psi }\left(s\right)=1-\Gamma (1-\beta ){\tau }_0^{\beta }{s}^{\beta }+o\left(s^{\beta }\right).$$

(24)

Substituting ${\tilde{S}}_0\left(s\right)$, Eq. (23), and $\tilde{\psi }\left(s\right)$, Eq. (24), into Eq. (18) yields

$$\tilde{M}\left(s\right)=K{\displaystyle \frac{\Gamma (1+{\gamma }_0-\beta )}{\Gamma (1-\beta )}}{s}^{-(1+{\gamma }_0)}.$$

(25)

Finally, taking the inverse Laplace transform of $\tilde{M}\left(s\right)$ in Eq. (25), we obtain an analytical expression for the asymptotic Lévy walk MSD under power-law resetting of the form anticipated in Eq. (6), with

$$\gamma ={\gamma }_0,0<\beta <1,$$

(26)

and

$$K^{\prime }=K{\displaystyle \frac{\Gamma (1+{\gamma }_0-\beta )}{\Gamma (1-\beta )\Gamma (1+{\gamma }_0)}}.$$

(27)

We immediately observe that, while the diffusion exponent for $0<\beta <1$ remains invariant, the prefactor $K^{\prime }$ strongly depends on $\beta$.

Next, we consider the special case of $\beta =1$. For $s\to 0$, the Laplace transform of the Pareto distribution with $\beta =1$ reads

$$\tilde{\psi }\left(s\right)=1-s{\tau }_0\left[ln\left({\displaystyle \frac{1}{s{\tau }_0}}\right)+{\gamma }_E-1\right]+O(s^{2}lns).$$

(28)

Inserting Eqs. (23) and (28) into Eq. (18), we obtain

$$\tilde{M}\left(s\right)\sim K\Gamma \left({\gamma }_0\right){\displaystyle \frac{s^{-({\gamma }_0+1)}}{ln\left(1/s\right)}},$$

(29)

whence, upon inverse Laplace transformation,

$$M\left(t\right)\sim K{\displaystyle \frac{{\gamma }_0{t}^{{\gamma }_0}}{lnt}},\beta =1.$$

(30)

In Figure 6 we compare our analytical predictions for the relative prefactor change, $f=K^{\prime }/K$, with extensive numerical simulation results. The overall agreement is satisfactorily close. We observe that as $\beta \to 0$, resetting events become increasingly rare, and f approaches unity, indicating that SR becomes ineffective. Conversely, as $\beta$ increases from 0 to 1, f decreases monotonically to 0, suggesting suppression of unbounded diffusion. The $\beta$ dependence of f varies with the index $\alpha$; in particular, for $\alpha =2$ (Brownian walk), f decreases linearly with $\beta$, as first reported in Ref. [31].

4.2. $1<\beta <{\gamma }_0+1$: Diffusion Exponent Attenuation

When $\beta >1$, the mean value $\langle \tau \rangle$, Eq. (8), of the Pareto distribution is finite. This means that, as discussed in Section 3.1, the renewal equation (9) can be solved in terms of the age distribution to obtain the approximate solution of Eq. (12). The age distribution for the Pareto distribution can be readily obtained by inserting Eqs. (8) and (22) into Eq. (11), namely,

$$f_{\mathrm{age}}\left(t\right)=\left\{\begin{array}{cc}{\displaystyle \frac{1}{\langle \tau \rangle }}{\left({\displaystyle \frac{{\tau }_0}{t}}\right)}^{\beta },\hfill & \tau \ge {\tau }_0,\hfill \\ {\displaystyle \frac{1}{\langle \tau \rangle }},\hfill & \tau <{\tau }_0,\hfill \end{array}\right.$$

(31)

Upon substituting $M_0\left(t\right)$ from Eq. (4) and $f_{\mathrm{age}}\left(t\right)$ from Eq. (31) into Eq. (12) for $\alpha \ne 1$, $M\left(t\right)$ can be written as

$$M\left(t\right)={\displaystyle \frac{K}{\langle \tau \rangle }}{\int }_0^{{\tau }_0}{y}^{{\gamma }_0}dy+{\displaystyle \frac{K{\tau }_0^{\beta }}{\langle \tau \rangle }}{\int }_{{\tau }_0}^t{y}^{{\gamma }_0-\beta }dy.$$

(32)

In the limit $t\to \infty$, this integral expression must be treated differently based on the value of $\beta$. When $1<\beta <{\gamma }_0+1$, the second integral does not converge. Neglecting the first integral, the asymptotic expression for $M\left(t\right)$ can thus be formulated again according to Eq. (6) with

$$\gamma ={\gamma }_0+1-\beta ,1<\beta <{\gamma }_0+1,$$

(33)

and

$$K^{\prime }=K{\tau }_0^{\beta -1}{\displaystyle \frac{\beta -1}{\beta ({\gamma }_0-\beta +1)}}.$$

(34)

In contrast, for $\beta ={\gamma }_0+1$, the second integral in Eq. (32) diverges logarithmically, so that

$$M\left(t\right)\sim K{\displaystyle \frac{\beta -1}{\beta }}{\tau }_0^{\beta -1}lnt,\beta ={\gamma }_0+1.$$

(35)

Unlike Section 4.1, in this $\beta$ regime, SR attenuates the diffusion exponent from ${\gamma }_0$ to $\gamma ={\gamma }_0+1-\beta$. In Eq. (6), the limit $\gamma \to 0+$ when $\beta \to {\gamma }_0+1$ from below, disguises a logarithmic divergence of $M\left(t\right)$. As demonstrated in the following subsection, when the same limit of $\beta$ is approached from above, $\gamma \to 0-$ indeed implies localization.

Figure 7 illustrates the dependence of diffusion exponent $\gamma$ on Pareto exponent $\beta$ for various $\alpha$ values. For $\beta <1$, all curves are horizontal, indicating that the diffusion exponent remains constant. When $1<\beta <{\gamma }_0+1$, the diffusion exponent decreases linearly, as predicted by Eq. (33). The transition from superdiffusion to localization is clearly evident as $\beta$ increases beyond ${\gamma }_0+1$, with MSD from numerical simulations saturating to a constant plateau for $\beta >{\gamma }_0+1$.

4.3. $\beta >{\gamma }_0+1$: Localization

When $\beta >{\gamma }_0+1$, both integrals in Eq. (32) converge as $t\to \infty$, so that the Lévy MSD under SR now approaches a plateau. Direct integration of Eq. (32) yields

$${\mathrm{M}}_{\mathrm{st}}=K{\tau }_0^{{\gamma }_0}{\displaystyle \frac{\beta -1}{({\gamma }_0+1)(\beta -{\gamma }_0-1)}},\beta >{\gamma }_0+1.$$

(36)

Figure 8 compares diffusion plateaus ${\mathrm{MSD}}_{\mathrm{st}}$ in the localization regime, $\beta >{\gamma }_0+1$, with the analytical estimate ${\mathrm{M}}_{\mathrm{st}}$ from Eq. (36). The agreement is excellent across the entire $\alpha$ range. It can be observed that the plateau appears to diverge as $\beta \to {\gamma }_0+1$ from above, which is consistent with the logarithmic divergence of $M\left(t\right)$, Eq. (35), when the same limit is approached from below. The inset shows ${\mathrm{M}}_{\mathrm{st}}$ for large $\beta$ and ${\tau }_0=0.1$. As $\beta \to \infty$, $\langle \tau \rangle \to {\tau }_0$ [Eq. (8], which implies a high resetting frequency with $\tau$ close to ${\tau }_0$. Consequently, the plateau onset time is also of order ${\tau }_0$. With $c=v=1$ in the Lévy step generator [Eq. (1)], frequent resetting events with $\tau <1$ confine the walker to the short-time ballistic diffusive regime $M\left(t\right)=t^2$ (independent of $\alpha$) described in Section 2. This technical limitation can be circumvented by setting ${\tau }_0\gg 1$; in Fig. Figure 8, consistently with the discussion of Fig. Figure 4, we set ${\tau }_0=100$. However, it is possible to validate Eq. (36) even for ${\tau }_0<1$. Indeed, adopting ballistic parameter values, $K=1$ and ${\gamma }_0=2$, in Eq. (4), Eq. (36) predicts ${\mathrm{M}}_{\mathrm{st}}$ = 0.0033, in close agreement with the simulation result 0.0035.

4.4. The Case $\alpha =1$

For $\alpha =1$, $M_0\left(t\right)$ exhibits logarithmic divergence as shown in Eq. (3). This implies that asymptotic power-law scaling is untenable for the Lévy walk MSD both in the absence [Eq. (4)] and presence [Eq. (6)] of SR. For this reason Lévy walks with $\alpha =1$ must be treated separately. Simulation data for $0<\beta <1$ in Figs. Figure 5(b) and Figure 6 suggest the convergence of the t function $M\left(t\right)/M_0\left(t\right)$ toward a finite ratio, f, also for $\alpha =1$. To explain this result, we repeated the analytical procedure of Eqs. (23) - (27) for $M_0\left(t\right)=t^2/lnt$ instead of $M_0\left(t\right)=K{t}^{{\gamma }_0}$. In the limit of small s, one obtains

$$\tilde{M}\left(s\right)={\displaystyle \frac{{\tilde{S}}_0\left(s\right)}{1-\tilde{\psi }\left(s\right)}}\approx {\displaystyle \frac{\Gamma (3-\beta )}{\Gamma (1-\beta )}}·{\displaystyle \frac{s^{-3}}{-lns}},$$

(37)

whose inverse Laplace transform,

$$M\left(t\right)\sim {\displaystyle \frac{\Gamma (3-\beta )}{2\Gamma (1-\beta )}}{\displaystyle \frac{t^2}{lnt}},\alpha =1,0<\beta <1,$$

(38)

scales with time exactly like $M_0\left(t\right)$. The relative change in the prefactor, $f=\Gamma (3-\beta )/\left[2\Gamma (1-\beta )\right]$ coincides with that obtained from Eq. (27) for $0<\alpha <1$, i.e., for ${\gamma }_0=2$, as shown in Fig. Figure 6.

Stochastic resetting diffusion attenuation for $\alpha =1$ is expected to occur in the Pareto exponent range $1<\beta <3$ - we remind here that ${\gamma }_0\to 2$ as $\alpha \to 1$ from both below and above. Adopting the procedure of Section 4.2 for $M_0\left(t\right)=t^2/lnt$ leads to replacing Eq. (32) with

$${\mathrm{M}}_{\mathrm{st}}={\displaystyle \frac{1}{\langle \tau \rangle }}{\int }_0^{{\tau }_0}{\displaystyle \frac{y^2}{lny}}dy+{\displaystyle \frac{{\tau }_0^{\beta }}{\langle \tau \rangle }}{\int }_{{\tau }_0}^t{\displaystyle \frac{y^{2-\beta }}{lny}}dy.$$

(39)

On following the integration procedure adopted in Section 3.1, one recognizes immediately that the asymptotic regime of $M\left(t\right)$ is governed by the diverging second integral, that is,

$$M\left(t\right)\approx {\tau }_0^{\beta -1}{\displaystyle \frac{\beta -1}{\beta (3-\beta )}}{\displaystyle \frac{t^{3-\beta }}{lnt}},\alpha =1,1<\beta <3.$$

(40)

Comparing this result with the input free MSD, $M_0\left(t\right)=t^2/lnt$, indicates that the t-exponent and the prefactor in Eq. (40) coincide with $\gamma$ and $K^{\prime }$, respectively, from Eqs. (33) and (34), for $K=1$ and ${\gamma }_0=2$.

Finally we consider SR localization in the range $\beta >3$. An estimate for the plateaus, MSD${\mathrm{MSD}}_{\mathrm{st}}$, is also analytically tractable. One simply needs to take the limit $t\to \infty$ of Eq. (39) to obtain

$${\mathrm{M}}_{\mathrm{st}}\simeq {\displaystyle \frac{{\tau }_0^2(\beta -1)}{3(\beta -3)|ln{\tau }_0|}},\alpha =1,\beta >3.$$

(41)

This estimate for ${\mathrm{M}}_{\mathrm{st}}$ is also in good agreement with the simulation, as demonstrated in Fig. Figure 8.

The special cases $\beta =1$ [Eq. (38)] and $\beta =3$ [Eq. (40)] for $\alpha =1$ required solving Eqs. (18) and (12), respectively, for the specific memory kernel $M_0\left(t\right)=t^2/lnt$. Through laborious analytical calculations, we obtained approximate expressions for $M\left(t\right)$ in the large-t regime,

$$M\left(t\right)\simeq {\displaystyle \frac{1}{2}}{\displaystyle \frac{t^2}{{(lnt)}^2}},\alpha =1,\beta =1,$$

(42)

and

$$M\left(t\right)\simeq {\displaystyle \frac{2{\tau }_0^2}{3}}ln(lnt),\alpha =1,\beta =3,$$

(43)

which are consistent with the data from numerical simulations displayed in Fig. Figure 5(b).

We conclude this subsection remarking that, as in the case of $\alpha \ne 1$, SR tends to weaken the superdiffusive behavior of the Lévy walk for $\beta \ge 1$, while for $0<\beta <1$ it only modifies the prefactor.

5. Conclusions

In this paper we investigated the diffusive dynamics of Lévy walks under stochastic resetting with Pareto distributions of the resetting times. Analytical expressions for the asymptotic mean squared displacement (MSD) of Lévy walks under power-law resetting are derived using renewal theory, and validated through numerical simulations.

Our main results are summarized in the 2D phase diagram of Figure 9. There we display the dependence of the diffusion exponent $\gamma$ as a function of the Lévy index $\alpha$, and resetting-time distribution exponent $\beta$. The horizontal axis $\alpha$, ranging from 0 to 2, is divided into two segments: for $0<\alpha <1$, the MSD of the free Lévy walk ($\beta =0$) exhibits ballistic growth, i.e., $\langle x^2\rangle \sim t^2$; for $1<\alpha <2$, the free MSD behaves as $\langle x^2\rangle \sim t^{3-\alpha }$; finally, $\alpha =2$, corresponds to the regular Brownian motion with $\langle x^2\rangle \sim t$. This defines the free diffusion exponent ${\gamma }_0$. On the vertical axis, the exponent $\beta$ characterizes resetting event frequency: smaller $\beta$ indicates fatter tails in the resetting time distribution and less frequent resetting; larger $\beta$ indicates thinner tails and more frequent resetting. As shown by the figure color code, in the $\beta <1$ region, regardless of $\alpha$ value, $\gamma ={\gamma }_0$. As $\beta$ enters the $1<\beta <{\gamma }_0+1$ range, the diffusion exponent decreases linearly as $\gamma ={\gamma }_0+1-\beta$, until it reaches zero at $\beta ={\gamma }_0+1$. When $\beta >{\gamma }_0+1$, the diffusion exponent vanishes, corresponding to the MSD converging to a constant plateau (or logarithmic divergence for singular values $\alpha =1$ and $\beta =1$).

The results presented herein demonstrate robustness against variations of the fundamental model outlined in Section 2. For instance, following Refs. [17,32], we examined the more realistic scenario of non-instantaneous stochastic resetting with ballistic return dynamics characterized by finite large speed. We found that short finite return times do not significantly modify the diffusion exponents $\gamma$, but slightly enhance the diffusion prefactors, $K^{\prime }$.

Beyond their fundamental significance in stochastic process theory, studies of Lévy walk diffusion under stochastic resetting provide a foundation for optimizing biological foraging strategies, analyzing population dynamics, and enhancing financial risk assessment methodologies.