Langevin picture of Lévy walk in a constant force field
Yao Chen Xudong Wang Weihua Deng












School of Mathematics and Statistics, Gansu Key Laboratory
of Applied Mathematics and Complex Systems, Lanzhou University, Lanzhou 730000,
P.R. China





Lévy walk is a practical model and has wide applications in various fields. Here we focus on the effect of an external constant force on the Lévy walk with the exponent of the power-law distributed flight time α∈(0,2). We add the term Fη(s) (η(s) is the Lévy noise) on a subordinated Langevin system to characterize such a constant force, being effective on the velocity process for all physical time after the subordination.
We clearly show the effect of the constant force F on this Langevin system and find this system is like the continuous limit of the collision model. The first moments of velocity processes for these two models are consistent. In particular, based on the velocity correlation function derived from our subordinated Langevin equation, we investigate more interesting statistical quantities, such as the ensemble- and time-averaged mean squared displacements.
Under the influence of constant force, the diffusion of particles becomes faster. Finally, the super-ballistic diffusion and the non-ergodic behavior are verified by the simulations with different α.












Langevin picture of Lévy walk in a constant force field
    Weihua Deng
    December 30, 2023
=======================================================





§ INTRODUCTION


The movements of particles are often in an external potential.
The particles in complex disordered systems influenced by the external force generally  exhibits anomalous diffusion behavior <cit.>, which is characterized by the nonlinear evolution in time of the ensemble-averaged mean squared displacement (EAMSD), i.e.,

    ⟨(Δ x(t))^2⟩=⟨[x(t)-⟨ x(t)⟩]^2⟩ ∝ t^β    (β≠1),

which represents subdiffusion for 0<β<1 and superdiffusion for β>1; for the case β=2, it is called ballistic diffusion.

Just as the brackets ⟨⋯⟩ in Eq. (<ref>) shows, the EAMSD is the statistical average over a large amount of stochastic realizations and thus it is not easy to be measured in experiments. In 1914, Nordlund determined the diffusion coefficients of the traced droplets from separate analysis of single trajectory <cit.>. After that, people started to pay attention to evaluate the recorded time series in terms of the time-averaged mean squared displacement (TAMSD) <cit.>. Especially, for the case that an external force makes a nonzero mean value of the displacement, the TAMSD is defined as
<cit.>:

    δ^2(Δ)   =1/T-Δ∫_0^T-Δ  [x(t+Δ)-x(t)  
              -⟨ x(t+Δ)-x(t)⟩]^2dt,

where Δ is the lag time; it separates the displacement between trajectory points and is much shorter than measurement time T to obtain good statistical properties. The single particle tracking techniques have been widely employed to study diffusion of particles in living cell <cit.>.

One representative model to describe diffusion phenomenon is Lévy walk <cit.>. This model characterizes particle motion with finite velocity, and seems to be more reasonable and originally characterized by coupled continuous time random walk (CTRW) <cit.>, in which the probability density functions (PDFs) of jump length and flight time are coupled through a constant velocity. Depending on the power exponent α of PDF of flight time, Lévy walk could describe ballistic diffusion (0<α<1), sub-ballistic superdiffusion (1<α<2), and normal diffusion (α>2). Later, the equivalent Langevin picture of Lévy walk is presented in Refs. <cit.>. Since the Lévy walk model is especially practical and can depict multiple types of diffusion behaviors, it has wide applications in various fields, not only in the tracking studies of animals or humans <cit.>, but also the anomalous superdiffusion of cold atoms in optical lattices <cit.>, endosomal active transport within living cells <cit.>, etc.

Based on the statistical quantities—EAMSD and TAMSD, the ergodic behavior, an important property of stochastic process, can be investigated. The free Lévy walk has shown a very special phenomenon about the ergodic property. It says that the free Lévy walk presents “ultraweak” non-ergodic behavior for the power exponent α∈(1,2) <cit.>. The word “ultraweak” means that the TAMSD and EAMSD only differ by a constant factor independent of the lag time Δ for a long measurement time. While for the case α∈(0,1) <cit.> with divergent first moment of flight times, the TAMSD is not self-averaged when the measurement time tends to infinity, although the ensemble-averaged TAMSD only differs with EAMSD by a constant factor.

In this paper, we pay attention to how the Lévy walk model reacts to an external constant force in physical time. There has been a collision model proposed in Ref. <cit.> describing such a motion. It assumes that the velocity of a free test particle is changed after the collision with the surrounding bath particles.
If the model is under the influence of a constant force F, it takes effect between two successive collisions; the first moments of velocity and displacement of the model are derived in Ref. <cit.>.

Here, we deal with this problem by building its Langevin picture, and detailedly discuss the way of adding the external force which can make an effect on the process for the whole physical time. We establish a Langevin system coupled with a subordinator to characterize this process and then find it looks like a continuous limit of the collision model.
The first moments of velocity and displacement are verified to be consistent in these two models. Besides, we can make use of the advantage of Langevin equation to calculate the velocity correlation function of the concerned process, which has not been derived before.

The method to obtain the velocity correlation function is technical, where the four-point PDF of the inverse subordinator has to be used.
We establish the relationship between the subordination <cit.> and the renewal theory <cit.>, which enhances the understanding of these theories and greatly simplifies the calculations.
Finally, the velocity correlation function is found to be a sum of two independent terms. Then the EAMSD and TAMSD can be directly obtained by applying the generalized Green-Kube formula <cit.>.





