L\'evy flights versus L\'evy walks in bounded domains

TL;DR

This paper compares Le9vy flights and Le9vy walks in bounded domains, analyzing their differences and similarities in terms of statistical properties like survival probability and mean first passage time, influenced by boundary conditions and jump distribution.

Contribution

It provides a detailed analytical and numerical comparison of Le9vy flights and Le9vy walks within bounded domains, highlighting conditions for their similarity.

Findings

Similarity depends on boundary conditions and stability index.

Le9vy flights exhibit pathological properties addressed by Le9vy walks.

Models show comparable results under certain conditions.

Abstract

L\'evy flights and L\'evy walks serve as two paradigms of random walks resembling common features but also bearing fundamental differences. One of the main dissimilarities are discontinuity versus continuity of their trajectories and infinite versus finite propagation velocity. In consequence, well developed theory of L\'evy flights is associated with their pathological physical properties, which in turn are resolved by the concept of L\'evy walks. Here, we explore L\'evy flights and L\'evy walks models on bounded domains examining their differences and analogies. We investigate analytically and numerically whether and under which conditions both approaches yield similar results in terms of selected statistical observables characterizing the motion: the survival probability, mean first passage time and stationary PDFs. It is demonstrated that similarity of models is affected by the type of boundary conditions and value of the stability index defining asymptotics of the jump length distribution.

Figures

Tables

Table 1. Table 1: Values of K α LW superscript subscript 𝐾 𝛼 LW K_{\alpha}^{\mathrm{LW}} calculated according to Eq. ( 55 ).

$\alpha$	1.1	1.3	1.5	1.7	1.9
$K_{\alpha }^{\mathrm{LW}}$	0.1495	0.3981	0.5863	0.7294	0.8401

Equations

$$\frac{dx}{dt}={\zeta }_{\alpha }(t),$$

$$\Delta x=x(t+\Delta t)-x(t)=\Delta t^{1/\alpha }{\zeta }_t$$

$$\varphi (k)=\exp \left[-{\sigma }_0^{\alpha }|k{|}^{\alpha }\right],$$

$$p_{\alpha }(\Delta x;{\sigma }_0)\sim \frac{{\sigma }_0^{\alpha }\Delta t\sin \frac{\pi \alpha }{2}\Gamma (\alpha +1)/\pi }{|\Delta x{|}^{\alpha +1}}.$$

$$\sigma (t)={\sigma }_0{t}^{1/\alpha },$$

$$\frac{\partial P(x,t|x_0,0)}{\partial t}={\sigma }_0^{\alpha }\frac{{\partial }^{\alpha }P(x,t|x_0,0)}{\partial |x{|}^{\alpha }}=K_{\alpha }\frac{{\partial }^{\alpha }P(x,t|x_0,0)}{\partial |x{|}^{\alpha }},$$

$$\tau (x_0)=\min \{t>0:x(0)=x_0\text{mbox}and|x(t)|⩾L\}.$$

$$⟨\tau (x_0)⟩=\frac{1}{\Gamma (1+\alpha )}\frac{(L^2-|x_0{|}^2{)}^{\alpha /2}}{{\sigma }_0^{\alpha }}$$

$$⟨\tau ⟩=\frac{\pi }{2^{1+\alpha }\Gamma \left[(1+\alpha )/2\right]\Gamma \left[(3+\alpha )/2\right]}\frac{L^{\alpha }}{{\sigma }_0^{\alpha }},$$

$$S(t|x_0)=1-\mathcal{F}(t|x_0)=1-{\int }_0^{t}p(s|x_0)ds,$$

$$S(t|x_0)=\sum _{i=1}^{\infty }{c}_i(x_0)\exp \left[-{\lambda }_i^{(\alpha )}t\right],$$

$$S(t|x_0)\approx \exp \left[-{\lambda }_1^{(\alpha )}t\right].$$

$${\lambda }_n^{(\alpha )}={\left[\frac{n\pi }{2}+\frac{(2-\alpha )\pi }{8}\right]}^{\alpha }.$$

$$S(t|x_0)\approx \exp \left[-t/⟨\tau (x_0)⟩\right],$$

$$⟨\tau (0)⟩=\frac{L^2}{2{\sigma }_0^2}.$$

$${\lambda }_1^{(2)}=\frac{{\pi }^2{\sigma }_0^2}{4L^2},$$

$$\tau (x_0)=\min \{t>0:x(0)=x_0\text{mbox}andx(t)⩾L\}.$$

$$\underset{n\to \infty }{\lim }{V}_n(x)=\underset{n\to \infty }{\lim }\frac{1}{n}\frac{|x{|}^n}{L^n},$$

$$\frac{dx}{dt}=-V_n^{\prime }(x)+{\zeta }_{\alpha }(t).$$

$$P_{st}(x)=\frac{\Gamma (\alpha )(2L{)}^{1-\alpha }(L^2-x^2{)}^{\alpha /2-1}}{{\Gamma }^2(\alpha /2)}.$$

$$\mathrm{Pr}\left(t,x_0\right)={\int }_{-L}^{L}P\left(x,t|x_0,0\right)dx.$$

$$P\left(x,0\right)=\delta \left(x-x_0\right).$$

$$\vartheta \left(k,t\right)={\int }_{-\infty }^{\infty }P\left(x,t\right)e^{ikx}dx$$

$$\frac{\partial \vartheta }{\partial t}=-{\sigma }_0^{\alpha }{\left|k\right|}^{\alpha }\vartheta .$$

$$\vartheta \left(k,0\right)=e^{ik{x}_0}.$$

$$\vartheta \left(k,t\right)=\exp \left[ik{x}_0-{\sigma }_0^{\alpha }{\left|k\right|}^{\alpha }t\right],$$

$$P\left(x,t|x_0,0\right)=\frac{1}{2\pi }{\int }_{-\infty }^{\infty }{e}^{-ik\left(x-x_0\right)-{\sigma }_0^{\alpha }{|k|}^{\alpha }t}dk.$$

$$\mathrm{Pr}\left(t,x_0\right)=\frac{2}{\pi }{\int }_0^{\infty }\frac{\sin kL\cos k{x}_0}{k}{e}^{-{\sigma }_0^{\alpha }{k}^{\alpha }t}dk.$$

$$⟨T\left(x_0\right)⟩={\int }_0^{\infty }\mathrm{Pr}\left(t,x_0\right)dt.$$

