Phase Diagram in Stored-Energy-Driven Lévy Flight

Phase diagram based on the mean square displacement (MSD) and the distribution of diffusion coefficients of the time-averaged MSD for the stored-energy-driven Lévy flight (SEDLF) is presented. In the SEDLF, a random walker cannot move while storing energy, and it jumps by the stored energy. The SEDLF shows a whole spectrum of anomalous diffusions including subdiffusion and superdiffusion, depending on the coupling parameter between storing time (trapping time) and stored energy. This stochastic process can be investigated analytically with the aid of renewal theory. Here, we consider two different renewal processes, i.e., ordinary renewal process and equilibrium renewal process, when the mean trapping time does not diverge. We analytically show the phase diagram according to the coupling parameter and the power exponent in the trapping-time distribution. In particular, we find that distributional behavior of time-averaged MSD intrinsically appears in superdiffusive as well as normal diffusive regime even when the mean trapping time does not diverge.


Introduction
In normal diffusion processes, the diffusivity can be characterized by the diffusion coefficient in the mean square displacement (MSD). However, many diffusion processes in nature show anomalous diffusion; that is, the MSD does not grow linearly with time but follows a sublinear or superlinear growth with time, where x t is a position in one-dimensional coordinate, t is time, and · means an ensemble average. In particular, anomalous diffusion in biological systems has been found by singleparticle tracking experiments [13,19,20,22,40,43,44]. Thus, the power-law exponent β is one of the most important quantities characterizingthe underlying diffusion process. Especially, anomalous diffusion with β < 1 is called subdiffusion and that with β > 1 is called superdiffusion.
Although anomalous diffusion can be characterized by the power-law exponent β in the MSD, the exponent cannot reveal the underlying physical nature in itself. This is because the same power-law exponent in the MSD does not imply that the physical mechanism in the anomalous diffusion is also the same. Therefore, clarifying the origin of anomalous diffusion is an important subject, and many researches on this issue have been conducted extensively [23,29,30,32,41]. One of the key properties characterizing anomalous diffusion is ergodicity, i.e., time-averaged observables being equal to a constant (the ensemble average). In some experiments [19,20,22,40,44], (generalized) diffusion coefficients for time-averaged MSDs show large fluctuations, suggesting that ergodicity breaks.
In stochastic models of anomalous diffusion, continuous-time random walk (CTRW) shows a prominent feature called distributional ergodicity [21,27,33,34]; that is, the timeaveraged observables obtained from single trajectories do not converge to a constant but the distribution of such time-averaged observables converges to a universal distribution (convergence in distribution). More precisely, the distribution of the time-averaged MSD (TAMSD), which is defined by converges to the Mittag-Leffler distribution of order α [5,21]. This statement can be represented by δ 2 ( ; t) for a fixed ( t), where M α is a random variable with the Mittag-Leffler distribution of order α. We note that the diffusion coefficients in the TAMSDs are also distributed according to the Mittag-Leffler distribution because the TAMSD shows normal diffusion, i.e, δ 2 ( ; t) D t [33,34], where we refer to D t as the diffusion coefficient. In stochastic models, such distributional behavior originates from the divergent mean trapping time. In diffusion in a random energy landscape such as a trap model [11], the trapping-time distribution follows a power law, w(τ ) ∝ τ −1−α , if heights of the energy barrier are distributed according to the exponential distribution. The exponent α smaller than 1 implies a divergence of the mean. In dynamical systems, this divergent mean brings an infinite measure [4]. Thus, such distributional behavior of time-averaged observables is also shown in infinite ergodic theory [1,5]. Moreover, large fluctuations of time-averaged observables, related to distributional ergodicity, have been observed in biological experiments [19,20,22,40,44] as well as quantum dot experiments [12,39].
Recently, we have shown that the distribution of the diffusion coefficients of the TAMSDs in the stored-energy-driven Lévy flight (SEDLF) is different from that in CTRW [6]. The SEDLF is a CTRW with jump lengths correlated with trapping times [24,26,28]. One of the most typical examples for such a correlated motion can be observed in Lévy walk [38]. However, Lévy walk and SEDLF are completely different stochastic processes in that a random walker cannot move while it is trapped in SEDLF, whereas it can move with constant velocity in Lévy walk. In other words, in SEDLF, a random walker does not move while storing a sort of energy, and it jumps using the stored energy (see Fig. 1). When the trappingtime distribution follows a power-law, the jump length distribution also follows a power-law, the same as in Lévy flight. Since we consider a power-law trapping-time distribution, we refer to our model as the stored-energy-driven Lévy flight. We note that the MSD does not diverge in SEDLF, whereas it always diverges in Lévy flight. Although the ensemble-averaged MSDs show subdiffusion as well as superdiffusion, the TAMSDs always increase linearly with time in SEDLF [6]. This behavior is completely different from that in Lévy walk [2,17]. Moreover, the distribution of the TAMSD with a fixed converge to a time-independent distribution which is not the same as a universal distribution in CTRW (Mittag-Leffler distribution): where Y α,γ is a random variable and γ is the coupling parameter between trapping time and jump length. We note that such coupling effects become physically important in turbulent diffusion [38], diffusion of cold atoms [10], and nonthermal systems such as cells [13,43]. In such enhanced diffusions, it has been known that the coupling between jump lengths and waiting times follows a power-law fashion like SEDLF [10,38], although a particle is always moving, which is different from SEDLF. Furthermore, it will also be important in complex systems such as finance [31] and earthquakes [14,25] because jump lengths are correlated with the waiting times in such systems.
In terms of an ensemble average, SEDLF exhibits a whole spectrum of diffusion: sub-, normal-, and super-diffusion, depending on the coupling parameter [6,24,26,28].Because distributional behavior of the time-averaged observables such as the diffusion coefficients in SEDLF is different from that in CTRW, it is important to construct a phase diagram in terms of the power-law exponent of the MSD as well as the form of the distribution function of the TAMSD. Here, we provide the phase diagram for the whole parameters range in SEDLF.