The structure of this paper is as follows. In Sec. <ref>, we review the subordinator and the inverse subordinator, and then build their relationship with the renewal theory. In Sec. <ref>, we establish the Langevin picture of Lévy walk in the presence of a constant force in physical times. Then we calculate the first moment, the velocity correlation function, EAMSD and TAMSD in Sec. <ref>–<ref>, respectively.
Finally, we make the summaries in Sec. <ref>. The detailed derivations of some results are presented in Appendix.






§ SUBORDINATOR AND INVERSE SUBORDINATOR

Subordinator is a non-decreasing Lévy process <cit.> and can be regarded as a stochastic model of time evolution. Especially, it helps to describe different kinds of subdiffusion <cit.> or superdiffusion <cit.> processes when coupled with an overdamped or underdamped Langevin equation. In order to characterize the power-law distributed flight time of Lévy walk, the subordinator t(s) in this paper is taken to be α-dependent (0<α<2) one with the characteristic function g(λ,s):=⟨e^-λ t(s)⟩=e^-sΦ(λ), where Φ(λ)= λ^α for 0<α<1 <cit.> and Φ(λ)= τ_0/(α-1)λ-τ_0|Γ(1-α)| λ^α for 1<α<2 <cit.>.  The two-point PDF of the subordinator t(s) can be expressed as

    g(t_1,s_1;t_2,s_2)= ⟨δ(t_1-t(s_1)) δ(t_2-t(s_2)) ⟩.

By virtue of the stationary and independent increments of subordinator t(s), this two-point PDF in Laplace space (t_1 →λ_1, t_2→λ_2) is given by <cit.>

    g(λ_1,s_1;λ_2,s_2)
           =Θ(s_2-s_1)  e^-(s_2-s_1)Φ(λ_2)  e^-s_1Φ(λ_1+λ_2)
        +Θ(s_1-s_2)  e^-(s_1-s_2)Φ(λ_1)  e^-s_2Φ(λ_1+λ_2),

where Θ(x) is the Heaviside step function: Θ(x)=1 for x>0, Θ(x)=0 for x<0, and Θ(x=0)=1/2.

The inverse α-dependent subordinator s(t):=inf_s>0{s:t(s)>t} is defined as the first-passage time of the subordinator t(s) with the two-point PDF

    h(s_1,t_1;s_2,t_2)= ⟨δ(s_1-s(t_1)) δ(s_2-s(t_2)) ⟩.

The specific expression of h(·) can be obtained by taking the partial derivatives toward s_1,s_2 in the following equation <cit.>:

    ⟨Θ(s_2-s(t_2))Θ(s_1-s(t_1))⟩
        = 1-⟨Θ(t_2-t(s_2))⟩
           -⟨Θ(t_1-t(s_1))⟩
     +⟨Θ(t_2-t(s_2))Θ(t_1-t(s_1))⟩,

which links the two-point PDF of inverse subordinator to the corresponding subordinator. Then performing the Laplace transform (t_1→λ_1,t_2→λ_2), we obtain the two-point PDF h(·) of the inverse subordinator s(t) in Laplace space as <cit.>

    h(s_1,λ_1;s_2,λ_2)  
       =∂/∂ s_1∂/∂ s_21/λ_1λ_2 g(λ_1,s_1;λ_2,s_2) 
       =δ(s_2-s_1)Φ(λ_1)+Φ(λ_2)-Φ(λ_1+λ_2)/λ_1λ_2  e^-s_1Φ(λ_1+λ_2)
          +Θ(s_2-s_1)Φ(λ_2)(Φ(λ_1+λ_2)-Φ(λ_2))/λ_1λ_2
          × e^-s_1Φ(λ_1+λ_2) e^-(s_2-s_1)Φ(λ_2)
          +Θ(s_1-s_2)Φ(λ_1)(Φ(λ_1+λ_2)-Φ(λ_1))/λ_1λ_2
          × e^-s_2Φ(λ_1+λ_2) e^-(s_1-s_2)Φ(λ_1).

The normalization of h(s_2,λ_2;s_1,λ_1) can be verified through the equality ∫_0^∞∫_0^∞  h(s_2,λ_2;s_1,λ_1)ds_1 ds_2=(λ_1λ_2)^-1.
Especially, the first term containing δ(s_2-s_1) in Eq. (<ref>) contributes to the part that no renewal happens between time t_1 and t_2. More precisely, let p_0(t_1,t_2) be the probability of no renewal happens between time t_1 and t_2, and ψ(λ) the Laplace transform of the waiting time distribution between two consecutive renewals. Then the double Laplace transform of p_0(t_1,t_2) is <cit.>

    p_0(λ_1,λ_2)
         ≃1+ψ(λ_1+λ_2)-ψ(λ_1)-ψ(λ_2)/λ_1λ_2(1-ψ(λ_1+λ_2))
       ≃Φ(λ_1)+Φ(λ_2)-Φ(λ_1+λ_2)/λ_1λ_2Φ(λ_1+λ_2),

where we have used the relation ψ(λ)=g(λ,s)|_s=1≃1-Φ(λ) for small λ in the second line.
On the other hand, the result in Eq. (<ref>) can be directly obtained by performing the double integrals with respect to s_1 and s_2 on the first term of the right-hand side of Eq. (<ref>).
Correspondingly, another two terms in Eq. (<ref>) contribute to the probability that renewal happens between time t_1 and t_2 for t_1<t_2 and t_1>t_2, respectively. Explicitly presenting the relationship between the arguments of subordinator and renewal theory helps us better understand these two kinds of theories. Besides, it will be the method for us to simplify some derivations in the following sections.