$$⟨T\left(x_0\right)⟩=\frac{1}{\pi {\sigma }_0^{\alpha }}{\int }_0^{\infty }\frac{\sin k\left(L+x_0\right)+\sin k\left(L-x_0\right)}{k^{1+\alpha }}dk.$$

Full text

I Introduction

Lévy flights Shlesinger et al. (1995) and Lévy walks Klafter et al. (1996); Zaburdaev et al. (2015) are two well known stochastic models of anomalous diffusion. Generally speaking, Lévy flights correspond to Markovian motions whose individual, uncorrelated random steps are drawn from a Lévy distribution, thus extending Brownian motion for which the step lengths are Gaussian. A resulting asymptotic Lévy diffusion is then characterized by infinite variance, indicating that the width of the diffusive “packet” must be understood in terms of some fractional moments or the interquartile distance Dybiec et al. (2008). This mathematical property of Lévy flights along with their instantaneous propagation are considered in many situations unphysical. In contrast, in Lévy walks Shlesinger and Klafter (1986); Klafter et al. (1996); Zaburdaev et al. (2015), the meandering particle has a finite velocity, so that long jumps take proportionally longer time. Still, in the absence of any boundary effect, the core of the Lévy walk packet disperses faster than linearly in time but slower than the ballistic front and is described by the Lévy distribution. This means that under free boundary conditions Lévy flights can serve as a good approximation to the Lévy walk, although with an improper prediction of the moments of the jump length distribution Zaburdaev et al. (2015).

In this paper we analyze these two popular models of stochastic motion in bounded domains. Such motions can e.g. represent foraging behaviors of animals and bacteria Zaburdaev et al. (2015); Metzler et al. (2012), the spreading of diseases Hufnagel et al. (2004) or particle transport along soft polymer chains Sokolov et al. (1997). The problem we aim to address is whether both aforementioned approaches yield similar results in terms of investigated kinetics (survival, occupation and first passage times) and long-term behavior (existence of stationary states and stationary probability densities). Depending on the value of the stability index $\alpha$ and used characteristics we find similarities but also large deviations between the two models. We furthermore compare explicit analytical results with numerical simulations of stochastic dynamics. Direct comparison of two models – motion of random walkers flying instantaneously between distinct sites and walkers performing motion at a constant speed – deepens our understanding of their behavior and relates characteristic properties amenable to measurements in real situations.

The article is organized as follows: Sec. II discuses the problem of boundary conditions, mean first passage time, mean residence time and stationary states for Lévy flights. In Sec. III problems of boundary conditions, mean first passage time and stationary states for Lévy walks are explored. The manuscript is closed with summary and discussion.

II Lévy flights in one dimensional intervals

Let us briefly reconsider the motion of a free overdamped particle described by the Langevin equation:

[TABLE]

where $\zeta_{\alpha}(t)$ is a symmetric white $\alpha$-stable noise, i.e. the formal time derivative of the symmetric $\alpha$-stable motion Janicki and Weron (1994a). Note, that in the Lévy flight (LF) scenario we do not take into consideration inertial effects, and similar to the Wiener process we are dealing with an overdamped kind of motion. In contrast, in the Lévy walk (LW) scheme, soon to be considered, we do include inertial effects and finite propagation velocity. Hence one origin of the difference between both models stems from neglecting the inertia.

Eq. (1) is supplemented with the initial condition $x(0)=x_{0}$. Accordingly, the stochastic process $\{X(t),t\geqslant 0\}$ governed by Eq. (1) has increments

[TABLE]

distributed according to the symmetric $\alpha$-stable density with the scale parameter depending on the discretization time step $\Delta t$. The discretization procedure assures that for a free particle with an arbitrary $\Delta t$, time dependent densities do not depend on the discretization time step but on $t$ (and remaining parameters) only. In Eq. (2) $\zeta_{t}$ represents independent, identically distributed (i.i.d) random variables following the symmetric $\alpha$-stable density Janicki and Weron (1994b); Janicki (1996) with the characteristic function $\phi(k)$

[TABLE]

where $\sigma_{0}>0$ is the scale parameter. The stability index $\alpha$ describes asymptotic (large $\Delta x$) behavior of the jump length density:

[TABLE]

Note that the parameter $\sigma_{0}$ scales the overall distribution width, hence its role is similar to a standard deviation for distributions with a finite second moment. From Eq. (1) and arithmetic properties of $\alpha$-stable densities Janicki (1996); Samorodnitsky and Taqqu (1994) it follows that the process $\{X(t),t\geqslant 0\}$ is distributed according to the $\alpha$-stable density with the time dependent parameter $\sigma(t)$

[TABLE]

where $\sigma_{0}$ is a fixed scale parameter associated with the underlying $\alpha$-stable white noise in Eq. (1). Consequently, its asymptotics is described by Eq. (4) after substitution of $\Delta t$ with $t$ and $\Delta x$ with $x$.

The Langevin equation (1) can be associated with the space-fractional Smoluchowski-Fokker-Planck equation

[TABLE]

which describes evolution of the probability density function (PDF) of finding a random walker in the vicinity of $x$ at time $t$ with the initial condition $P(x,0|x_{0},0)=\delta(x-x_{0})$. The fractional operator $\frac{\partial^{\alpha}}{\partial|x|^{\alpha}}$ stands for the fractional Riesz-Weil derivative, defined by the Fourier transform Podlubny (1999); Samko et al. (1993) $\mathcal{F}_{k}(\frac{\partial^{\alpha}f(x)}{\partial|x|^{\alpha}})=-|k|^{\alpha}\mathcal{F}_{k}(f(x)).$ In what follows, we interpret $K_{\alpha}=\sigma_{0}^{\alpha}$ in Eq. (6) as the generalized diffusion constant.

In the next subsections the main scope is to investigate interrelationship between formulation of boundary conditions for Lévy flights and properties of escape kinetics and stationary states. More precisely, in order, to assess various formulation of boundary conditions we explore two scenarios of escape kinetics from finite intervals (a) restricted by two absorbing boundaries and (b) restricted by reflecting and absorbing boundaries. We compare exact results (when applicable) with numerically estimated mean first passage times. Moreover, for a finite interval restricted by two reflecting boundaries we verify if numerically constructed stationary densities agree with analytical predictions. Finally, we study the various finite interval setups not only for Lévy flights but also for Lévy walks, which contrary to Lévy flights have continuous trajectories and finite propagation velocity. Comparison between behavior of the two models (LF and LW), with respect to a class of important observables, define the main scope of current research. Depending on the value of the stability index $\alpha$ and observable type we find similarities but also large deviations between the two models.