Model
SEDLF is a cumulative process, which is a generalization of a renewal process [15]. Equivalently, SEDLF is a CTRW with a non-separable joint probability of trapping time and jump length. Therefore, the SEDLF can be defined through the joint probability density function (PDF) ψ(x, t), where ψ(x, t)dxdt is the probability that a random walker jumps with length [x, x + dx) just after it is trapped for period [t, t + dt) since its previous jump [11,37]. Here, we consider the following joint PDF where w(t) is the PDF of trapping times and 0 ≤ γ ≤ 1 is a coupling strength. This kind of coupling has been introduced in [24,37]. The SEDLF with γ = 0 is just a separable CTRW.
In addition, we consider that the PDF of trapping times follows a power law: as t → ∞. Here, α ∈ (0, 2) is the stable index, a constant c 0 is defined by c 0 = c/| (−α)| with a scale factor c. We note that the mean trapping time diverges for α ≤ 1. For γ > 0, the PDF of jump length also follows a power law: Thus, the second moment of the jump length diverges for 2γ ≥ α. Because Lévy flight also has a power law distribution of jump length, we call this random walk the stored-energy-driven Lévy flight. In numerical simulations, we set the PDF of the trapping time as w(t) = αt −1−α (t ≥ 1). Thus, the jump length PDF is given by l(x) = α/(2γ |x| 1+α/γ ) from Eq. (7), and Because the mean trapping time is finite for α > 1, we consider two typical renewal processes; ordinary renewal process and equilibrium renewal process [15]. Equilibrium renewal process is assumed to start −t a = −∞ (see Fig. 1). The PDF of the first jump length x and the first apparent trapping time (the forward recurrence time) t, ψ 0 (x, t), is given by where μ is the mean trapping time. See Appendix for the derivation. Generally, the first (true) trapping time is longer than the first apparent trapping time. Integrating Eq. (8) in terms of x, we have the PDF of the first apparent trapping time: Note that the Eq. (9) is consistent with the result obtained in renewal theory, i.e., the PDF of the forward recurrence time [15]. Thus, the joint PDF of the first jump length and the first apparent trapping time, ψ 0 (x, t), is not given by the form in Eq. (5). This is because the first jump length x is determined by the time elapsed since a random walker's last jump (the first true trapping time) and thus it is not directly related to the first apparent trapping time, i.e., the time elapsed since the beginning of the measurement at t = 0. For ordinary renewal process, we just set w 0 (t) = w(t) and ψ 0 (x, t) = ψ(x, t). For α ≤ 1, we only consider an ordinary renewal process because there is no equilibrium ensemble due to divergent mean trapping time which causes aging [7,9,36].

