L\'{e}vy-walk-like Langevin dynamics

Continuous time random walks and Langevin equations are two classes of stochastic models for describing the dynamics of particles in the natural world. While some of the processes can be conveniently characterized by both of them, more often one model has significant advantages (or has to be used) compared with the other one. In this paper, we consider the weakly damped Langevin system coupled with a new subordinator|$\alpha$-dependent subordinator with $1<\alpha<2$. We pay attention to the diffusion behaviour of the stochastic process described by this coupled Langevin system, and find the super-ballistic diffusion phenomena for the system with an unconfined potential on velocity but sub-ballistic superdiffusion phenomenon with a confined potential, which is like L\'{e}vy walk for long times. One can further note that the two-point distribution of inverse subordinator affects mean square displacement of this coupled weakly damped Langevin system in essential.


Introduction
Introduced just 100 years ago for describing Brownian motion [1], the Langevin equation, nowadays, has more wide applications and plays a central role in modeling the dynamical systems coupled with a fluctuating environment [2]. One of the main advantages of this equation is that it builds a relation between physically transparent and mathematically tractable description for a complex stochastic dynamical system. Assuming that a particle moves in a fluid without friction, it receives a blow due to a random collision with a molecule, then the velocity of the particle changes. This procedure can be well modeled by weakly damped Langevin system. However, if the fluid is very viscous, the change of velocity is quickly dissipated and the net result of an impact is a change in the displacement of the particle. The overdamped Langevin system becomes more suitable to describe the motion of particles in this case.
Another kind of popular microscopic model is continuous time random walk (CTRW), originally introduced by Montroll and Weiss in 1965 [3]. It is a powerful mathematical framework to model complex dynamical behaviors, especially anomalous diffusion phenomena characterized by nonlinear time dependence of mean squared displacement (MSD); see the reviews [4,5,6] and references therein. In the CTRW framework, the motion of particles is described through consecutive waiting times and jumps between them. The two random variables, waiting times and jump lengths, are drawn from some associated distributions, which could decide the diffusive behaviour of the particles. Fogedby [7] proposed that an overdamped Langevin equation in operation time s coupled with a physical time process t(s) (named as a subordinator) can model the same process as CTRWs in scaling limits. There is an advantage in the Langevin system that the external force field can be included naturally with clear physical meaning, and the system is given as d ds where the position x is penalized by the operation time s. When f (x) = 0 and ξ(s) is Gaussian white noise, model (1) yields subdiffusion. Nowadays, subordinator is a very effective tool to characterize various complex dynamical systems, especially some reallife data in biology [8], financial time series [9], ecology [10], and physics [11]. Besides, overdamped Langevin equation together with its generalizations are well-developed in recent years, e.g., the heterogeneous diffusion processes [12,13,14] and Brownian yet non-Gaussian diffusion [15,16]. At the same time, there are also many processes in practice which could not be well characterized by model (1), since the particles may be weakly damped, e.g., cold atoms diffusing in optical lattices [17,18], and the class of viscoelastic diffusion described by the generalized Langevin equation with (tempered) power-law friction kernel [19,20,21,22] and of (tempered) fractional Brownian motion [23,24,25]. Another main class of generalized weakly damped Langevin equations are coupled with a subordinator not a friction kernel. Eule et al. [26] presented three kinds of weakly damped Langevin system coupled with the α 0 -stable subordinator (0 < α 0 < 1), which are related to three kinds of fractional Klein-Kramers equations [27,28,29], respectively. In this paper, we also consider the weakly damped Langevin system, but extend the subordinator to be αdependent with 1 < α < 2. To the best of our knowledge, the subordinator with 1 < α < 2 has never been considered in Langevin system before. One possible difficulty is that the original α-stable Lévy process with 1 < α < 2 is not a non-decreasing random process while the one-sided α 0 -stable with 0 < α 0 < 1 is. The condition of nondecreasing must be guaranteed from a physical point of view since the subordinator t(s) denotes the waiting time process in CTRWs [30]. Fortunately, through Lévy-Khinchin representation [31], an appropriate subordinator can be designed by specifying a Lévy measure. The method of generating this subordinator for numerical simulations is given in the last part.
Based on the designed α-dependent subordinator with 1 < α < 2, we mainly discuss the diffusive behaviour of the weakly damped Langevin system coupled with this subordinator and with two different potentials, i.e., unconfined and confined ones on velocity. The harmonic potential U(v) = γv 2 /2 is chosen with γ = 0 and γ = 0, respectively, corresponding to the unconfined and confined case. For long times, in the former case, the particles spread like Richardson-Obukhov diffusion in turbulence, where the velocity follows a simple Brownian motion [32,33]. In the latter case, the particle motion is like Lévy walk [34,35,36,37] in the sub-ballistic superdiffusion regime. In this way, the mathematical description of Lévy walk confined to an external force field could be constructed naturally. It is discovered that the essential difference made by the new subordinator on the diffusive behavior comes from the two-point probability density function (PDF) of inverse subordinator.
The remainder of this paper is organized as follows. In section 2, we define the α-dependent subordinator (1 < α < 2), and discuss its properties as well as its corresponding inverse subordinator. In section 3, we present the weakly damped Langevin system coupled with this subordinator, derive its corresponding fractional Klein-Kramers equation, and explicitly investigate the diffusive behavior of the stochastic process described by this Langevin system for two cases of friction factor γ = 0 and γ = 0. Then we discuss the relations between CTRWs and Langevin system with different subordinators; the Langevin system with α-dependent subordinator (1 < α < 2) is presented in section 4 and another kind of Langevin system in section 5. We give the method of generating α-dependent subordinator with 1 < α < 2 and thus the subordinated processes for numerical simulations in section 6. A summary of the key results is made in section 7. In the appendices some mathematical details are collected.
We call the pair (b, ν) the characteristics of the subordinator t(s). If it is taken to be b = 0 and with 0 < α < 1, then Φ(λ) = λ α is the Laplace exponent of the one-sided α-dependent subordinator for 0 < α < 1, which has been fully discussed in Langevin systems [7,26,38,39]. Here, we would like to specify the pair (b, ν) to form a subordinator for 1 < α < 2. For this purpose, considering the requirements of ν, we take b = 0 and where τ 0 is the characteristic time. Then its Laplace exponent reads The asymptotic behavior λ → 0 in Laplace space corresponds to t → ∞ in the time domain, which is that people always pay attention to in physical experiments. The two-point PDF of the subordinator t(s) can be expressed as Its corresponding Laplace transform can be directly derived due to the Markovian character of this process. If s 1 < s 2 , considering the stationary and independent increments of the Lévy process, the characteristic function of g(t 2 , s 2 ; t 1 , s 1 ) is [40] g(λ 2 , s 2 ; As for general s 1 and s 2 , we usually write the characteristic function aŝ where Θ(x) denotes the Heaviside step function: Θ(x) = 1 for x > 0, Θ(x) = 0 for x < 0, and Θ(x = 0) = 1/2. By the same approach, n-point joint PDF for the subordinator t(s) can also be obtained. The multiple-point PDFs of (inverse) subordinator have been considered in the pioneering works [40,41]. Here, we just extend the results to more general subordinator, e.g., the α-dependent subordinator t(s) with 1 < α < 2.
The first-passage time of a subordinator {t(s), s ≥ 0} is called inverse subordinator {s(t), t ≥ 0} [42,43], defined as Denote the (multiple-point) PDF of inverse subordinator s(t) as The specific expressions of PDF h of inverse subordinator s(t) can be derived through the intimate links with subordinator t(s) [40]: .
Considering the formula that dΘ(x)/dx = δ(x), taking the partial derivatives of s or s 1 , s 2 in (8) and (9), together with Laplace transform (t → λ, t 1 → λ 1 , t 2 → λ 2 ), we obtain the PDF of h: and h(s 2 , λ 2 ; s 1 , λ 1 ) = ∂ 2 ∂s 1 ∂s 2 Note that the PDFs in (10) and (11) are both normalized, i.e., These PDFs of inverse subordinator play an important role in bridging the PDFs of the subordinated processes and original processes in a Langevin system.

Model
We consider the following set of Langevin equations: d dt where γ is the friction coefficient, ξ(s) is the Gaussian white noise satisfying ξ(s 1 )ξ(s 2 ) = 2D v δ(s 1 − s 2 ), and t(s) is the α-dependent subordinator (1 < α < 2) with Laplace exponent Φ(λ) ≃ µ 1 λ − µ α λ α introduced in section 2. This model will be investigated in three aspects in the following three subsections: firstly, we derive the Klein-Kramers equation based on Feynman-Kac equation; then we discuss the diffusive behavior of x(t) in two cases of γ = 0 and γ = 0.

Fractional Klein-Kramers equation
The fractional Klein-Kramers equation corresponding to (12) can be directly obtained from the forward Feynman-Kac equation in [30,49] (see also [50]), since x(t) = t 0 v(t ′ )dt ′ could be interpreted as a functional of v(t). The only difference with the one in [30] is that a new subordinator t(s) with 1 < α < 2 is considered here. Fortunately, the method in [30,49] can be applied to any subordinator with Laplace exponent Φ(λ), which only makes a difference in fractional substantial derivative operator proposed by [29]. In [29,51], t(s) is an α 0 -stable subordinator (0 < α 0 < 1) and the corresponding Laplace exponent is Φ 0 (λ) = λ α 0 . In this case, the fractional Klein-Kramers equation governing the joint PDF p(x, v, t) of position x and velocity v at time t is ∂ ∂t where L FP is the Fokker-Planck collision operator is the fractional substantial derivative operator defined as [29] as λ → 0 and ρ → 0. Taking the inverse Fourier-Laplace transform, we get the operator and obtain the fractional Klein-Kramers equation in the case of 1 < α < 2 Integrating over the position x, or making the Fourier transform (x → ρ) together with letting ρ = 0, the fractional equation governing the PDF of velocity is obtained, where D α−1 t is the fractional Riemann-Liouville derivative operator [52] with Laplace symbol λ α−1 , defined as The corresponding equation governing the PDF of positive x cannot be easily obtained by the similar procedure, since v is embedded into the fractional substantial derivative operatorD α−1 t , where the time t and position x are coupled with each other. Hence, it seems not easy to get the PDF p(x, t) and the moments of position x from the Fokker-Planck equation of x. Instead, we will calculate the moments straightly from the Langevin system (12).

Moments for the case γ = 0
In the case of γ = 0, the Langevin system (12) reduces to d dt which shows that v is a standard Brownian motion with respect to operation time s.
Denote v(s) as the velocity in operation time and v(t) := v(s(t)) in physical time. For convenience, we assume that the initial conditions are x 0 = v 0 = 0. So the odd-order moments of v and x are all zero. For the even-order moments of v and x, we can firstly calculate the correlation function v(s 1 )v(s 2 ) of velocity in (17) as The corresponding correlation function of v(t) in physical time t is given by [40,41] v This relationship is due to the fact that the process v(s) and subordinator t(s) are statistically independent. For convenience, we always make the calculations in Laplace space and obtain v( Substituting (11) and (18) into (20) gives as λ 1 , λ 2 → 0. Taking the inverse Laplace transform provides the second moment of v(t): Using the relation that dx(t)/dt = v(t) in (17), we get the correlation function of x(t): and its expression in Laplace space: after taking the inverse Laplace transform, which results in the second moment of x(t): Next, we calculate the fourth moment of x(t). Similarly, the four-point correlation function of v(t) should be presented firstly. In operation time s, where the integrand equals to [53] Similarly to (19), it seems that the four-point distribution h of inverse subordinator which might be too complicated or even unavailable. But following (24), it could be directly given as The formula (25) provides a shortcut and reduces the four-point distribution h to twopoint. But that (25) holds has the preconditional hypothesis s 1 < s 2 < s 3 < s 4 . We claim that (25) is valid on the condition that t 1 < t 2 < t 3 < t 4 and provide the detailed derivations in Appendix A. The techniques used will also work in other places. Two-point correlation function s(t 1 )s(t 2 ) can be directly calculated using (11). After some lengthly calculations in Laplace space, we get For the fourth moment of x(t), it can be written as Substituting the result (27) into above formula, we obtain The low-order moments of velocity v(t) and position x(t) have been obtained in (21), (23), (28) and (29), with their numerical simulations presented in figure 1. Note that these results only differ with the Langevin system without subordinator t(s) by a prefactor. This maybe understandable since the mean value of subordinator t(s) exists. In this sense, our subordinator t(s) just changes the time scale in a linear way by the parameter µ 1 for long times. But this is just a special case for γ = 0. In the following section, we consider γ = 0 and show the non-trivial moments.

Moments for the case γ = 0
In the case of γ = 0, the velocity v(s) in (12) is not a Brownian motion, but an Ornstein-Uhlenbeck process [53], which ensures a steady state of the diffusivity dynamics with respect to velocity for long times. The velocity process v(s) can be analytically expressed as which implies that the mean of v(s) is v(s) = v 0 e −γs , tending to zero for long times, and the correlation function Then the second moment of v in operation time s reads The second moment of v in physical time t can be obtained by using the relation [40] p where p(v, t) and p 0 (v, s) denote the PDFs of v(t) and v(s), respectively. Multiplying v 2 on both sides and integrating over v, together with Laplace transform and (10), we get for long times.
For the second moment of x, we resort to (11) and (20) and obtain v( and as λ 1 , λ 2 → 0. Inversing (33) and (34), the correlation function of v(t) and x(t) can be obtained (presented in Appendix B). Letting t 1 = t 2 = t there, we get for long times. The simulation results of the second moments of v(t) and x(t) with α = 1.3 and α = 1.7 are shown in figure 2, which are consistent to the theoretical results in solid lines for long times. It can be seen that the second moment of x(t) depends on α. This result is different from the case without the subordinator, in which [53] x 2 (t) ∝ t.
If we pay attention to the correlation function of x(t), the main difference from the case of γ = 0 in the previous subsection (22) is that γ = 0 here makes the asymptotic expression of the correlation function of v(t) (33) depend on α. More essential reasons about the difference the new subordinator brings in will be discussed in the next section. By adding a harmonic potential on v (i.e., γ = 0) in (12), the correlation function of v(t) is obtained in (32). This result can be extended to a system within an arbitrary confined potential U(v), where the steady state on v can be achieved. In this case, we denote the average of an observable O(v) on the Boltzmann distribution as is the normalizing function, and k B T the thermal energy. By imitating the method in [54], the correlation function of v(t) in confined potential U(v) can be presented in Laplace space as which recovers (32) when v 2 B = D v /γ and v B = 0 in model (12). When constructing single particle tracking experiments, the process x(t) is evaluated in terms of the time averaged MSD, defined via ∆ denoting the lag time. Typically, δx 2 (∆) is considered in the limit ∆ ≪ T to obtain good statistics. The correlation function of x(t) in the integrand depends on the correlation function of Alternatively, the time averaged MSD can also be obtained from the corresponding time averaged velocity correlation function [55] C and the Green-Kubo formula [56] Substituting the correlation function of v(t) in (39) into (40) and (41), we obtain the mean of the time averaged MSD for ∆ < < T . This result has been simulated in figure 3 with different α. It shows that the averaged quantity δx 2 (∆) experiences ∆ 2 in short time and ∆ 3−α in long time, consistent to the result in [55] about the Lévy walk. Define the ergodicity-breaking parameter as the ratio of time versus ensemble averaged MSD. Combining (35) and (42) shows which implies that the MSD is ultraweak ergodicity breaking, consistent to the results of Lévy walk in [55,57]. From the MSD of x(t) (35) and the time averaged MSD (42), we deem that the model (12) with γ = 0 describes the motion like Lévy walk. Especially for α-dependent subordinator with 1 < α < 2, it corresponds to the Lévy walk of sub-ballistic superdiffusion regime.  We observe that the regime of δx 2 (∆) changes from ∆ 2 to ∆ 3−α , which looks more obvious in (b) for larger α.

Relation with CTRWs and Lévy walk
In CTRWs, the motion of a particle is described by consecutive random waiting times between random jumps. The particle may undergo normal or anomalous diffusion, depending on whether the distributions are heavy-tailed or not. One special case is Lévy flight [58,59,60], where the waiting times have finite mean value but the jump lengths have infinite second moment. The possible disadvantages of Lévy flight are the diverging mean square displacement and the infinite velocity, which may be lack of physical meaning for a particle with finite mass. But Lévy walk avoids these drawbacks, where waiting time and jump length are coupled with each other. The standard Lévy walk says a particle moves ballistically for a random time and then randomly changes direction but keeps the same magnitude of velocity [34]. Therefore, in Lévy walk, much time penalizes a large jump and this balances the velocity to be finite. Fogedby [7] proposed the coupled Langevin equation (1) to describe the process in CTRWs, where ξ(s) and η(s) are independent with each other, characterizing the jump lengths and waiting times, respectively. Commonly, η(s) is taken to be one-sided α 0stable (0 < α 0 < 1) for describing the heavy-tailed waiting times distribution, and ξ(s) might be β-stable (0 < β < 2) for characterizing heavy-tailed jump lengths distribution in CTRWs. But for Lévy walk, its corresponding Langevin picture should be presented like (12), where the derivative of position x with respect to physical time t is velocity v and the subordinator t(s) characterizes the distribution of duration of each flight. The second equation in (12) gives the distribution of velocity v. One special case that η(s) is a one-sided α 0 -stable distribution (0 < α 0 < 1) and v(s) = γ −1 ξ d (s) has been pointed out in [44], where ξ d (s) is a dichotomous noise source, i.e., a random sequence of the values −1 and 1. It is just a one-to-one correspondence to the standard Lévy walk with the exponent of waiting time distribution less than 1. In general, the distribution of velocity v could be various, such as, Gaussian distribution, exponential distribution, and uniform distribution [36,61]. In more general cases, velocity v may be fluctuant due to a random force [62] and thus its distribution becomes time-dependent. All in all, velocity v can be described by a Langevin equation, i.e., the second equation of (12). The nonzero constant γ makes sure a steady state of velocity v could be reached for long times, analogously to the finite moments of v in Lévy walk. In a word, the overdamped Langevin equation with a subordinator (1) corresponds to CTRWs, while the weakly damped Langevin equation coupled with a subordinator (12) corresponds to Lévy walks. More generally, formula (37) in some sense implies that Lévy walk can also be modeled by arbitrary confined potential U(v) in velocity not only harmonic potential; the harmonic petential together with Gaussian white noise may be the simplest choice. For an asymmetric confined potential U(v), the biased Lévy walk together with the correlation function of v(t) can also be obtained.
For the coupled Langevin equation (12), one-sided α 0 -stable subordinator (0 < α 0 < 1) has been considered in [29,51], where the second moment of position x is and The result with γ = 0 is consistent to the standard Lévy walk in ballistic regime. Furthermore, we extend the subordinator to be α-dependent (1 < α < 2) and obtain the sub-ballistic superdiffusion regime (35) consistent to the corresponding Lévy walk. These two cases confirm the statement that the Langevin system (12) models the Lévy walk in long times. This subordinator of 1 < α < 2 has never been considered, but it is important to give rise to strong anomalous diffusion in this system, where the moments |x(t)| q exhibit different diffusion scales for different ranges of q [36,61]. We find an intriguing phenomenon that, compared to the second moments of position x with 0 < α 0 < 1 in (44) and (45), the diffusive behavior in (23) and (35) with 1 < α < 2 is enhanced for γ = 0 but suppressed for γ = 0. The Langevin system (12) with γ = 0 or γ = 0 are completely different models. Since α-dependent subordinator characterizes the heavy-tailed distribution of waiting times in CTRWs, it may yield longer waiting time for 0 < α < 1 than 1 < α < 2. For γ = 0, the subordinator suppresses the diffusion of velocity v and thus x due to the occasionally long waiting time, which implies the diffusion with 0 < α < 1 is suppressed more seriously. But for γ = 0, velocity v can reach a steady state for long time and the subordinator suppresses the rate of changing direction of particles and thus enhances the diffusion of displacement x, which results in a contrary result compared with γ = 0.
It is worth to note that the α-dependent subordinator does not always contribute to the strong anomalous diffusion phenomenon. Sometimes it makes a trivial result, like the case of γ = 0 in (17), where the moments of x(t) exhibit single diffusion scale. Actually, the position x(t) is a stochastic process with self-similarity, which can be briefly demonstrated. The α-dependent subordinator t(s) is 1/α self-similar [31], wherẽ α = α for 0 < α < 1 andα = 1 for 1 < α < 2. And then the inverse subordinator s(t) isα self-similar [63]. Therefore, the coupled velocity process where d = denotes identical distribution. Formula (46) implies v(t) isα/2 self-similar and thus x(t) isα/2 + 1 self-similar. In this way, This single diffusion scale indicates that there is no strong anomalous diffusion. But for γ = 0 or a more general nonlinear external force, which can be naturally added into the Langevin system, this subordinator might introduce a multiple diffusion scales and a different diffusion phenomenon.

Comparison with another Langevin system
Different from the Langevin system (12), another kind of commonly considered coupled Langevin system is d ds where position x and velocity v are both subordinated. If t(s) is the α-dependent subordinator (0 < α < 1), its corresponding fractional Klein-Kramers equation is proposed in [27]. Here we consider the case of 1 < α < 2, and the fractional Klein-Kramers equation will be different from (14). Denote the joint PDF of position x and velocity v in operation time as p 0 (x, v, s) and the one in physical time p(x, v, t). Then p 0 (x, v, s) solves the Klein-Kramers equation [2] ∂ ∂s Using the relation This is the fractional Klein-Kramers equation governing the joint PDF of positionvelocity of the Langevin system (48). Note that in this case, the Newton relation does not hold between x(t) and v(t) and Galilean invariance is violated [6,64].
Integrating over the position x on (49), we get the same equation governing the PDF of velocity v(t) as (15). This is reasonable since the only difference between the Langevin system (12) and (48) is the position x(t). But here, we can also derive the equation governing the PDF of position x(t) by integrating (49) over dv and vdv, and combining the two resulted equations. With v 2 (t) ≃ D v /γ in (31) for the case of γ = 0, this procedure yields the fractional diffusion equation of p(x, t): which becomes, in the long time or high-friction limit, It can be seen that in the long time limit, the Langevin system (48) undergoes normal diffusion, with the odd-order moments vanishing and even-order moments as Another way to derive the moments of x(t) (52) is based on the Gaussian distribution of the original process of x(s) in operation time. For a Gaussian process, its PDF can be completely determined from the knowledge of its variance and mean. Based on the correlation function of v in operation time (30), we calculate the second moment of x(s) in operation time for model (48): 2D v γ 2 s. Since the mean of v(s) is zero, the motion is unbiased and the odd-order moments of x(s) are zero; the even-order moments are Then using the relation and the asymptotic expression Φ(λ) ≃ µ 1 λ, as λ → 0, one can also get the result (52). In the long time, the Langevin system (48) coupled with α-dependent subordinator (1 < α < 2) still exhibits normal diffusion (52) as in the operation time (53), although this subordinator might change the PDF of x(t) and v(t) in the Langevin system. At first glance, the α-dependent subordinator (1 < α < 2) has finite mean, and might make no difference with the exponential distribution or simply without any subordinator. This recognition is correct just in some special cases, e.g., the coupled Langevin system (48). But for most of complex system or various statistical quantities, this subordinator may still bring in some new interesting phenomena or diffusion behavior, which essentially depend on whether the observed statistical quantities are related to the multiple-point distribution of the inverse subordinator. For the simple cases, some quantities of the subordinated processes might only depend on the singlepoint distribution of inverse subordinatorĥ(s, λ) in (10), where Φ(λ) ≃ µ 1 λ, as λ → 0, for long times, just like the procedure (54). But if it depends on the two-point distribution of inverse subordinatorĥ(s 2 , λ 2 ; s 1 , λ 1 ) in (11), where , the result will be different.
In this sense, it is not hard to understand why Lévy walk exhibits a special subballistic superdiffusion regime when the exponent of waiting times is 1 < α < 2 while a trivial phenomenon is observed for CTRWs because of the boundedness of the first moment of the waiting time distribution. The second moment x 2 (t) in Lévy walk depends on the correlation function of velocity v(t 1 )v(t 2 ) and thus the two-point distributionĥ(s 2 , λ 2 ; s 1 , λ 1 ). But x 2 (t) in the case of CTRWs only depends on x 2 (s) and thus depends on single-time distributionĥ(s, λ). So if we consider the overdamped Langevin equation (1) with such a subordinator, the diffusion behaviours for long times will be the same as the ones of the original process.

Numerical simulations of subordinator
Below, we show how to numerically approximate the sample paths of the process x(t) of (12). In the first step, we numerically approximate the α-dependent subordinator t(s) with 1 < α < 2 on the lattice {τ k = k∆τ : k = 1, · · · , N} where ∆τ = T /N. For making some preparations, let us give a brief introduction of the idea of two time scales in [65].
Suppose that X 1 , X 2 , X 3 , · · · are the sequence of independent identically distributed positive random variables representing the waiting times between consecutive jumps of the walker, with the distribution [34] which is consistent to the Lévy measure ν(dy) defined in (4). The Laplace transform of the PDF φ(t) iŝ as λ → 0. Note that X i has a positive mean µ 1 . Consider the total time with the scale factor c, and [ct] denotes an integer number satisfying [ct] ≤ ct < [ct] + 1.
Note that the first sum grows like c 1/α while the second grows like c as c → ∞. Hence, we cannot get a convergence by normalizing only at one scale. So the Lévy (and Gaussian) central limit theorem [66] is not valid here. Instead, we use the technique in [65] of normalizing T [ct] at two scales, and get the centered and normalized sum Note that T c (t) cannot represent the time of [ct]-th jump for large c, since it is not non-decreasing. This can be verified from the increment of T c (t) that for 1 < α < 2, Taking c → ∞ in T c (t), we obtain the Lévy process T (t), but it is not nondecreasing. Therefore, we consider the non-decreasing supremum processT (t) defined as [65]T Since the first-passage times (i.e., inverse subordinator defined in (7)) of the process T (t) and its supremum processT (t) are the same, we will generate the inverse subordinator s(t) based on T (t) in numerical simulations. This is appropriate and can be verified from the double Laplace transform of the PDF ofT (t) given in [65,67] that as λ → 0 and q is a holomorphic function. Taking the inverse Laplace transform u → t of (58), the characteristic function ofT (t) is obtained as e −t(µ 1 λ−µαλ α ) , consistent to the Laplace exponent Φ(λ) in (5).
Following the discussions above, we can generate the random variable X i drawn from the distribution (55) by where U is uniformly distributed between 0 and 1. Then we get the mean µ 1 of X i : and thus obtain the centred and normalized sum as (56) which is the approximation value of the subordinator t(s) at lattice τ k , k = 1, · · · , N. Based on the subordinator of (59), we use the methods in [68] to generate the inverse subordinator process s(t) in (7) and the subordinated process v(t) = v(s(t)) in (12). Supposing that we have got the velocity v(t i ) on another set of lattices {t i = i∆t : i = 0, 1, · · · , M}, the position x(t) can be obtained directly by For the inverse subordinator s(t), we simulate its first two moments. From (10) and Φ(λ) ≃ µ 1 λ for small λ, we get h(s, λ) ≃ µ 1 e −µ 1 λs , and thus Taking inverse Laplace transform (λ → t) gives Figure 4 shows the numerical simulations of the first and second moments of s(t). It can be seen that the simulation results (circle markers and square markers) are consistent to the theoretical results (63) (solid lines) for long times, which also verifies the long time asymptotic approximation in (58) for subordinator.

Summary and conclusions
Lévy walk is an important model for describing random walk with finite velocity, which exhibits anomalous superdiffusion phenomenon. For the standard Lévy walk, it can be divided into three categories, depending on the value of the power-law exponent β of the waiting times distribution: ballistic diffusion for 0 < β < 1, sub-ballistic superdiffusion for 1 < β < 2, and normal diffusion for β > 2. Based on the feature of finite velocity, we claim that the weakly damped Langevin system in a confined potential together with a subordinator on velocity v can model the dynamics (almost the same as the  [29,51], and presented its corresponding Langevin picture in [26], where a weakly damped Langevin system is coupled with one-sided α-stable subordinator (0 < α < 1). In [29], the second moment x 2 (t) ∝ t 2 was obtained, which is consistent to the ballistic regime of Lévy walk. Another way to characterise Lévy walk from overdamped Langevin equation is to assume that jump sizes are some functions of waiting times in [69].
In this paper, we define a new α-dependent subordinator (1 < α < 2) and provide its simulation method when applied to Langevin systems. The weakly damped coupled Langevin equation with this subordinator is build. The lower-order moments of velocity v(t) and position x(t) in this Langevin system are calculated. Especially for γ = 0, the diffusion behaviours of the Langevin system are the same as the ones of the Lévy walk in sub-ballistic superdiffusion regime for long times, where strong anomalous diffusion can be observed. There is a relatively intuitive interpretation for this regime in Lévy walk [34]. Compared with the case of 0 < α < 1, there are less particles still in their very first flights to form the ballistic fronts when 1 < α < 2. So they are slower than the ones of 0 < α < 1, but still faster than normal diffusion. Here we present the interpretation from the perspective of Langevin system. For general simple cases, this subordinator might be trivial, working like a linear transform for long times. This is because that the exponent of the single-point distribution of inverse subordinator in Laplace space reduces to be linear with λ. But for the model (12) with γ = 0, the second moment of position x(t) depends on the two-point distribution of inverse subordinator, where λ α plays an important role, and eventually contributes to the sub-ballistic superdiffusion regime.
The overlooked α-dependent subordinator with 1 < α < 2 helps to model the motion of Lévy walk in sub-ballistic superdiffusion regime. Keeping this essential/potential mechanism in mind, it will be helpful to characterize more complex stochastic processes with this subordinator, e.g., turbulent in fluids, complex liquids and various biological system. Besides, more complex Langevin system with this subordinator can be considered, such as, the system with a nonlinear external force field, and even the functional distribution of the particle trajectory in the weakly damped or overdamped system. Appendix A. Derivation of (25) For simplicity, we denote h 4 (s, t) and g 4 (t, s) as the four-point distribution of inverse subordinator s(t) and subordinator t(s), respectively. So the relation of correlation function between operation time s and physical time t is
Appendix B. Correlation functions of v(t) and x(t) from (33) and (34) Here we derive the correlation functions of v(t) and x(t) by inversing (33) and (34), respectively. Since λ 1 , λ 2 → 0, we only give the results for long times and further assume t 1 < t 2 without loss of generality. Taking the inverse Laplace transform (λ 1 → t 1 , λ 2 → t 2 ) of the three terms in (33), respectively, yields Note that t 1 < t 2 has been used in the inverse Laplace transform of these three terms.