II.1 First escape problem

In the presence of boundary conditions imposed for Eq. (6) the translational invariance is broken and the resulting evolution equation for PDFs becomes a non-trivial integro-differential equation with non-local boundary conditions Dybiec et al. (2006); Zoia et al. (2007). In order to avoid the problem, analysis of first escape is studied here by use of the Langevin methods, for which – contrary to the methods based on the fractional Smoluchowski-Fokker-Planck equation – the implementation of boundary conditions is significantly simpler, although not fully resolved.

II.1.1 Absorbing boundaries at both ends

We consider a first escape problem from the $[-L,L]$ interval with both boundaries being absorbing, see Fig. 1. The evolution of $x(t)$ is determined by the Langevin Eq. (1) and the first passage time $\tau(x_{0})$ ($|x_{0}|\leqslant L$) is defined as

[TABLE]

In this case, the formula for the mean first passage time (MFPT, $\langle\tau(x_{0})\rangle$) reads Getoor (1961); Zoia et al. (2007)

[TABLE]

demonstrating that the MFPT asymptotically scales as $(L/\sigma_{0})^{\alpha}$, what is especially well visible for $x_{0}=0$ when such a scaling is recorded for the whole $L$ range. Note that, similar formula can be also found in Buldyrev et al. (2001). Eq. (8) can be averaged over initial conditions. For example assuming that $x_{0}$ is uniformly distributed over $[-L,L]$ the mean exit time reads

[TABLE]

which has exactly the same $(L/\sigma_{0})^{\alpha}$ dependence like Eq. (8) for $x_{0}=0$. Therefore, in forthcoming considerations the fixed initial $x_{0}=0$ condition is used.

In numerical simulations of the corresponding Langevin equation the absorption condition is realized by assuming that the whole exterior of the prescribed interval is absorbing, i.e. each time when a trajectory crosses the absorbing boundary, it is removed and the first passage time is recorded. Results of numerical simulations and formula (8) are presented in Fig. 1 showing perfect agreement with the theoretical curve.

In Fig. 1, the scale parameter and the interval half-width have been arbitrary preset to $\sigma_{0}=1$, $L=1$ and the initial condition to $x_{0}=0$. Consequently, the $\langle\tau(x_{0})\rangle(\alpha)$ curve attains one of its possible forms. In more general cases of $\sigma_{0}\neq 1$ and $L\neq 1$, or more precisely the ratio $L/\sigma_{0}\neq 1$, this curve can be of very different type, see Szczepaniec and Dybiec (2015).

The cumulative distribution of first passage times for this problem defines the survival probability

[TABLE]

derived with the corresponding PDF $p(t|x_{0})$ and depicted in Fig. 2. Clearly, the survival probability denotes the probability that a process starting at $x(0)=x_{0}=0$ has not reached or crossed up to time $t$ the levels $\pm L$. Note that, by construction, the process described by Eq. (1) is Markovian, which remains in line with observation of exponential asymptotics in Fig. 2. The behavior is well documented in simulations of Lévy flights Dybiec et al. (2006); Dybiec and Sokolov (2015) and can be inferred by estimation of lower and upper bounds Katzav and Adda-Bedia (2008); Chen et al. (2010); Kim et al. (2015) for tails of $S(t|x_{0})$ or from the master equation Drysdale and Robinson (1998); Buldyrev et al. (2001); Araújo and Raposo (2016). It can be also deduced by separation of variables Zoia et al. (2007); Dybiec and Sokolov (2015)

[TABLE]

with $\lambda_{i}^{(\alpha)}$ denoting eigenvalues of the fractional Laplacian on bounded domains Kwaśnicki (2012). Accordingly, the decay of the survival probability $S(t|x_{0})$ is determined by the smallest eigenvalue of the fractional Laplacian Zoia et al. (2007); Katzav and Adda-Bedia (2008); Kwaśnicki (2012); Dybiec and Sokolov (2015), prompting a long time approximation

[TABLE]

The smallest eigenvalue $\lambda_{1}^{(\alpha)}$ can be calculated according to (Kwaśnicki, 2012, Eq. (2)) which for $L=1$ reads

[TABLE]

Eq. (12) along with properties of the survival probability, i.e. $\langle\tau(x_{0})\rangle=\int_{0}^{\infty}S(t|x_{0})dt$, suggest another possible approximation to the survival probability

[TABLE]

where $\langle\tau(x_{0})\rangle$ is the mean first passage time given by Eq. (8).

Figure 2 compares both approximations, see Eqs. (12) and (14).

Solid lines represent Eq. (14): Lévy motion on a confined interval between two absorbing boundaries decays with the steepness parameter depending on the stability index $\alpha$. At the same time, with increasing $\alpha$ the deviations from a single exponential approximation given by Eq. (14) become more pronounced, as more and more terms from Eq. (11) have to be retained Katzav and Adda-Bedia (2008) in order to reconstruct an initial part of the survival probability. Therefore, approximation (14) does not reproduce correct decay rate of the survival probability.

Additional dotted lines in Fig. 2 depict single exponential, smallest eigenvalue approximation (12), which does not work perfectly, but it predicts the right exponent characterizing the asymptotic slope of the survival probability.

In particular, due to pedagogical reasons, approximations (12) and (14) can be compared for $\alpha=2$. In such a case, the MFPT can be calculated from Eq. (8)

[TABLE]

The smallest eigenvalue of the Laplacian is Cox and Miller (1965); Dybiec and Gudowska-Nowak (2012)

[TABLE]

leading to $\langle\tau(x_{0})\rangle=1/\lambda_{1}^{(2)}=4L^{2}/\pi^{2}\sigma_{0}^{2}\approx 0.405L^{2}/\sigma_{0}$, which differs by $24\%$ from the exact value, see Eq. (15). The quality of approximation (14) depends on the exact value of the stability index $\alpha$ what can be inferred from Fig. 3, which presents the ratio of the exact value of the MFPT, see Eq. (8), to its approximation $\langle\tau(x_{0})\rangle\approx 1/\lambda_{1}^{(\alpha)}$. The dependence of $\langle\tau(x_{0})\rangle\times\lambda_{1}^{(\alpha)}$ is a non monotonous function of the stability index $\alpha$ with the maximum located around $\alpha\approx 1$. For small $\alpha$ Eqs. (12) and (14) provide reasonable approximations to the survival probability. Finally, it is worthy to underline, that quality of the approximation (14) becomes worse with the increasing number of spatial dimensions Dybiec and Szczepaniec (2015).