Generalized Renewal Equation
The spacial distribution P(x, t) of SEDLFs with starting from the origin satisfies the generized renewal equations: where Q(x, t)dtdx is the probability of a random walker reaching an interval is the probability of being trapped for longer than time t just after a renewal (after the measurement starts at t = 0). (t) and 0 (t) are defined as follows: Then, the Laplace transforms of these functions are given bŷ In an ordinary renewal process, ψ 0 (x, t) is the same as ψ(x, t). On the other hand, the first jump length x is not determined by the first apparent trapping time t in an equilibrium renewal process, while the first jump length is not independent of the first apparent trapping time.
Fourier-Laplace transform with respect to space and time (x → k and t → s, respectively), defined byP givesP In what follows, we use the notations P o (x, t) andP o (k, s) for the ordinary renewal process, and P eq (x, t) andP eq (k, s) for the equilibrium renewal process. For ordinary renewal process, i.e., ψ 0 (t) = ψ(t) and 0 (t) = (t), we have the following generalized renewal equation in the Fourier and Laplace space: where we used Eq. (14). For equilibrium renewal process (α > 1), we havê where we used Eq. (14) andψ The derivation of Eq. (19) is shown in Appendix. (17) and (18)]. Now, we derive the explicit forms of these functions. From Eq. (5),ψ(k, s) is given bŷ Note thatψ(0, s) =ŵ(s). In addition, from Eq. (6), the asymptotic behavior of the Laplace transformŵ(s) for s → 0 is given bŷ

Mean Square Displacement
The asymptotic behavior of the moments of position x t for t → ∞ can be obtained using = 0, x t = 0 for both renewal processes.
In ordinary renewal process, the Laplace transform of the second moment, i.e., the ensemble-averaged MSD (EAMSD), is given by where the ensemble average . . . o is taken with respect to an ordinary renewal process. For α ∈ (0, 1), we obtain the EAMSD [6]: where we used t 2γ = l 2 = ∞ −∞ x 2 l(x)dx when l 2 < ∞. For α = 1, the EAMSD is given by Finally, for α ∈ (1, 2), the EAMSD is given by These results are consistent with a previous study [24]. We note that the EAMSD for γ = 1 is smaller than that in Lévy walk, whereas the scaling exponent 3 − α is the same as that in Lévy walk [45]. This is because SEDLF is a wait and jump model, while Lévy walk is a moving model. In equilibrium renewal process (α > 1), Eq. (25) is valid only for 0 < 2γ < α − 1, otherwise the second moment of the first jump length diverges, i.e., the EAMSD diverges. This is very different from Lévy walk process because there exists an equilibrium renewal process in Lévy walk with α > 1. Becausê ψ (0, 0) = l 2 for 0 < 2γ < α − 1, the EAMSD is given by In SEDLF, the leading order of the EAMSD in an ordinary renewal process is the same as that in an equilibrium renewal process (α > 2γ + 1). On the other hand, in Lévy walk, the proportional constant of the EAMSD in a non-equilibrium ensemble such as an ordinary renewal process differs from that in an equilibrium one [16][17][18]45]. We note that the TAMSD coincides with the EAMSD in an equilibrium ensemble as the measurement time goes to infinity. Significant initial ensemble dependence of statistical quantity has been also observed in non-hyperbolic dynamical systems [3].