Performing inverse Laplace transform of Eq. (<ref>) yields <cit.>

    p_0(t_1,t_2)≃{[ sin(πα)/πB(t_1/t_2;α,1-α),                      0<α<1;     (t_2-t_1)^1-α-t_2^1-α,                      1<α<2 ].

for t_1<t_2. The results in Eq. (<ref>) are the asymptotic forms for large difference between t_1 and t_2.
More detailed messages of the quantities of the stochastic processes are commonly determined by the multi-point PDFs of the inverse subordinator s(t), especially the two-point one in Eq. (<ref>). In the following sections, we will establish the Langevin picture to describe the Lévy walk in the presence of a constant force and evaluate some concerned statistical quantities through the two- and four-point PDFs of inverse subordinator.





§ LÉVY WALK IN A CONSTANT FORCE FIELD


The Lévy walk model with constant force F can be described as a set of Langevin equations coupled with a subordinator as

    d/d tx(t)   =v(t),
    d/d sv(s)   =-γ v(s) +F η(s)+ξ(s),
    d/d st(s)   = η(s),

where γ is the friction coefficient, ξ(s) is a Gaussian white noise with null mean value ⟨ξ(s)⟩=0 and ⟨ξ(s_1)ξ(s_2)⟩=2Dδ(s_1-s_2). The Lévy noise η(s), regarded as the formal derivative of the α-dependent subordinator t(s), is independent with the Gaussian white noise ξ(s). The derivative of position x with respect to physical time t is velocity v and the subordinator t(s) is aimed to characterize the distribution of duration of each flight of Lévy walk. The initial position and velocity are assumed to be x(0)=v(0)=0. When F=0, the Langevin picture (<ref>) is coincided with the free-force case <cit.>.

The constant force F multiplied by Lévy noise η(s) is meant to affect the stochastic process for all physical time t after making the subordination <cit.>. Otherwise, the force F is only effective over operational time s, similar to the case that the external force only affects on the instant of jump in CTRW model, and it is invalid during the trap event with the constant s <cit.>.

To be more intuitive, the Langevin equation for the velocity process in Eq. (<ref>) can be transformed to the one in physical time t by using the identity v(t):=v(s(t)), i.e.,

    d/d tv(t)=-γ v(t) d/d ts(t)+F+ξ(s(t))d/d ts(t),

where ξ(s(t)) can be seen as the time-changed Gaussian white noise <cit.>. Equation (<ref>) implies that the external force F affects the velocity process for all physical time t. Especially when the inverse subordinator s(t) remains a constant, this equation reduces to dv(t)/dt=F, characterizing the fact that the constant force provides an acceleration F to the unidirectional motion of Lévy walk.

On the other hand, the velocity process v(s) can be solved from the second equation of Eq. (<ref>), by using the Laplace transform technique, i.e.,

    v(s)=F∫_0^s e^-γ(s-s')η(s')ds'+∫_0^s e^-γ(s-s')ξ(s')ds'.

Equation (<ref>) shows that the velocity v(s) is contributed by two parts; one comes from the external force F and another one from the random force ξ(s) which corresponds to the free Lévy walk. Due to the independence of noises η(s) and ξ(s), the two terms in Eq. (<ref>) are independent, and the subordinated process v(t):=v(s(t)) in physical time t is given by

    v(t)   =F∫_0^t e^-γ(s(t)-s(t'))dt'  
          +∫_0^t e^-γ(s(t)-s(t'))ξ(s(t'))ds(t').

See Appendix <ref> for detailed derivations of Eqs. (<ref>) and (<ref>). Equations (<ref>) and (<ref>) are complementary for further quantitative calculations in the following sections. Here, the first term of velocity in Eq. (<ref>) is resulted from the contribution of the constant force F. Instead of the naive expectation that this contribution is Ft during time t, it is found to be not larger than Ft due to the non-decreasing property of s(t).
The expression of v(t) in Eq. (<ref>) means that the evolution of velocity process starts from time zero, and the force's contribution reaches the maximum Ft only if s(t')≡ 0 when t'∈[0,t], which implies no renewal happens during the measurement time t.

Once the renewal happens, the system will lose some of the extra momentums provided by the constant force, and then the net increased velocity is less than Ft.

In particular, let us focus on the value of velocity between any two successive renewal points.
As the graph of one sample of inverse subordinator s(t) in Fig. <ref> shows, if no renewal happens between the renewal points t_1 and t_2, the s(t) will remain a constant for t∈[t_1,t_2]. By using this property of s(t), we obtain the increment of velocity process in Eq. (<ref>) between two successive renewal points

    v(t_2)-v(t_1)   =F∫_t_1^t_2 e^-γ(s(t_1)-s(t'))dt'
       =F(t_2-t_1),

where the second term in Eq. (<ref>) vanishes for constant s(t). Equation (<ref>) implies that the acceleration in the no-renewal period is constant F.



An alternative description of Lévy walk in a constant force that we have to mention is the collision model in Ref. <cit.>. It says that the change of the velocity of a free test particle is due to the collision with the surrounding bath particles. If a constant force F is added, the acceleration is F between two successive collisions. However, some of the added momentum in the history will be transmitted to the surrounding bath particles when the collision happens according to the energy and momentum conservation laws.

Now, the coupled Langevin equation in Eq. (<ref>) is a continuous limit of the collision model in some sense. The Gaussian white noise ξ(s) in operation time s describes the random force coming from the collision with the surrounding small bath particles, and the subordinator t(s) corresponds to the random time distributions between two successive collisions by taking s=1. Besides the intuitive analyses, we will show that the coupled Langevin system (<ref>) and the collision model are quantitatively consistent in some important quantities—the means of velocity v(t) and displacement x(t).
 
 Furthermore, based on the
explicit expressions of velocity process in Eqs. (<ref>) and (<ref>), the correlation function of velocity, EAMSD and TAMSD, which have not been obtained before, can also be evaluated.







§ FIRST MOMENT

This section is aimed to derive the mean value of v(t) and x(t) in Eq. (<ref>) by using the velocity process in Eqs. (<ref>) and (<ref>), together with the subordination method <cit.>.
The mean value of velocity process in physical time can be obtained by directly taking the average on Eq. (<ref>),

    ⟨ v(t)⟩=F∫_0^t ⟨ e^-γ(s(t)-s(t'))⟩ dt'.

The ensemble average of the second term in Eq. (<ref>) vanishes due to the zero mean value of ξ(s).  It is interesting to find that the first moment of the velocity process in Eq. (<ref>) depends on the two-point PDF of the inverse subordinator.
By using the two-point PDF h(s,t;s',t') in (<ref>) and the technique of Laplace transform (t→λ,t'→λ'), the integrand in Eq. (<ref>) can be obtained, for small λ and λ', as

    ℒ [⟨ e^-γ s(t)e^γ s(t')⟩]
       =∫_0^∞∫_0^∞ e^-γ (s-s')h(s,λ;s',λ')ds'ds
       ≃Φ(λ)+Φ(λ')-Φ(λ+λ')/λλ'Φ(λ+λ'),

which is the same as p_0(λ',λ) in Eq. (<ref>).
The rationality of small λ and λ' is that we pay attention to the large-t behavior of this system and the integrand of Eq. (<ref>) is dominated by the part with large t'. Equation (<ref>) implies that the integrand of Eq. (<ref>) is asymptotically equal to p_0(t',t) with t'<t.
This result indicates that the first term of two-point PDF h(s_1,t_1;s_2,t_2) in Eq. (<ref>) containing δ-function  plays a leading role when integrated together with an exponential kernel e^-γ(s_1-s_2) in Eq. (<ref>).

Substituting the expression p_0(t',t) in Eq. (<ref>) and calculating the integral in Eq. (<ref>) yield

    ⟨ v(t)⟩≃{[       F(1-α)t,         0<α<1,; Fα-1/2-αt^2-α,         1<α<2. ].

The mean value of x(t) affected by the constant force F is the integral of ⟨ v(t)⟩, i.e.,

    ⟨ x(t)⟩_F ≃{[            F(1-α)/2t^2,                  0<α<1,; F(α-1)/(2-α)(3-α)t^3-α,                  1<α<2. ].

Considering the EAMSD ⟨ x^2(t)⟩_0 of free Lévy walk <cit.> and the first moment in presence of the constant force F in Eq. (<ref>), the generalized Einstein relation <cit.> is satisfied by Lévy walk <cit.>

    ⟨ x(t)⟩_F = ⟨ x^2(t)⟩_0/2k_B𝒯F,

where the effective kinetic temperature k_B𝒯 is equal to D/γ for the Langevin system in Eq. (<ref>). However, we find that the generalized Einstein relation for time averages does not hold, i.e.,

    ⟨δ^1(Δ)⟩_F ≠⟨δ^2(Δ)⟩_0/2k_B 𝒯F

with ⟨δ^1(Δ)⟩_F=1/T-Δ∫_0^T-Δ⟨ x(t+Δ)-x(t)⟩_F dt and  ⟨δ^2(Δ)⟩_0 being the TAMSD of free Lévy walk. Actually, for Lévy walk, there is <cit.>

    ⟨δ^1(Δ)⟩_F/⟨δ^2(Δ)⟩_0=F|α-1|/2k_B𝒯(T/Δ)^α̃-1,

which is consistent with the result in Ref. <cit.>. Here, α̃=2 when 0<α<1 and α̃=3-α when 1<α<2.




§ VELOCITY CORRELATION FUNCTION

Both the EAMSD and TAMSD of x(t) depend on the velocity correlation function ⟨ v(t_1)v(t_2)⟩, which is much more difficult to calculate than the mean value in Section <ref>. Considering the independence of the two terms in Eq. (<ref>), the velocity correlation function only consists of two parts:

    ⟨ v(t_1)v(t_2)⟩=⟨ v(t_1)v(t_2)⟩_1+⟨ v(t_1)v(t_2)⟩_2.

The first term comes from the F-dependent part in Eq. (<ref>),

    ⟨ v(t_1)v(t_2)⟩_1= F^2∫_0^t_1∫_0^t_2⟨ W(t_1,t_1',t_2,t_2') ⟩ dt_2'dt_1'.

We use the symbol

    W(t_1,t_1',t_2,t_2')=e^-γ s(t_1)e^γ s(t_1')e^-γ s(t_2)e^γ s(t_2')

to denote the double exponential kernel, the ensemble-average of which depends on the four-point PDF of inverse subordinator s(t).
As for the second term in Eq. (<ref>), it is the same as the free-force case and is more convenient to be derived in operation time from Eq. (<ref>) and then apply the subordination method, i.e.,

    ⟨ v(t_1)v(t_2)⟩_2
    
        =∫_0^∞∫_0^∞⟨ v(s_1)v(s_2)⟩_2 h(s_1, t_1; s_2, t_2)ds_1ds_2,

where the velocity correlation function ⟨ v(s_1)v(s_2)⟩_2 in operation time can be obtained through the second part of v(s) in Eq. (<ref>)

    ⟨ v(s_1)v(s_2)⟩_2= D/γ (e^-γ |s_1- s_2|-e^-γ(s_1+s_2)).

Then by use of the two-point PDF of the inverse subordinator in Eq. (<ref>), we obtain

    ⟨ v(t_1)v(t_2)⟩_2 ≃D/γp_0(t_1,t_2)

for large t_1, t_2 and t_1<t_2, which is the same as the velocity correlation function of free Lévy walk in CTRW framework when the constant D/γ=v_0^2 <cit.>.

The remaining question is how to evaluate ⟨ W(t_1,t_1',t_2,t_2') ⟩ in the first part of velocity correlation function in Eq. (<ref>), which looks quite complex since it depends on the four-point PDF of inverse subordinator s(t). However, based on the
analyses of Eq. (<ref>), we find that the leading term of four-point PDF should contain δ(t_1-t_1')δ(t_2-t_2') due to the exponential kernel W, and the value of ⟨ W(t_1,t_1',t_2,t_2') ⟩ might be closely related to p_0(t_1',t_1) and p_0(t_2',t_2).

More precisely, the value of ⟨ W(t_1,t_1',t_2,t_2') ⟩ depends on the magnitude relation between the four time points t_1,t_1',t_2,t_2'. Since t_1'<t_1 and t_2'<t_2 in Eq. (<ref>) and it is assumed that t_1<t_2, there are totally three different cases for different range of t_2'. 
For the first two cases, there are

    ⟨ W(t_1,t_1',t_2,t_2') ⟩≃{[     p_0(t_2',t_2), t_2'<t_1'<t_1<t_2,;     p_0(t_1',t_2), t_1'<t_2'<t_1<t_2. ].

The situation becomes more complex for the case t_1'<t_1<t_2'<t_2. In this case, ⟨ W(t_1,t_1',t_2,t_2') ⟩ consists of two parts; one is p_0(t_1',t_2) denoting no renewal happens throughout [t_1',t_2] and another one denoting some renewals happen in the middle interval [t_1,t_2']. The former one is easy to be calculated as those in Eq. (<ref>),  while the latter one depends on the four-point PDF of inverse subordinator s(t) and it should be contributed from the term containing δ(t_2-t_2')Θ(t_2'-t_1)δ(t_1-t_1') in the four-point PDF of inverse subordinator. With some detailed derivations in Appendix <ref>, we obtain the result for t_1'<t_1<t_2'<t_2,

    ⟨ W(t_1,t_1',t_2,t_2') ⟩≃{[ p_0(t_1',t_1)p_0(t_2',t_2),                      0<α<1,;              p_0(t_1',t_2),                      1<α<2. ].


For the case t_1<t_2, splitting the integral Eq. (<ref>) into three parts according to the range of t_2' as

    ∫_0^t_1dt_1' ∫_0^t_2dt_2' =   ∫_0^t_1dt_1' ∫_0^t_1'dt_2' +∫_0^t_1dt_1' ∫_t_1'^t_1dt_2' 
       +∫_0^t_1dt_1' ∫_t_1^t_2dt_2'

and using the corresponding forms of ⟨ W(t_1,t_1',t_2,t_2') ⟩, we obtain, for large t_1 and t_2,


    ⟨ v(t_1)v(t_2)⟩_1 ≃{[ F^2sin(πα)/π[ t_2^2B(t_1/t_2;α+2,1-α) + α t_1^2B(t_1/t_2;α,1-α) .;         . -(1+α)t_1t_2B(t_1/t_2;α+1,1-α) ]  + F^2 (1-α)^2 t_1t_2,                                                            0<α<1,;         F^2(α-1)/2-α(t_1t_2^2-α-1/3-αt_2^3-α+1/3-α(t_2-t_1)^3-α),                                                            1<α<2.;                                                                   ].







§ EAMSD AND TAMSD

Once the velocity correlation function is obtained, we can rewrite the summation of Eqs. (<ref>) and (<ref>) as the scaling form with two separate parts 

    ⟨ v(t_1)v(t_2)⟩≃ C_1t^ν_1-2ϕ_1(τ/t)+C_2t^ν_2-2ϕ_2(τ/t)

for large t and τ. See Appendix <ref> for the expressions of C_i, ν_i, ϕ_i (i=1,2). We use two different scaling forms since the velocity correlation function contains two independent parts in Eq. (<ref>).
The dependence of EAMSD and TAMSD on the scaled velocity correlation function has been revealed through the generalized Green-Kubo formula <cit.>.



Here what we should pay attention to is to use the generalized Green-Kubo formula to handle the two parts in Eq. (<ref>) individually and then combine them together. After the lengthy calculations in Appendix <ref>, we obtain the EAMSD of Lévy walk in the presence of a constant force as

    ⟨(Δ x(t))^2⟩   =⟨ x^2(t)⟩-⟨ x(t)⟩^2  
       ≃{[ F^2α(1-α)(2-α)(3+α)/72t^4,                     0<α<1,; F^2 (α-1)/(4-α)(5-α)t^5-α,                     1<α<2,;                            ].

both of which present the super-ballistic diffusion behavior resulting from the acceleration. For 0<α<1, the frequency of the collision is too low so that the ensemble-averaged acceleration keeps proportional to F in Eq. (<ref>). While for 1<α<2, the collision happens more frequently then the ensemble-averaged acceleration decays as t^1-α in Eq. (<ref>) and the diffusion behavior is slower than t^4.
The corresponding TAMSD is

    ⟨δ^2(Δ)⟩   ≃{[            F^2(1-α)α/6T^2Δ^2,                        0<α<1,; F^2 (α-1)/(3-α)(4-α)T^3-αΔ^2,                        1<α<2,;                               ].

for large T and Δ≪ T. Different from the free Lévy walk, both of the TAMSDs here are increasing as ∝Δ^2 whatever 0<α<1 or 1<α<2. This system is weakly nonergodic since the EAMSD and TAMSD are different in the exponents of t and Δ. In Fig. <ref>, the simulations of EAMSDs and TAMSDs for different α agree with the theoretical results very well.






§ SUMMARY

Lévy walk is a typical and practical model and it has various applications in the natural world.
In this paper, we mainly investigate the anomalous and nonergodic behavior of the Lévy walk under the effects of a constant force in physical times. The Langevin picture of such a model has been established and compared with a previously proposed collision model in <cit.>.
The consistency between the two models has been illustrated/verified through the intuitive descriptions of the models and the quantitative comparison of the first moment of velocity process.

Compared with the collision model, one of the advantages of Langevin equation is that the velocity correlation function can be conveniently obtained by using the techniques of subordination. As for the method of subordination, we get some new findings in this paper. The first one is the relationship between the PDF of inverse subordinator and the PDF of the number of renewals. The second one is that by adding the term Fη(s) in the Langevin equation on operation time s (in Eq. (<ref>)), we obtain an effective force being active for all physical times. With some changes of the formulation of the Langevin equation, we clearly show the effects of the external force on the Langevin system and connect it with the collision model. The third one is an interesting and rare phenomenon that the first moments of the velocity and displacement depend on the two-point PDF of the inverse subordinator. More complicated, the velocity correlation function depends on the four-point PDF of the inverse subordinator. By using the relationship between the method of subordinator and the renewal theory, we simplify the derivations and obtain the velocity correlation function.

After getting the velocity correlation function, we can obtain the EAMSD and TAMSD by using the generalized Green-Kubo formula. The super-ballistic diffusion and weakly non-ergodic behavior are observed for this Langevin system. Apart from the interesting phenomena for Lévy walk in the presence of a constant force, we are also interested in the Lévy walk under the effects of other kinds of external forces, such as a harmonic potential or a time-dependent force. How to build an effective Langevin equation to describe these forces and to analyze some important statistical quantities will be investigated in our future work.






§ ACKNOWLEDGMENTS

This work was supported by the National Natural Science Foundation of China under grant no. 11671182, and the Fundamental Research Funds for the Central Universities under grants no. lzujbky-2018-ot03 and no. lzujbky-2019-it17.




§ DERIVATIONS OF EQS. (<REF>) AND (<REF>)

Let us first derive Eq. (<ref>), the Langevin equation of velocity process v(t) evolving in physical time t. Replacing s by s(t) in the second equation of Eq. (<ref>), one arrives at

    d/d s(t)v(s(t))=-γ v(s(t))+F η(s(t))+ξ(s(t)).

It is equivalent to the equation

    d/dtv(t)   =-γ v(s(t))d/dts(t)+F η(s(t))d/dts(t)+ξ(s(t))d/dts(t)
       =-γ v(t)d/dts(t)+F+ξ(s(t))d/dts(t),

where we have used the identity v(t):=v(s(t)) and

    η(s(t))=dt(s(t))/ds(t)= dt/ds(t)

in the last line.

As for Eq. (<ref>), replacing s in Eq. (<ref>) with s(t) gives

    v(t):   =v(s(t))
       =F∫_0^s(t) e^-γ(s(t)-s')η(s')ds'+∫_0^s(t) e^-γ(s(t)-s')ξ(s')ds'
       =F∫_0^t e^-γ(s(t)-s(t'))η(s(t'))ds(t')+∫_0^t e^-γ(s(t)-s(t'))ξ(s(t'))ds(t')
       =F∫_0^t e^-γ(s(t)-s(t'))dt'+∫_0^t e^-γ(s(t)-s(t'))ξ(s(t'))ds(t'),

where we use the variable substitution s'→ s(t') in the third line.






§ DERIVATIONS OF EQ. (<REF>)

Let us consider the symbol in Eq. (<ref>):

    W(t_1,t_2,t_3,t_4)= e^-γ s(t_1)e^γ s(t_2)e^-γ s(t_3)e^γ s(t_4)

for t_4<t_3<t_2<t_1. As Sec. <ref> discusses, the dominating role in ⟨ W(t_1,t_2,t_3,t_4)⟩ comes from two parts; one is p_0(t_4,t_1) denoting no renewal happens during [t_4,t_1] and another one denoting some renewals happen in the middle interval [t_3,t_2]. The latter one depends on the four-point PDF h(·) containing δ(s_1-s_2)Θ(s_2-s_3)δ(s_3-s_4) of inverse subordinator s(t).
Similarly to the two-point case, the expression ⟨ W(t_1,t_2,t_3,t_4)⟩ can be calculated through Laplace transform (t_1→λ_1,t_2→λ_2,t_3→λ_3,t_4→λ_4)

    ⟨ W(λ_1,λ_2,λ_3,λ_4)⟩   =∫_0^∞∫_0^∞∫_0^∞∫_0^∞  e^-γ s_1e^γ s_2e^-γ s_3e^γ s_4h(s_1,λ_1;s_2,λ_2;s_3,λ_3;s_4,λ_4)ds_1ds_2ds_3ds_4.

The key in Eq. (<ref>) is the four-point PDF h(·) of inverse subordinator s(t), which is related to the corresponding four-point PDF Z(·) of subordinator t(s) as

    h(s_1,λ_1;s_2,λ_2;s_3,λ_3;s_4,λ_4) = ∂/∂ s_1∂/∂ s_2∂/∂ s_3∂/∂ s_41/λ_1λ_2λ_3λ_4 Z(λ_1,s_1;λ_2,s_2;λ_3,s_3;λ_4,s_4).

Considering the independence and stationarity of the increments of subordinator t(s), the four-point PDF Z(·) can be obtained as

    Z(λ_1,s_1;λ_2,s_2;λ_3,s_3;λ_4,s_4) = e^-s_4Φ(λ_1+λ_2+λ_3+λ_4)e^-(s_3-s_4)Φ(λ_1+λ_2+λ_3)
      e^-(s_2-s_3)Φ(λ_1+λ_2)e^-(s_1-s_2)Φ(λ_1)

under the condition s_4<s_3<s_2<s_1. There are twenty-four different orders for s_1,⋯,s_4, which makes the exact expression of h(s_1,λ_1;s_2,λ_2;s_3,λ_3;s_4,λ_4) quite complicated. Fortunately, here we only need to know the term containing δ(s_1-s_2)Θ(s_2-s_3)δ(s_3-s_4) in h(·) which denotes some renewals happening in the middle interval [t_3,t_2]. This term can be obtained from totally four kinds of terms in Z(·), which are

    Θ(s_1-s_2)Θ(s_2-s_3)Θ(s_3-s_4)e^-s_4Φ(λ_1+λ_2+λ_3+λ_4)e^-(s_3-s_4)Φ(λ_1+λ_2+λ_3)
      e^-(s_2-s_3)Φ(λ_1+λ_2)e^-(s_1-s_2)Φ(λ_1)
    
    +   Θ(s_1-s_2)Θ(s_2-s_4)Θ(s_4-s_3)e^-s_3Φ(λ_1+λ_2+λ_3+λ_4)e^-(s_4-s_3)Φ(λ_1+λ_2+λ_4)
      e^-(s_2-s_4)Φ(λ_1+λ_2)e^-(s_1-s_2)Φ(λ_1)
    
    +   Θ(s_2-s_1)Θ(s_1-s_3)Θ(s_3-s_4)e^-s_4Φ(λ_1+λ_2+λ_3+λ_4)e^-(s_3-s_4)Φ(λ_1+λ_2+λ_3)
      e^-(s_1-s_3)Φ(λ_1+λ_2)e^-(s_2-s_1)Φ(λ_2)
    
    +   Θ(s_2-s_1)Θ(s_1-s_4)Θ(s_4-s_3)e^-s_3Φ(λ_1+λ_2+λ_3+λ_4)e^-(s_4-s_3)Φ(λ_1+λ_2+λ_4)
      e^-(s_1-s_4)Φ(λ_1+λ_2)e^-(s_2-s_1)Φ(λ_2).

Taking derivatives with respect to s_4,s_3,s_2,s_1 on these four terms, together with some technical calculations, we find the factor multiplied by δ(s_1-s_2)Θ(s_2-s_3)δ(s_3-s_4) in h(·) is

    [Φ(λ_1)+Φ(λ_2)-Φ(λ_1+λ_2)]· [Φ(λ_1+λ_2+λ_3)+Φ(λ_1+λ_2+λ_4)-Φ(λ_1+λ_2)-Φ(λ_1+λ_2+λ_3+λ_4)]/λ_1λ_2λ_3λ_4  e^(s_2-s_3)Φ(λ_1+λ_2)e^s_3Φ(λ_1+λ_2+λ_3+λ_4).

Substituting it into (<ref>) gives

    [Φ(λ_1)+Φ(λ_2)-Φ(λ_1+λ_2)]· [Φ(λ_1+λ_2+λ_3)+Φ(λ_1+λ_2+λ_4)-Φ(λ_1+λ_2)-Φ(λ_1+λ_2+λ_3+λ_4)]/λ_1λ_2λ_3λ_4 Φ(λ_1+λ_2)Φ(λ_1+λ_2+λ_3+λ_4).

Performing inverse Laplace transform of (<ref>), we can obtain the term of ⟨ W(t_1,t_2,t_3,t_4) ⟩ which represents some renewals happening in the middle interval [t_3,t_2].
For 1<α<2, it can be seen that the term (<ref>) is far less than another term p_0(t_4,t_1) in ⟨ W(t_1,t_2,t_3,t_4) ⟩ for small λ_i, i=1,⋯,4, while for 0<α<1 the two terms are of the same order and the sum is

    ⟨ W(t_1,t_2,t_3,t_4) ⟩≃ p_0(t_4,t_3)p_0(t_2,t_1),

which is the result in Eq. (<ref>) after changing the symbol of time (t_4⟶ t_1',t_3⟶ t_1,t_2⟶ t_2',t_1⟶ t_2).











§ DERIVATIONS OF EAMSD IN EQ. (<REF>) AND TAMSD IN EQ. (<REF>)

According to the generalized Green-Kubo formula <cit.>, one should first rewrite the velocity correlation function in Eqs. (<ref>) and (<ref>) into the scaling form with two separate parts to obtain the EAMSD and TAMSD, i.e.,

    ⟨ v(t)v(t+τ)⟩   =⟨ v(t)v(t+τ)⟩_1+⟨ v(t)v(t+τ)⟩_2
       ≃ C_1t^ν_1-2ϕ_1(τ/t)+C_2t^ν_2-2ϕ_2(τ/t);

and we assume that ϕ_1(q)→ c_1q^-δ_1 and ϕ_2(q)→ c_2q^-δ_2 as q→ 0.

For the case 0<α<1, the velocity correlation function can be written as

    ⟨ v(t)v(t+τ)⟩
    =F^2/Γ(1-α)Γ(α)t^2ϕ_1(τ/t)+D/γ1/Γ(1-α)Γ(α)ϕ_2(τ/t),

where

    ϕ_1(q)    = (1+q)^2 B(1/1+q;α+2,1-α) + α B(1/1+q;α,1-α) 
           - (1+α)(1+q)B(1/1+q;α+1,1-α) + Γ(α)Γ(1-α)(1-α)^2(1+q)

and

    ϕ_2(q)=B(1/1+q;α,1-α).

When q→0,  ϕ_1(q) converges to the constant c_1=Γ(α)Γ(3-α)/2 and ϕ_2(q) to c_2=Γ(α)Γ(1-α). Then we can obtain the parameters needed: C_1=F^2/Γ(1-α)Γ(α), ν_1=4, C_2=D/γ1/Γ(1-α)Γ(α), and ν_2=2.



After getting the necessary parameters, according to the generalized Green-Kubo formula, the second moment of x(t) can be obtained as

    ⟨ x^2(t)⟩=⟨ x^2(t)⟩_1+⟨ x^2(t)⟩_2=2K_1t^ν_1+2K_2t^ν_2

with ⟨ x^2(t)⟩_i=∫_0^t∫_0^t⟨ v(t_1')v(t_2')⟩_idt_1' dt_2' and K_i=C_i/ν_i∫_0^∞ (1+q)^-ν_iϕ_i(q)dq, i=1,2. Further  substituting the expressions of C_i, ν_i and ϕ_i into it, one has
K_1=F^2(α(1-α)(2-α)(3+α)/144+(1-α)^2/8) and K_2=D(1-α)/2γ.
Then the second moment of the Lévy walk with constant force

    ⟨ x^2(t)⟩=2K_1t^4+D(1-α)/γt^2

and the EAMSD

    ⟨ x^2(t)⟩-⟨ x(t)⟩^2
    ≃F^2α(1-α)(2-α)(3+α)/72t^4

for large time t.

As for the TAMSD of the Lévy walk with constant force, we split it into three parts

    ⟨δ^2(Δ)⟩   =1/T-Δ∫_0^T-Δ⟨ (x(t+Δ)-x(t))^2⟩ dt -1/T-Δ∫_0^T-Δ⟨ x(t+Δ)-x(t)⟩^2dt
       =⟨δ^2(Δ)⟩_1+⟨δ^2(Δ)⟩_2-⟨δ^2(Δ)⟩_3,

where

    ⟨δ^2(Δ)⟩_i=1/T-Δ∫_0^T-Δ⟨ (x(t+Δ)-x(t))^2⟩_idt

and

    ⟨δ^2(Δ)⟩_3=1/T-Δ∫_0^T-Δ⟨ x(t+Δ)-x(t)⟩^2dt.

By virtue of the generalized Green-Kubo formula, the first two terms can be obtained as

    ⟨δ^2(Δ)⟩_i
    =2c_iC_i/(β_i+1)(ν_i-β_i-1)(ν_i-β_i)T^β_iΔ^ν_i-β_i,

where β_i is the exponent of second moment of velocity ⟨ v^2(t)⟩_i∝ t^β_i with β_1=2 and β_2=0.
After substituting the parameters into Eq. (<ref>), one has

    ⟨δ^2(Δ)⟩_1=F^2(1-α)(2-α)/6T^2Δ^2,     ⟨δ^2(Δ)⟩_2=D/γΔ^2.

Besides,  substituting the first moment of position ⟨ x(t)⟩_F= F(1-α)/2t^2 into Eq. (<ref>) leads to

    ⟨δ^2(Δ)⟩_3=F^2(1-α)^2/3T^2Δ^2.

Finally, the TAMSD for large T and Δ≪ T is
⟨δ^2(Δ)⟩≃F^2(1-α)α/6T^2Δ^2.



For the case 1<α<2, the velocity correlation function is

    ⟨ v(t)v(t+τ)⟩   =⟨ v(t)v(t+τ)⟩_1+⟨ v(t)v(t+τ)⟩_2
       =F^2(α-1)/2-αt^3-αϕ_1(τ/t)+D/γt^1-αϕ_2(τ/t),

where

    ϕ_1(q)= (1+q)^2-α-1/3-α(1+q)^3-α+1/3-αq^3-α

and

    ϕ_2(q)=q^1-α-(1+q)^1-α.

When q→0, we have ϕ_1(q)→2-α/3-α and ϕ_2(q)→ q^1-α. Then we obtain the parameters needed: C_1=F^2(α-1)/2-α, ν_1=5-α, c_1=2-α/3-α, β_1=3-α, K_1=F^2(α-1)/2(4-α)(5-α), C_2=D/γ, ν_2=3-α, c_2=1, β_2=0, and K_2=D(α-1)/γ(2-α)(3-α).
Finally, the second moment and the EAMSD of x(t) are

    ⟨ x^2(t)⟩=F^2(α-1)/(4-α)(5-α)t^5-α+2D(α-1)/γ(2-α)(3-α)t^3-α

and

    ⟨ x^2(t)⟩-⟨ x(t)⟩^2
    ≃F^2(α-1)/(4-α)(5-α)t^5-α

for large time t.
Similarly, the TAMSD for large T and Δ≪ T is

    ⟨δ^2(Δ)⟩   =F^2(α-1)/(3-α)(4-α)T^3-αΔ^2+2D/γ (2-α)(3-α)Δ^3-α-F^2(α-1)^2/(2-α)^2(5-2α)T^4-2αΔ^2
       ≃F^2(α-1)/(3-α)(4-α)T^3-αΔ^2.



















































§ REFERENCES

apsrev4-1