II.1.2 Reflecting-absorbing boundary conditions

Next, we consider a first escape problem from the $[-L,L]$ interval with reflecting (left) and absorbing (right) boundaries, see Fig. 4. The first passage time $\tau(x_{0})$ ($|x_{0}|\leqslant L$) is then defined as

[TABLE]

Analogously, like in Sec. II.1.1 we use $x_{0}=0$.

Imposing a reflecting boundary at $x=-L$ requires some additional care in numerical simulations. Here we consider three different realizations of the reflecting condition:

• (i)

motion reversal: a trajectory which ends at $x<-L$ is wrapped around the left boundary, i.e. $x\to-L+|x+L|$;

• (ii)

motion stopping: a trajectory which crosses $-L$ is paused at $-L+\varepsilon$, where $\varepsilon$ is a small and positive parameter. The point $-L+\varepsilon$ is used as starting point for a next jump;

• (iii)

motion confined within a potential: the reflecting boundary can be implemented by considering the motion in a bounding potential

[TABLE]

which is described by the following Langevin equation

[TABLE]

For $n\geqslant 2$, the potential (18) and its first derivative (force) are continuous. In the limit of $n\to\infty$, the potential $V_{n}(x)$ reduces to the infinite steep rectangular potential well with (reflecting) boundaries at $\pm L$.

In all above scenarios we use Langevin equations (1) or (19) to simulate the Lévy motion. The absorbing boundary is executed in a standard way, i.e. whenever $x\geqslant L$, a particle becomes absorbed, both for a free particle motion and motion in the potential well (18).

As it can be inferred from the middle panel of Fig. 4, the stopping and potential scenarios yield same values of estimated MFPT, conditioned that $n$ and $\varepsilon$ are sufficiently large and small, respectively. In contrast, the implementation of the wrapping method underestimates MFPT in comparison to the other two cases. All scenarios are equivalent for $\alpha=2$ indicating that the source of discrepancy lies in discontinuity of trajectories for Lévy motion with $\alpha<2$. In fact, this conclusion can be drawn by a more accurate analysis of wrapping trajectories: the bottom panel of Fig. 4 presents the fraction of escape events (wrapping scenario) in which a particle staring at $x_{0}=0$ escaped from the $[-1,1]$ interval by a single long jump to the left and has been reversed around the reflecting boundary. This fraction decreases with increase of the stability index $\alpha$ and tends to be arbitrarily small in the Gaussian limit ($\alpha\to 2$), when the trajectory $x(t)$ becomes continuous. In the opposite limit of $\alpha\to 0$ almost half of escape events are due to trajectory wrapping, as might be expected for extremely wide and heavy tailed PDF of increments $\Delta x$.

Additionally to the analysis of escape events from the $[-L,L]$ interval (reflecting-absorbing), the escape from the $[-2L,2L]$ interval (absorbing-absorbing) has been investigated. Such a scenario gives the same results as wrapping, since the escape from the reflecting-absorbing interval is equivalent to the escape from two times wider absorbing-absorbing interval, see central panel of Fig. 4 where $\langle\tau(0)\rangle$ given by Eq. (8) for a system $[-2L,2L]$ is depicted. The theoretical formula perfectly matches the simulations on the half sized system with the reversing (wrapping) condition. For $\alpha=2$ and various types of boundary conditions, this correspondence can be trivially read off from the analytic formula for the MFPT Gardiner (2009). In more general settings with $\alpha<2$, the equivalence relies on infinite propagation of the trajectory $x(t)$: from any point $x$ the distance to the absorbing boundary is either $L-x$ or $x+3L$ (when wrapped along a reflecting boundary). The sum of these two distances is the only relevant model parameter and equals exactly the sum of distances to absorbing boundaries located at $\pm 2L$.

Figure 5 presents survival probabilities $S(t|x_{0})$ corresponding to motions analyzed in Fig. 4. The top panel presents results for motion in the potential $V_{800}(x)$ and in a rectangular potential well with stopping. These two scenarios result in the same first passage time PDFs and, accordingly, in the same mean first passage times. The bottom panel of Fig. 5 exemplifies results for the survival probability in the wrapping scenario which leads to faster escape kinetics. Analogously to the absorbing-absorbing setup, the first passage time densities display exponential asymptotics reflecting Markovian characteristics of the process $x(t)$, cf. Eqs. (1) and (19).

Solid lines in Fig. 5 represent $S(t|x_{0})\approx\exp\left[-t/\langle\tau(x_{0})\rangle\right]$ approximation to the survival probability. Analogously, like in Fig. 2 this kind of approximation works better for small values of the stability index $\alpha$. Therefore, the largest discrepancies are observed for $\alpha=2$. For the wrapping scenario discrepancies are smaller than for remaining scenarios, because wrapping of trajectories accelerates decay of the survival probability.

II.2 Stationary states

In case of Lévy flights in an infinitely deep rectangular well, a particle executing the motion becomes confined in a domain restricted by two impenetrable boundaries (cf. top panel of Fig. 6). With these conditions, stationary states can be observed. In what follows, we compare simulations of Eq. (19) in the limit of large $n$ with the results of simulated free diffusion bounded by the stopping or wrapping conditions at reflecting boundaries at $x=\pm L$. The stopping (empty points) and potential (filled points) scenarios produce stationary PDFs displayed in Fig. 6, in full agreement with analytical result derived by Denisov et al. Denisov et al. (2008):

[TABLE]

Bottom panel of Fig. 6 depicts cumulative densities $\mathcal{F}_{x}(x)=\int_{-L}^{x}P_{st}(x^{\prime})dx^{\prime}$ corresponding to histograms presented in Fig. 6 and calculated (solid lines) from Eq. (20). For clarity of the presentation, the cumulative densities have been plotted for $0\leqslant x\leqslant L$ with $L=1$ only.

II.3 Mean Residence time