Time-Averaged Mean Square Displacement
In normal Brownian motion, the TAMSD defined by Eq. (2) converges to the MSD with an equilibrium ensemble: Such ergodic property does not hold in various stochastic models of anomalous diffusion such as CTRW and Lévy walk [17,21,27]. Here, we derive the TAMSD in the SEDLF. In wait and jump random walks with random waiting times such as CTRW and SEDLF, the TAMSD can be represented using the total number of jumps, denoted by N t , and h k = l 2 k + 2 k−1 m=1 l k l m θ( − t k + t m ) [6,34,35]: where l k is the k-th jump length, t k is the time when the k-th jump occurs, and θ(t) is a step function, defined by θ(t) = 0 for t < 0 and t otherwise. For γ > 0, one can show that . Therefore, the TAMSD can be written as where D t = N t k=0 l 2 k /t. We note that the relation (28) does not hold if the random walker moves with constant speed as in Lévy walk because (x t+ − x t ) 2 is simply zero in SEDLF but not in Lévy walk. In fact, the TAMSD does not increase linearly with time in Lévy walk [16][17][18].
To investigate an ergodic property of the time-averaged diffusion coefficient D t , we derive the PDF P 2 (z, t) of Z t ≡ N t k=0 l 2 k . Because l 2 k and N t are mutually correlated, we use the generalized renewal equation for Z t : where the joint PDF φ(z, t) is given by The joint PDF of the first squared jump length z = l 2 0 and apparent trapping time t, φ 0 (z, t), is given by φ 0 (z, t) = φ(z, t) for the ordinary ensemble, whereas for equilibrium ensemble, which can be derived in the same way as the derivation of ψ 0 (x, t) given in Appendix. The double Laplace transform with respect to z and time t is defined bŷ From the generalized renewal equations (30) and (31), we obtain where the double Laplace transform of φ(z, t) is given bŷ and that of φ 0 (z, t) is given by φ 0 (z, t) = φ(z, t) for the ordinary process, and bŷ for the equilibrium process. Note thatφ(0, s) =ŵ(s).

Ordinary Renewal Process
In ordinary renewal process, the Laplace transformP o 2 (k, s) is given bŷ Thus, we have the Laplace transform of Z t o as follows: where we used −φ (0, s) = −ψ (0, s) and Eq. (21). Then, averaging Eq. (29) over an ordinary ensemble, we have Therefore, using Eqs. (22)- (24), we have the leading terms of the mean diffusion coefficient, for 0 < α < 1, for α = 1, and for 1 < α < 2. It follows that the mean diffusion coefficient diverges as t goes to infinity for 1 < α < 2 and α < 2γ < 2. Similarly, the Laplace transform of Z 2 t is given by It follows that the leading order of the second moment of D t is given by for 0 < α < 1, for α = 1, and μ . Now, we study the relative standard deviation (RSD) of D t , σ EB (t) ≡ D 2 t − D t 2 / D t [21,33,34], to measure the ergodicity breaking. First, for 0 < α < 1, RSD σ EB (t) does not converge to zero but to a finite value as t → ∞. Therefore, the diffusion coefficients remain random even when the measurement time goes to infinity [6]. Second, for α = 1, we expect usual ergodic behavior for 0 < 2γ ≤ 1, because the RSD goes to zero as t → ∞. However, for 1 < 2γ < 2, the RSD diverges as t → ∞: Finally, for 1 < α < 2, we have Thus, TAMSDs show ergodic behavior when the parameters satisfy 0 < 2γ < α+1 2 because the RSD goes to zero as t → ∞. On the other hand, the RSD converges to a finite value for γ = α+1 4 , and diverges for α+1 2 < 2γ ≤ 2. Numerical simulations suggest that this divergence of the RSD for the case of α+1 2 < 2γ ≤ 2 will be attributed to a power-law with divergent second moment in the PDF of D t / D t o (see Fig. 2). Because the RSD is defined using the second moment of D t , it diverges, whereas the PDF converges to a power-law distribution. Thus, the RSD,σ EB (t), is not helpful to characterize the ergodicity breaking in this case. Using the relative fluctuation defined by R(t) ≡ |D t − D t | / D t [8,42] instead of the RSD, we can clearly see the ergodicity breaking in the case of α+1 2 < 2γ ≤ 2: we numerically found that the relative fluctuation, R(t) = |D t / D t o − 1| o , converges to a constant as t → ∞ because D t / D t o does not converge to one but converges in distribution. Thus, the ergodicity in TAMSD breaks down for α+1 2 ≤ 2γ ≤ 2. For α > 1, the asymptotic behavior of the Laplace transform of Z n t o at s → 0 (n > 1) is given by It follows that the leading order of the nth moment of D t is given by By numerical simulations, we confirm that the scaled diffusion coefficient D t / D t o converges in distribution to a random variable S α,γ : for α+1 2 < 2γ ≤ 2. The distribution depends on γ and α (see Fig. 2). We note that the distribution obeys a power-law with divergent second moment because all the n-th moments (n > 1) of D t / D t o diverge as t goes to infinity. In fact, as shown in Fig. 2, the power-law exponents are smaller than 3.
For α ≤ 1, we obtained all the higher moments of D t [6]. In particular, for 2γ < α, all the moments are given by Therefore, the distribution of the scaled diffusion coefficient converges to the Mittag-Leffler distribution: where Moreover, for 2γ > α, the distribution of D t D t o also converges to a time-independent nontrivial distribution as t → ∞ [6]:

Equilibrium Renewal Process
In an equilibrium renewal process, i.e., 1 < α < 2 and 0 < 2γ < α − 1, the Laplace transformP eq 2 (k, s) is given bŷ Thus, we have the Laplace transform of Z t eq as follows: Averaging Eq. (29) over equilibrium ensemble, we have where the second moment of the first jump length is finite for 0 < γ < α − 1, otherwise the ensemble average of the TAMSD diverges. The second moment of the diffusion constant is also derived in the similar way: and thus we have for 0 < γ < α − 1. From these results, RSD is given by

Discussion
We have shown the phase diagram based on the power-law exponent of anomalous diffusion and the distribution of TAMSDs in SEDLF. The results are summarized in Fig. 3. Although SEDLF is closely related to CTRW, Lévy walk, and Lévy flight, its statistical properties on anomalous diffusion are different form them. In particular, while the visit points in SEDLF are the same as the turning points of a random walker in Lévy walk [38], a random walker cannot move while it is trapped in SEDLF, which is completely different from Lévy walk. This discrepancy makes the scaling of TAMSD different. In fact, TAMSDs in SEDLF increase linearly with time even when the EAMSD shows superdiffusion. On the other hand, TAMSDs show superlinear scaling (superdiffusion) in Lévy walk. Because a particle is always moving and there is a power-law coupling between waiting times and moving length in turbulent diffusion and diffusion of cold atoms [10,38], a moving model of SEDLF can be applied to them. On the other hand, a wait and jump model (SEDLF) will be important in finance and earthquakes. In particular, SEDLF will be applied to describe dynamics of energy released in earthquake because energy is gradually accumulated and released in earthquake. Since TAMSDs can not be represented by Eq. (28) in a moving model such as Lévy walk and a moving model of SEDLF, to investigate ergodic properties of TAMSDs in a moving model of SEDLF is left for future work.

Conclusion
In conclusion, we have shown the phase diagram in SEDLF for a wide range of parameters, where the EAMSD shows normal diffusion, subdiffusion and superdiffusion, and the distribution of the TAMSD depends on the power-law exponent α of the trapping-time distribution as well as the coupling parameter γ . We consider two typical processes: ordinary renewal process and equilibrium renewal process. An equilibrium distribution for the first renewal time (the forward recurrence time) exists in renewal processes when the mean of interoccurrence time between successive renewals does not diverge. However, even when the mean does not diverge, an equilibrium distribution does not exist in the SEDLF because of divergence of the second moment of the first jump length. Therefore, we have found that the TAMSDs remain random in some parameter region even when the mean trapping time does not diverge. In particular, it is interesting to note that this distributional ergodicity is observed even when the EAMSD shows a normal diffusion, i.e., α+1 2 < 2γ < α and 1 < α < 2. In this regime, both the mean trapping time and the second moment of jump length are finite. Therefore, this result provides a novel route to the distributional ergodicity, because so far the distributional ergodicity has been found only in systems with the divergent mean trapping time or the divergent second moment of jump length, which break down the law of large numbers.
Acknowledgments This work was partially supported by Grant-in-Aid for Young Scientists (B) (Grant no. 26800204).
Open Access This article is distributed under the terms of the Creative Commons Attribution License which permits any use, distribution, and reproduction in any medium, provided the original author(s) and the source are credited.

Appendix: Derivation of ψ 0 (x, t)
Here, we derive the joint PDF ψ 0 (x, t) of the first jump length x and the first apparent trapping time t for an equilibrium renewal process (α > 1). Let ψ 0,n (x, t; t a ) be the joint PDF that the first jump and the first apparent trapping time after time t a given that the number of jumps in [0, t a ] is n. The joint PDF can be represented by ψ 0,n (x, t; t a ) = 1 2 δ(x − (t n+1 − t n ) γ )δ(t − (t n+1 − t a ))I (t n < t a < t n+1 ) The Fourier and double Laplace transform with respect to x, t and t a is given bŷ Through the inverse Fourier-Laplace transforms ofψ 0 (k, s), we obtain Eq. (8).