It is instructive to compare results of the above MFPT analysis with the mean residence time (MRT). MRT represents the average time that a freely diffusing particle, moving on $-\infty<x(t)<\infty$ resides in a given region (say, in the interval $[-L,L]$) in a measurement process of duration $t$. Consequently, the MRT is always shorter than the measurement time $t$. Moreover, due to possible multiple returns to the $[-L,L]$ interval MRT can be larger than MFPT. We will consider two limits, when measurement time is infinite, and also the time dependence of the MRT. We start with evaluating the probability to find a particle in the interval $[-L,L]$ at the time $t$, provided it has started from some internal point $x_{0}\in\left(-L,L\right)$

[TABLE]

The probability density of transitions $P\left(\left.x,t\right|x_{0},0\right)$ in Eq. (21) is the solution of Eq. (6) with the initial condition

[TABLE]

If we make the Fourier transform of Eq. (6) and introduce the characteristic function which is defined as

[TABLE]

we arrive at

[TABLE]

From Eqs. (22) and (23) we find the initial condition for Eq. (24)

[TABLE]

The exact solution of Eq. (24) with the initial condition (25) reads

[TABLE]

which is the $\alpha$-stable density with the same stability index $\alpha$ like the underlying noise and the time dependent scale parameter $\sigma_{0}t^{1/\alpha}$, compare Eqs. (26) and (3). A non-zero initial condition introduce additional shift to the time dependent density, which is centered at $x_{0}$.

Using the inverse Fourier transform we obtain from Eq. (26)

[TABLE]

Substituting Eq. (27) in Eq. (21) and integrating over $x$ we arrive at

[TABLE]

According to the definition, MRT $\langle T(x_{0})\rangle$ in the interval $[-L,L]$ can be calculated as

[TABLE]

This MRT is the mean residence time in the limit where the total measurement time is large, below we will obtain the time dependence of the MRT. Substituting Eq. (28) in Eq. (29) and integrating over $t$ we get

[TABLE]

As seen from Eq. (30), the integral diverges for $1\leqslant\alpha<2$, i.e. MRT goes to infinity as for standard Brownian diffusion ($\alpha=2$). Divergence of the mean residence time is the consequence of multiple returns to the $[-L,L]$ interval. In the case $0<\alpha<1$ Eq. (30) yields

[TABLE]

where $\Gamma(x)$ is the Euler gamma function. In particular, for $x_{0}=0$ from Eq. (31) we have

[TABLE]

MRT given by Eqs. (31) and (32) is finite, but larger than the MFPT. In the limit of $\alpha\to 0$, the mean residence time in Eq. (32) tends to $1$, which is the mean first passage time from the $[-L,L]$ interval given by Eq. (8). Equality of MRT and MFPT signals that the particle after escaping the domain $[-L,L]$, does not return to the interval any more. This is in accordance with properties of Lévy flights, which for $\alpha<1$ are transient Blumenthal et al. (1961), i.e. for LF with $\alpha<1$ starting outside a finite interval, there is nonzero probability of not visiting this interval at all. In the limit of $\alpha\to 0$ this probability tends to 1. In the intermediate range of $0<\alpha<1$, returns to the $[-L,L]$ interval are possible, but their probability is smaller than 1. Therefore, MFPT is smaller than (finite) MRT.

The long time asymptotics of the probability $\mathrm{Pr}\left(t,x_{0}\right)$ can be estimated from Eq. (28). In the limit $t\rightarrow\infty$ the function $\exp\left[-\sigma_{0}^{\alpha}k^{\alpha}t\right]$ under the integral becomes very narrow near the point $k=0$ and we can approximately estimate this integral. As a result, we arrive at

[TABLE]

As seen from Eq. (33) the integral (29) diverges for $\alpha\geqslant 1$.

It is interesting to compare the result (33) with a free Brownian diffusion described by the following Smoluchowski-Fokker-Planck equation

[TABLE]

The solution of Eq. (34) with the initial condition (22) is well-known and has the form

[TABLE]

Substituting Eq. (35) in Eq. (21) in the limit $t\rightarrow\infty$ we approximately have

[TABLE]

The asymptotics (33) transforms into asymptotics (36) when $\alpha=2$ because of $\Gamma(1/2)=\sqrt{\pi}$ and $\sigma_{0}=\sqrt{K_{2}}$.

Equations (30) and (32) provide formulas for the mean residence time when $\alpha<1$ and for infinitely long measurement times. Depending on the observation time $t$, the amount of time spent in the $[-L,L]$ interval changes. Therefore, a related question is to estimate how the mean residence (occupation) time grows with the measurement time $t$. The scaling of the average occupation time as a function of the measurement time can be calculated by general properties of Lévy flights. By the definition, the occupation time $T$ is

[TABLE]

where $x(t)$ is a random walker location, $t$ is the measurement time and $\Theta(\cdot)$ is the characteristic function

[TABLE]

The average residence time $\langle T(x_{0})\rangle$ is (for simplicity $x_{0}=0$)

[TABLE]

where $P(x,t|0,0)$ is

[TABLE]

see Eq. (27).

For infinite measurement time ($t\to\infty$) and $0<\alpha<1$ the integral in Eq. (39) is convergent

[TABLE]

and is equal to Eq. (32) representing the fact that for $\alpha<1$ a random walker spends a constant amount of time in the $[-L,L]$ interval, what is the consequence of already discussed transient character of $\alpha$-stable motions with $\alpha<1$, Dybiec et al. (2016).

For finite measurement time $t$ and any $0<\alpha\leqslant 2$

[TABLE]

Eq. (II.3) can be approximated for short and large measurement times. For the short measurement time one has

[TABLE]

In order to find large measurement time asymptotics, the general formula given by Eq. (II.3) needs to be rewritten as

[TABLE]

where $\xi=L/(\sigma_{0}t^{1/\alpha})$ and $\kappa^{\alpha}=k^{\alpha}\sigma_{0}^{\alpha}t$. From Eq. (44) one gets the large $t$ (small $\xi$) asymptotics

[TABLE]

In summary, for $1<\alpha\leqslant 2$ the residence time depends on the measurement time $t$ like

[TABLE]

while for $\alpha<1$

[TABLE]

For $\alpha=1$, after further calculations, one gets

[TABLE]

which for $t\to\infty$ gives logarithmic dependence

[TABLE]

Figure 7 presents the mean residence time $\langle T(x_{0}=0)\rangle$ for $\alpha\in\{0.7,1.5,2.0\}$ (from top to bottom) as a function of the measurement time $t$. Various curves correspond to various interval half-width $L$. The time a particle spends in the interval $[-L,L]$ is always less than the measurement time $t$. Dashed (short time) and solid (long time) lines in Fig. 7 present scaling of the MRT with the measurement time given by Eq. (46) (middle and bottom panels) and Eq. (47) (top panel). In both cases simulation results agree well with theoretical predictions. Altogether, we observe that for Lévy motions with $1<\alpha<2$, the exponent of the asymptotic scaling of MRT is smaller than $1/2$, the value characteristic for a free Brownian diffusion (cf. Barkai (2006) and Eq. (46). This dependence on the stability index $\alpha$ originates from discontinuity of Lévy flight trajectories and asymptotics of the first arrival (hitting) time distribution which is of $t^{1/\alpha-2}$ type, cf. Blumenthal et al. (1961); Chechkin et al. (2003). For $\alpha<1$ the Lévy flight process after leaving the $[-L,L]$ interval does not need to return to its interior, see Blumenthal et al. (1961), but it can be jumping back and forth above the $[-L,L]$ interval. Due to infinite propagation velocity of Lévy flights such jumps above the interval do not contribute to the residence time. Consequently, for $\alpha<1$, the mean residence time saturates at the finite value, as in the $\alpha\to 0$ limit discussed above.

III Lévy Walks

Results presented so far have been constructed for Lévy flights, when after each time-step $\Delta t$ a random walker performs a jump whose length is distributed according to the $\alpha$-stable density. In this type of the random motion model, the integration time step $\Delta t$ scales the jump length distribution in such a way that for sufficiently small $\Delta t$ resulting mean first passage time and stationary states are invariant with respect to the actual value of $\Delta t$. On the other hand, propagation of a random walker performing Lévy flights is characterized by unphysical, infinite velocity. This apparent drawback of the Lévy flights scenario can be resolved by Lévy walk models, which in contrast, assume finite velocity of a jump Shlesinger and Klafter (1986); Zumofen and Klafter (1993).

Here we use a one-dimensional version of the Lévy walk model Zumofen and Klafter (1993); Froemberg et al. (2015); Rebenshtok et al. (2014) for which position is a continuous variable: the jump durations $\mathcal{T}_{i}$ are set to $\mathcal{T}_{i}=|\varsigma_{i}|$ with $\varsigma_{i}$ being i.i.d random variables drawn from the symmetric $\alpha$-stable density, see Eq. (3) and the jump velocity is characterized by a two-state PDF:

[TABLE]

Accordingly, during each jump a particle position’s changes continuously. The jump is finished when the particle travels the total distance $v_{0}\mathcal{T}_{i}$. After completion of a jump, immediately a new jump duration and a jump velocity are generated. In order to calculate the position at time $t$, the whole procedure is repeated $n$ times, where $n$ satisfies $\sum_{i=1}^{n-1}\mathcal{T}_{i}<t<\sum_{i=1}^{n}\mathcal{T}_{i}$. Finally, the position at time $t$ is calculated by adding to the position at time $\sum_{i=1}^{n-1}\mathcal{T}_{i}$, the velocity $v_{0}$ multiplied by the time interval $t-\sum_{i=1}^{n-1}\mathcal{T}_{i}$. Equivalently, the employed Lévy walk model can be defined by drawing the jump lengths $\Delta x_{i}$ from the symmetric $\alpha$-stable density, see Eq. (3). Negative increments $\Delta x_{i}$ correspond to jumps taken to the left, while positive $\Delta x_{i}$ represents jumps to the right. A particle is assumed to move with a constant velocity $v_{0}$ resulting in the jump duration $\mathcal{T}_{i}=|\Delta x_{i}|/v_{0}$ and after finishing the jump, a new jump is immediately generated. Without loss of generality we further set $v_{0}=1$. Note, that in the considered Lévy walk model, the distribution of $\mathcal{T}$ is one-sided, according to the definition $p(\mathcal{T})=p(\varsigma)|d\varsigma/d\mathcal{T}|$ and following Eq. (4) assumes asymptotically the form $p(\mathcal{T})\approx 2\sigma_{0}^{\alpha}\sin(\frac{\pi\alpha}{2})\Gamma(\alpha+1)\pi^{-1}\mathcal{T}^{-(1+\alpha)}$.

Unlike for Lévy flights, the implementation of boundary conditions for Lévy walks is natural thanks to the continuity of their trajectories. Here, every time a random particle crosses the absorbing boundary, it is removed from the system. At the reflecting boundary, every time a particle hits the point, its trajectory is wrapped, i.e. its motion becomes reversed at the boundary and the random walker continues movement in the opposite direction.

III.1 First escape problem

First we study the problem of the escape from the domain restricted by two absorbing boundaries located at $\pm L$. In order to verify how the mean first passage time scales with the system size the MFPT has been numerically estimated for a series of increasing interval half-widths $L$. Fig. 8 shows the MFPT as a function of the interval half width $L$ for $\alpha<1$ (top panel) and $\alpha>1$ (bottom panel). From these plots one can conclude that for large $L$

[TABLE]

Therefore, for Lévy walks with $\alpha>1$, the mean first passage time scales on the interval half-width $L$ in the same manner like for Lévy flights, see Eq. (8). Different scaling is recorded for $\alpha<1$, i.e. in the situation when average jump duration/length is infinite, see Eq. (51) and top panel of Fig. 8.

In order to improve the clarity of the presented figure, in the top panel of Fig. 8, curves corresponding to increasing values of $\alpha$ have been shifted upwards by multiplying MFPTs by a constant factor, otherwise they are superimposed. Furthermore, from obtained $\langle\tau(x_{0})\rangle$ the ratio of mean first passage times $\langle\tau(x_{0})\rangle$s has been calculated

[TABLE]

Fig. 9 presents $R$ for $\alpha<1$ and $R^{-\alpha}$ for $\alpha>1$ (with $x_{0}=0$) which in more details demonstrate how scaling from Eq. (51) is reached. From Eq. (51) it implies that $R\approx 10$ (for $\alpha<1$) and $R^{-\alpha}\approx 10$ (for $\alpha>1$). Indeed such a behavior is visible in Fig. 9 especially for large $L$.

Equation (51) and Fig. 8 suggest that the mean first passage time for Lévy walks with $1<\alpha\leqslant 2$ scales similarly like for Lévy flights, i.e. $\langle\tau(0)\rangle\sim L^{\alpha}$, see Eq. (8). As mentioned above, a main difference between these two categories of free motion is in the finite propagation velocity $v_{0}$ and continuity of trajectories for LW, versus infinite propagation velocity and discontinuity of trajectories for LF. In both cases however, jumps are distributed with an $\alpha$-stable density. Therefore, one can expect that for large interval half-width $L$, Lévy walks can be effectively approximated by Lévy flights with some unknown scale parameter $\sigma_{0}^{\mathrm{LW}}$ or anomalous diffusion constant $K_{\alpha}^{\mathrm{LW}}=\left[\sigma_{0}^{\mathrm{LW}}\right]^{\alpha}$, see Eq. (6). At the same time the finite propagation introduces cut-off for time dependent densities for LW model. Their support is restricted to $[-v_{0}t,v_{0}t]$.

The anomalous diffusion constant $K_{\alpha}^{\mathrm{LW}}$ can be estimated by (Rebenshtok et al., 2014, Eqs. (4) and (42)) and Eqs. (4) and (50) resulting in

[TABLE]

where $\langle\mathcal{T}\rangle$ is the mean waiting time for a next jump (flight duration) and $\sigma_{0}$ is the scale parameter in Eq. (3) and $v_{0}=1$. $\langle\mathcal{T}\rangle$ is finite for $1<\alpha\leqslant 2$. Application of Eq. (53) requires knowledge of $\langle\mathcal{T}\rangle$ which can be estimated numerically or with the help of (Rebenshtok et al., 2014, Eq. (77))

[TABLE]

Formula (54) perfectly approximates $\langle\mathcal{T}\rangle$. For $\alpha\geqslant 1.2$ errors are smaller than $0.3\%$ of the exact value. For $\alpha\to 1$, approximation breaks down reflecting the fact that $\langle\mathcal{T}\rangle$ diverges. Finally one gets the formula for $K_{\alpha}^{\mathrm{LW}}$

[TABLE]

If the Lévy flight approximation to Lévy walks works, the mean first passage time for LW could be estimated from the formula analogous to Eq. (8)

[TABLE]

Estimated values of $K_{\alpha}^{\mathrm{LW}}$, see Eq. (55), are given in Tab. 1. These values have been used to calculate the MFPT according to Eq. (56), see solid lines in the bottom panel of Fig. 8. Therefore, solid lines in the bottom panel of Fig. 8 not only display the $L^{\alpha}$ scaling, but also the MFPT values estimated by the LF approximation to LW, see Eq. (56). The latter works nicely for large $L$ and $\alpha>1$.

In order to measure the quality of LF approximation to LW, the ratio of numerically estimated MFPT for LW, $\langle\tau(0)\rangle$, and value of the MFPT evaluated according to Eq. (56), $\langle\tau(0)\rangle^{\mathrm{LW}}$, has been introduced:

[TABLE]

and analyzed for various ranges of $L$ and $\alpha$, see Fig. 10. With the increasing interval half width $L$ the quality of LF approximation improves resulting in $\tilde{R}\approx 1$. Moreover, approximation of Lévy walks by Lévy flights works better for larger values of the stability index $\alpha$.

Finite propagation velocity results in finite support for LW which is restricted to the interval $[-v_{0}t,v_{0}t]$. At the same time Lévy flights are unconstrained and located at any point on the real line as they propagate with the infinite velocity. Despite this fundamental difference in the propagation velocity Lévy flights seem to well approximate Lévy walks in the central part of the respective propagator (time-dependent PDF). For free Lévy flights the scale parameter, which defines central part of the distribution, grows like $\sigma_{0}^{\mathrm{LF}}t^{1/\alpha}$, where $\sigma_{0}^{\mathrm{LF}}$ is the scale parameter of the jump length density, see Eq. (5). The same dependence of the scale parameter is observed for free Lévy walks, i.e. $\sigma_{0}^{\mathrm{LW}}t^{1/\alpha}$.

Figure 11 presents MFPT for Lévy walks with $L=10$ and $L=100$ as a function of the stability index $\alpha$. Middle and bottom panels present survival probabilities $S(t|x_{0})$ with $x_{0}=0$ for $L=10$ and $L=100$, respectively. From the top panel it is clearly visible that approximation (56) works better for large interval half-width $L$. This is in line with previous considerations. Analogously like for Lévy flights, for a fixed interval half width $L$, the mean first passage time for reflecting-absorbing scheme is larger than for absorbing-absorbing scenario because in the former scenario distance to an absorbing boundary is larger than in the latter scheme. Due to finite particle velocity, for a fixed $L$, the mean first passage time for reflecting-absorbing scenario is not equal to the MFPT for absorbing-absorbing scenario with two times wider interval half-width (results not shown) as it was observed for Lévy flights in the wrapping scenario. Like for Lévy flights, see Fig. 2, survival probabilities have exponential tails. Empty symbols in middle and bottom panels of Fig. 11 correspond to the absorbing-absorbing scenario, while full symbols to the reflecting-absorbing scenario.

At the same time lines present $S(t|x_{0})\approx\exp\left[-t/\langle\tau(x_{0})\rangle\right]$ approximation to the survival probability for absorbing-absorbing (solid lines) and reflecting-absorbing (dotted lines) boundary conditions. Quality of the approximation depends in non trivial way on the stability index $\alpha$ and type of boundary conditions. In particular, such an approximation works better for reflecting-absorbing than absorbing-absorbing setup. The discrepancy between the actual decay and $\exp\left[-t/\langle\tau(x_{0})\rangle\right]$ means that the exponential decay is only asymptotic, since if we had an exponential decay for all times, the law (14) should hold.

Finally, the distribution of number of jumps performed until leaving of the $[-L,L]$ interval has exponential asymptotics (results not shown).

III.2 Stationary states

Figure 12 presents histograms for Lévy walk confined by two reflecting boundaries located at $\pm 10$ (top panel) and corresponding cumulative densities $\mathcal{F}_{x}(x)$ (bottom panel) at $t=10^{3}$. For $\alpha>1$, the system already reached its stationary state which is a uniform probability density $P_{st}(x)=\frac{1}{2L}$, while for $\alpha<1$ a persistent cusp at the $x_{0}=0$ (origin) is still visible. With increasing time $t$ the hight of this peak is decreasing and time dependent density become more and more uniform. In the probability density, the peak visible for $\alpha=0.5$ is responsible for the jump of the cumulative density at $x\approx 0$.

Figure 13 present sample time dependent densities $P(x,t|x_{0},0)$ as heat-maps (top panel), location of maxima of $P(x,t|x_{0},0)$ (middle panel) and maximal values of $P(x,t)$ (bottom panel). Maxima of probability densities originate due to the initial condition, i.e. $x(0)=x_{0}$. Their height is a decaying function of time with the decay rate dependent on the stability index $\alpha$. The slowest decay is observed for $\alpha<1$, when the average jump length is infinite, because the initial peak has smallest number of chances to bifurcate. Putting it differently, height of the maximum of probability density decreases every time a jump direction and jump length are generated. Therefore, the number of jumps (determined by $\alpha$) defines the decay rate of the initial condition. Fig. 13 presents results for $\alpha=0.5$ with $x_{0}=0$ and $x_{0}=6$ because such a choice of parameters nicely shows general properties of time dependent densities $P(x,t|x_{0},0)$. Analogously, positions of maxima are not fixed but they constantly, ballistically move and bifurcate due to finite propagation velocity $v_{0}=1$ and equal probability of jumps to the left and right. The initial cusp splits into two parts moving to the right and left. For $\alpha<1$ average jump duration is infinite, therefore these parts return to the initial position after $\Delta t=2L/v_{0}$ (for $x_{0}=0$) or $\Delta t=4L/v_{0}$ (for $x_{0}\neq 0$). This effect is nicely visible in the middle panel of Fig. 13. Multiple points for the same $t$ in the bottom panel correspond to bimodal time dependent densities. The non-monotonous, point-like amplification of peaks takes place when two peaks moving in opposite directions meet in one point. At this point height of the peak is twice the background. In the limit of large times $t$ time dependent densities converge to the uniform stationary density, i.e. $P_{st}(x)=1/2L$.

IV Summary and conclusions

Both Lévy walks and Lévy flights assume that a random walker performs long jumps distributed according to a heavy-tailed, power-law density. The main difference between both scenarios is in continuous trajectories and finite propagation velocity of Lévy walks versus discontinuous trajectories and infinite propagation velocity of Lévy flights. Nevertheless, due to the same type of the jump length density both scenarios are deeply related and need to be compared in great detail. Such a comparison between behavior of LF and LW, with respect to a class of important observables like mean first passage time, survival probabilities and stationary states provided the main motive of current research. We have verified when both models are similar and when they differ. In addition to comparison of two classical random walk schemes we studied the problem of posing boundary conditions for Lévy flights and Lévy walks.

The mean first passage time for Lévy flights scales asymptotically as $L^{\alpha}$ with the interval-half width, what is especially well visible for $x_{0}=0$. The same scaling is observed for Lévy walks with $\alpha>1$ and large interval half-width $L$. For $\alpha<1$ the mean first passage time is proportional to $L$. Eq. (51), which is the main outcome of this part of the manuscript, indicates that for $\alpha>1$ Lévy flights with properly adjusted scale parameter can approximate Lévy walks with respect to the analysis of mean first passage time.

In addition to the escape of Lévy flights from finite intervals we have analyzed the problem of the residence time, i.e. the fraction of time which unbounded process spends in the prescribed part of the line, e.g. a finite interval. The mean residence time of LF as a function of the measurement time $t$ has a universal short time scaling $\langle T(0)\rangle\sim t$ and a non-universal long time scaling $t^{1-1/\alpha}$ (for $1<\alpha\leqslant 2$). For $\alpha<1$, due to discontinuity of trajectories of Lévy flights and their transient character, the mean residence time saturates at a constant value as measurement time goes to infinity. Results regarding the mean residence time, Eq. (32), and dependence of the mean residence time on the measurement time $t$, Eqs. (46)–(48), constitute key results of this part of the manuscript.

Problem of formulation of boundary conditions on a single trajectory level for Lévy flights is not fully resolved. While it is known how to treat the absorbing boundaries, it is not uniquely defined how to implement reflecting boundary conditions. We have compared two scenarios: wrapping (reversing) and stopping (pausing) showing that the former results in uniform stationary states for a motion restricted by two reflecting boundaries, while the latter scenario gives the same stationary states as impenetrable boundary conditions Denisov et al. (2008) and steep potential wells Dubkov and Spagnolo (2007). In the case of Lévy walks, the problem of posing reflecting boundary conditions is more apparent. Finite propagation velocity and continuous paths suggest that trajectory should be wrapped along reflecting boundaries. As a consequence stationary states for a motion restricted by two reflecting boundaries are uniform. Nevertheless, Lévy walks display slowly decaying memory about initial conditions especially in situations when mean jump length (flight time) is divergent.

Our work clarifies the cases where Lévy flights can be used as an approximation to Lévy walks. Finite propagation velocity of LW makes them inertial, i.e. a test particle moves with finite kinetic energy, and collisions leading to instantaneous alternations of velocities are fully elastic. In consequence LW trivially fulfills Newton equations. This is not however the case of LFs: Inertial Lévy flights Sokolov et al. (2011); Lü and Bao (2011) break the equipartition theorem and lead to inter-dependence between the position and velocity, which is manifested by the nontrivial joint distribution of both variables. Moreover, in view of thermodynamic interpretation, the Langevin equation with non-Gaussian Lévy noise, Eq. (19) describes change of the position of a particle under the action of nonequilibrated external forcing which does not fulfill the standard fluctuation-dissipation theorem Dybiec et al. (2012); Kuśmierz et al. (2014).

Within the effort of explaining real life phenomena exhibiting anomalous diffusion, Langevin or generalized Langevin description is frequently a natural choice Zaburdaev et al. (2015); Lubashevsky et al. (2009); Ferreira et al. (2012); Kneller (2011); Bao et al. (2006). However, a full correspondence of that picture with Lévy walks and Lévy flights dynamics still calls for further systematic studies.

V Acknowledgments

This project has been supported in part (BD and EGN) by the grant from National Science Center (2014/13/B/ST2/02014). EB thanks the Israel Science Foundation for support. Computer simulations have been performed at the Academic Computer Center Cyfronet, Akademia Górniczo-Hutnicza (Kraków, Poland).