Novel anomalous diffusion phenomena of underdamped Langevin equation with random parameters

The diffusion behavior of particles moving in complex heterogeneous environment is a very topical issue. We characterize particle's trajectory via an underdamped Langevin system driven by a Gaussian white noise with a time dependent diffusivity of velocity, together with a random relaxation timescale $\tau$ to parameterize the effect of complex medium. We mainly concern how the random parameter $\tau$ influences the diffusion behavior and ergodic property of this Langevin system. Besides, the comparison between the fixed and random initial velocity $v_0$ is conducted to show the effect of different initial ensembles. The heavy-tailed distribution of $\tau$ with finite mean is found to suppress the decay rate of the velocity correlation function and promote the diffusion behavior, playing a competition role to the time dependent diffusivity. More interestingly, a random $v_0$ with a specific distribution depending on random $\tau$ also enhances the diffusion. Both the random parameters $\tau$ and $v_0$ influence the dynamics of the Langevin system in an non-obvious way, which cannot be ignored even they has finite moments.


Introduction
Anomalous diffusion phenomenon is very common in both living and abiotic systems [1,2,3,4,5]. It is characterized by the significant deviation from Brownian motion with the nonlinear time dependence of the ensemble-averaged mean-squared displacement (EAMSD) [6,7]: Different range of the anomalous exponent µ determines a more detailed division of the anomalous diffusion phenomenon: µ < 1 representing the subdiffusion and µ > 1 the superdiffusion. The subdiffusion phenomena of passive particles abound in the cytoplasm of living biological cells [8,9,10] and in artificially crowded fluids [11,12,13], as well as in quasi-two-dimensional systems such as lipid bilayer membranes [14,15,16,17]. The superdiffusion is typically associated with active matter and can also be observed in living cells [18,19,20]. Modern microscopy techniques have made it possible to track the random motions of tracer particles in a variety of living and abiotic systems. The universal existence of the anomalous diffusion phenomena in the natural world inspires the scientists to devise the possible physical models and explore the underlying mechanisms of them. The most typical model of describing particle's motion is continuous-time random walk (CTRW) [6,21,22,23], which is characterized by waiting time and jump length. By assuming the specific distributions of waiting time and jump length, together with the possible correlation between them, the CTRW can yield many kinds of diffusion processes, which present anomalous diffusion phenomena with different anomalous exponents [24,25,26,27,28]. However, more and more new classes of diffusive dynamics have been observed, especially in the field of soft matter, biological, and other complex systems. One typical example is the Brownian non-Gaussian phenomenon, where the EAMSD is normal but the probability density function is non-Gaussian [29,30,31,32,33]. The effective idea of interpreting this novel phenomenon is to introduce a random diffusivity which is analogous to the concept of superstatistics [34,35,36]. It says that each one of particles moving in a complex heterogeneous environment has its own diffusivity [37,38,39,40]. In general, the idea of superstatistics has been applied to many dynamics: turbulent dispersion [34], renewal critical events in intermittent systems [41,42], and different effective statistical mechanics with fluctuating intensive quantities [34,35].
This paper considers an underdamped Langevin system driven by a Gaussian white noise with a time dependent diffusivity of velocity. The effect of the complex medium on particles is parameterized via a random relaxation timescale τ . Each particle has its own fixed relaxation timescale τ as it moves in the complex medium, but the τ might be discrepant between different particles. For the classical Brownian motion with a constant relaxation timescale τ , the EAMSD evolves from ballistic diffusion to normal diffusion, with the crossover point at time τ . One can imagine that a random relaxation timescale τ might yield novel dynamics especially when τ has a heavy-tailed distribution. On the other hand, to maintain the good statistics of the Langevin system with constant τ , we suppose that the first moment τ is finite. Considering the two conditions, the relaxation timescale τ is assumed to obey the modified Lévy distribution as (4) shows.
We mainly concern the ergodic property of the underdamped Langevin system, which is embodied by the time-averaged mean-squared displacement (TAMSD), defined as [43,44,45] Here, ∆ is the lag time, T is the total measurement time. To obtain good statistical properties, the lag time is assumed to be much shorter than the measurement time, i.e., ∆ ≪ T . Compared with EAMSD, the TAMSD places emphasis on particle-toparticle diffusion properties and contains more information of particle's motion in a single trajectory. Based on the advance of single-particle tracking techniques, scientists often evaluate the recorded time series in terms of TAMSD to study the diffusion behavior of particles in living cells [8,46,47]. A stochastic process is called ergodic if its TAMSD and EAMSD are equivalent, i.e., δ 2 (∆) = x 2 (∆) as the measurement time T → ∞, such as Brownian motion and (tempered) fractional Brownian motion [48,49,50]. But for anomalous diffusion processes, especially for the molecules diffusing in living cells, the time average often becomes a random variable and irreproducible for the individual trajectories, implying the ergodicity breaking, such as the continuous-time random walk [9,51,52], Lévy walk [53,54,55], heterogeneous diffusion process [56,57,58,59], and random diffusivity process [60,61,62], etc.
When investigating the ergodic property of a diffusion process, the effect of the initial condition usually cannot be ignored. The time average is in some sense an equilibrium measure, while the ensemble average is not. The different diffusion behavior resulting from the discrepant initial ensemble has been observed in the context of single file diffusion [63], diffusing diffusivity model [64], Lévy walk [65,66,67], biased random walk [45,68]. For the underdamped Langevin equation with random relaxation timescale τ , we also study the effects of different initial velocity v 0 , which is also the main part of this paper.
The key of evaluating the EAMSD and TAMSD is to obtain the velocity correlation function first, which has been deeply influenced by the random relaxation timescale τ . Hence, the diffusion behavior and the ergodic property are closely related to the parameter τ . Fixing the initial velocity v 0 , we find that the random τ enhances the diffusion behavior while the time-dependent diffusivity slows down the diffusion behavior. They play the competition roles and present different diffusion phenomena. On the other hand, when the initial velocity becomes random, especially depends on the random parameter τ , the random initial velocity also enhances the diffusion behavior as the random relaxation timescale τ .
The structure of this paper is as follows. In section 2, we establish the underdamped Langevin system driven by the Gaussian white noise with a time dependent diffusivity. Then we evaluate the EAMSD and TAMSD, and make comparison between the cases with fixed and random relaxation timescale τ in section 3, and between the cases with fixed and random initial velocity v 0 in section 4. Some discussions and summaries are provided in section 5. For readability, some mathematical details are present in Appendix.

Model
The underdamped Langevin equation driven by a Gaussian white noise with a timedependent diffusivity is written aṡ where τ is the relaxation timescale, νt β−1 is the time-dependent diffusivity with 0 < β < 1, ξ(t) is Gaussian white noise with zero mean value and correlation . This model can be regarded as the underdamped form of scaled Brownian motion, which describes the motion of a particle in a viscous medium with time-dependent temperature. Another kind of underdamped scaled Brownian motion described by Langevin equation with internal noise is studied in [69], where the time dependent diffusivity is related to the damping coefficient according to the (time local) fluctuation-dissipation theorem [70]. The novelty of the Langevin equation (3) is that the relaxation timescale τ is not a constant any more, but a random variable, due to the complexity of the media. Each particle has its own relaxation timescale τ , which might discrepant for different particles. That is to say, the parameter τ is random for ensemble of particles, but fixed for each particle or trajectory. Our aim is to find how the random parameter τ influence the diffusion behavior and ergodicity of the Langevin system (3) by evaluating the EAMSD and TAMSD of particles. Besides, the initial velocity v 0 of particles can also be random, which might lead to a different phenomenon of the ergodicity from the cases with a fixed v 0 .
Throughout this paper, we consider the relaxation timescale τ obeying the distribution [71] where L α (·) is the one-sided fully skewed α-stable Lévy distribution with 0 < α < 1 [22,72] and τ is the mean value of the random parameter τ . The advantage of choosing such a kind of distribution is utilizing its heavy-tailed property g(τ ) ∝ τ −2−α for large τ . This property can yield a very large value of variable τ , which may change the diffusion behavior of the Langevin equation (3). The prefactor 1/τ in front of the Lévy distribution is aimed to guarantee a finite mean τ of the relaxation timescale. Otherwise, the variance of velocity would be infinite [71]. When simulating the random variable τ , the parameter τ controls the width of the Lévy distribution and can be chosen arbitrary. For completeness on the description of the distribution in (4), we present its normalization and the mean value τ , together with the algorithm of generating τ in Appendix A.
The key of obtaining the EAMSD and TAMSD is the velocity correlation function v(t 1 )v(t 2 ) . To obtain it, we firstly solve (3) in Laplace domain, and achieve the analytic expression of velocity process which is valid for both fixed and random v 0 and τ . When v 0 and τ are random variables, instead of v(t 1 )v(t 2 ) , we express the velocity correlation function in the condition of parameters τ and v 0 , i.e., the conditional velocity correlation function v(t 1 )v(t 2 )|v 0 , τ = e − t 1 +t 2 for t 1 ≤ t 2 , where the ensemble average is only performed over the Gaussian white noise ξ(t), and 1 F 1 (a, b, z) is the confluent hypergeometric function.

Fixed and random relaxation timescale τ
In this section, we fix constant initial velocity v 0 and compare the two cases with fixed and random parameters τ by investigating the diffusion behavior and ergodicity of the Langevin equation (3). Let us firstly consider the constant relaxation timescale τ and initial velocity v 0 . The velocity correlation function, as well as the variance of the velocity process, can be evaluated directly based on (6) for large t 1 , i.e., where the asymptotic formula of hypergeometric for large z has been used. Fixing t 1 and putting t 2 to the infinity in the first equation of (7), the velocity correlation function decays exponentially as the Ornstein-Uhlenbeck process, although the diffusivity becomes time dependent in our model (3). For 0 < β < 1, however, the velocity variance v 2 (t) tends to zero at a power law rate in the second equation of (7), while the one of the Ornstein-Uhlenbeck process converges to a constant. Further based on the velocity correlation function, the EAMSD of the stochastic process with fixed parameter τ is which is independent of the initial velocity v 0 in the long time asymptotics. The stochastic process described by the Langevin equation (3) with fixed parameter τ exhibits subdiffusion behavior since 0 < β < 1. The simulation results of the EAMSD with different β are shown in figure 1, which tend to the theoretical lines for large t.
When evaluating the ensemble-averaged TAMSD, we should additionally calculate the aging EAMSD (x(t + ∆) − x(t)) 2 with the initial observation time t which is much longer than the subsequent evolutionary time ∆. More precisely, for fixed parameter τ , the aging EAMSD can be obtained by integrating the velocity correlation function (6) as where the detailed derivations are present in Appendix B. The result in (9) shows that the aging EAMSD exhibits the transition from ballistic diffusion at short time to normal diffusion at long time. The dependence of the initial observation time t implies the aging behavior of the process x(t). Based on the aging EAMSD, the ensemble-averaged TAMSD of the stochastic process for ∆ ≪ T can be obtained by use of the definition of TAMSD in (2) as which exhibits normal diffusion behavior with respect to the lag time ∆. It shows the pronounced difference from the EAMSD in (8), and thus implies the nonergodicity behavior of the Langevin equation (3). Figure 2 shows the simulation results of the ensemble-averaged TAMSDs with different β, which coincide with the theoretical results well. Now we turn attention to the case with random parameter τ obeying the distribution (4). The velocity correlation function averaged over random parameter τ can be expressed as Using the residue theorem (with detailed derivations in Appendix C), we obtain the asymptotic behavior of velocity correlation function for large t 1 , t 2 and t 1 < t 2 : where The Gaussian hypergeometric function is defined as where (q) n is the (rising) Pochhammer symbol, and the asymptotic expression for small z is: 2 F 1 (a, b, c, z) ≃ 1 + ab c z. Therefore, for fixed time t 1 and t 2 → ∞, the velocity correlation function in (12) scales as By observing the velocity correlation function (14), one finds that the correlation of velocity process at two different times decays at a power law rate, presenting the significant difference from the exponential decaying rate in the case with fixed parameter τ and time-dependent diffusivity in (7). The decay rate is controlled by the Lévy exponent α in the distribution g(τ ) in (4), while the coefficient is effected by the timedependent diffusivity with power β. Based on the velocity correlation function (12), the EAMSD can be obtained through double integrals as (8) shows: where The detailed derivations are present in Appendix D. Note that the EAMSD (15) for random parameter τ cannot be got by performing the emsemble average over parameter τ on the EAMSD (8) directly, since the latter is obtained with a series of asymptotic calculations conditioned on fixed τ . More analyses on the ineffectivity of the direct ensemble average over the random variable τ will be presented in the next section.
Considering the relationship between the power law exponents in two terms of (15), which is 1 − α < 1 − α + β, the second term plays the dominating role at long time and the effect of the initial velocity v 0 can be omitted, in other words, the EAMSD of the Langevin equation (3) behaves as The coincidence between the theoretical and simulated results of the EAMSD is shown in figure 3 with different α and β. Since the diffusion exponent satisfies 0 < 1−α+β < 2, the stochastic process described by Langevin equation (3) with both random relaxation timescale τ and time-dependent diffusivity t β−1 can show various types of diffusion behavior, including normal diffusion (α = β), subdiffusion (α > β) and superdiffusion (α < β).  The Langevin equation (3) recovers the Brownian motion and the EAMSD (16) becomes normal when α = β = 1 . Compared with this normal case, although the parameters α and β are both in the range of (0, 1), they provide the opposite effects on the diffusion behavior, which can be inferred from the opposite signs in front of α and β in the exponent of t in (16). More precisely, the random relaxation timescale τ with parameter α enhances the diffusion while the time-dependent diffusivity with parameter β suppresses the diffusion. The promoting to the diffusion of the random parameter τ can also be found by comparing the EAMSDs between (16) and (8) for the cases with random and deterministic parameter τ , respectively. Since 0 < α, β < 1, the one, which is smaller or more deviates from 1, plays the dominating role in the diffusion behavior.
Based on the velocity correlation function v(t 1 )v(t 2 ) (12), similar to (9), the aging EAMSD is, for ∆ ≪ t, which depends on the time t and indicates the aging behavior of the process x(t). The interesting thing is that the parameter α of random variable τ only acts on the lag time ∆ and does not changing the aging behavior. The aging EAMSD (17) shows the superdiffusion behavior, which is faster than the case with fixed parameter τ showing normal diffusion (9). Substituting the aging EAMSD into the definition of TAMSD and calculating the integral, we arrive at the ensemble-averaged TAMSD The corresponding simulation results are presented in figure 4. Comparing the ensembleaveraged TAMSDs in (10) and (18), one can see that the randomness of parameter τ accelerates the diffusion behavior in the sense of TAMSD (from order 1 to order 2 − α). The inequality of the EAMSD (16) and TAMSD (18) means the ergodicity breaking of the Langevin equation (3) with respect to the MSD.

Fixed and random initial velocity v 0
The results in section 2 are obtained in the condition of a fixed initial velocity v 0 , where the two terms in the EAMSD (15) correspond to those in the velocity correlation function v(t 1 )v(t 2 ) (12), respectively, one containing v 2 0 and another one not. The first terms in (12) and (15) can be neglected compared with the corresponding second terms when v 0 is a constant. Actually, for the case with random initial velocity v 0 , as long as it is independent of the relaxation timescale τ and has a finite second moment, the first terms in (12) and (15) can be omitted compared with the second terms after performing the ensemble average with respect to v 0 .
On the other hand, if the initial velocity v 0 is not only random, but also depends on relaxation timescale τ , then the ensemble average with respect to τ on (11) will also work on v 2 0 and the result will be different. Let us consider a special case v 2 0 = ντ , which implies that the two random variable v 0 and τ are closely related and v 0 has a finite second moment. For the particle's motion described by Langevin equation with random v 0 and τ satisfying v 2 0 = ντ , each particle has its own fixed initial velocity v 0 and relaxation timescale τ throughout its whole trajectory, and the pair of parameters (v 0 , τ ) satisfy the quantitative relation v 2 0 = ντ . But different particles have discrepant parameter pair (v 0 , τ ).
Before performing the calculations of EAMSD and TAMSD for the case of v 2 0 = ντ , let us see the reason of choosing such a v 0 . Taking β = 1 in (6) leads to the conditional velocity correlation function of the Langevin system with time-dependent diffusivity: where the average is only made over the Gaussian white noise ξ(t). For fixed v 0 and τ , the second term on the right-hand side dominates the velocity correlation function at long time, that is, the velocity correlation function is asymptotic stationary. When τ is random, however, the distribution of τ in (4) implies that τ can take an arbitrary large value, and the first term cannot be ignored even for large t 1 and t 2 . But when v 2 0 = ντ , the first term vanishes and the conditional velocity correlation function is stationary from the beginning of the evolution whether τ is fixed or random. Although we consider 0 < β < 1 in this paper, the first term containing v 2 0 in velocity correlation function (11) is independent of parameter β. That is why we still choose the initial condition satisfying v 2 0 = ντ . In the case of v 2 0 = ντ for any 0 < β < 1, when performing the average over random parameter τ in (11), the first term containing v 2 0 in (12) should be modified as which is obtained by using residue theorem with the similar procedures in Appendix C. When evaluating the EAMSD by performing the double integrals over velocity correlation function v(t 1 )v(t 2 ) , the term in (20) contributes to Similarly, the term in (20) has another contribution to the ensemble-averaged TAMSD. Here, we use the subscript "0" to denote the part coming from random initial velocity v 0 . Comparing the EAMSDs between (16) and (21), we find that the latter contributed by the random initial velocity v 0 plays the dominating role, i.e., the total EAMSD is However, the magnitude relation between ensemble-averaged TAMSDs in (18) and (22) cannot be easily determined. To verify the correctness of the partial ensemble-averaged TAMSD (22) contributed by the random initial velocity, we make the simulations on the quantity where the first term on the right-hand side is the real ensemble-averaged TAMSD and the second term is the part coming from Gaussian white noise obtained in previous section as (18) shows. The corresponding simulations for EAMSD and ensemble-averaged TAMSD are present in figure 5, which are consistent to the theoretical results (23) and (22), respectively.  The ensemble-averaged TAMSD becomes simple for the case with time-independent diffusivity (β = 1). In this case, both the EAMSDs in (16) and (21) behaves as t 2−α , and the random initial velocity v 0 does not change the diffusion behavior. As to the ensemble-averaged TAMSD, the second term ∆ 2−α on the right-hand side of (24) plays the leading role compared with the term ∆ 2 T α in δ 2 (∆) 0 . The difference between the EAMSD (23) and the ensemble-averaged TAMSD (22) or (18) implies the nonergodicity of the model (3) with random initial velocity v 0 and random relaxation timescale τ . However, the special case β = 1 may lead to a misunderstanding, where the EAMSD and the ensemble-averaged TAMSD both behave as It seems that the Langevin equation (3) is ergodic due to the equivalence between the EAMSD and the ensemble-averaged TAMSD as (25) shows. However, the TAMSD is usually a random variable for anomalous diffusion processes. To detect the ergodicity accurately, we should evaluate the TAMSD to judge whether the TAMSD converges to its ensemble average as the measurement time T → ∞. For the Langevin equation (3) with random parameter τ , each particle (or trajectory) has a fixed τ , namely, τ i for the i-th particle. Therefore, each particle undergoes the Brownian motion with a fixed relaxation timescale τ i . But different particles own different parameter τ i , the specific value of which is yielded from the distribution g(τ ) in (4). Further due to the ergodicity of Brownian motion, the TAMSD of each particle converges to the ensemble-averaged TAMSD with fixed relaxation timescale τ i , i.e., which is the classical result of Brownian motion [73]. It can also be obtained from the aging EAMSD (9) by taking β = 1 and performing the integral with respect to t. We make the sample of two trajectories with τ 1 = 0.22 and τ 2 = 0.57, and present the corresponding TAMSDs in figure 6, which agree with the theoretical lines very well. For both short time and long time, the TAMSD of two trajectories are parallel, since the difference between them is just the relaxation timescale τ i . The inconsistence of two trajectories implies the nonergodicity of the Langevin equation (3) when β = 1. When we observe the ensemble-averaged TAMSD for large ∆, it seems strange that the result ∆ 2−α in (25) cannot be obtained by directly performing average over the result of asymptotic TAMSD (∆ 2 or ∆) in (26). Actually, due to the randomness of parameter τ , each trajectory owns a different but fixed τ i . Therefore, given ∆, the large amounts of trajectories can be divided into two categories, one with τ i larger than ∆ and another one smaller than ∆. Thus, the ensemble average should not be performed solely on any one of the terms in (26), but where the two integrands come from the different asymptotics in (26). The two terms in (27) can both be verified to scale as ∆ 2−α as (25) shows. The detailed calculations are presented in Appendix E. The analyses around (27) are also valid for the case β < 1, where the TAMSD still has the good property of self-averaging. According to (9), the corresponding TAMSD for each trajectory is which is different from the result (26) by a constant. The self-averaging property of the TAMSD and its equivalence to the theoretical result (28) has been verified through the simulations as figure 7 shows. The aim of making different assumptions about the initial velocity v 0 is to detect how the initial ensemble influences the diffusion behavior. For the underdamped Langevin equation describing the classical Brownian motion with fixed τ , a fixed v 0 is named non-equilibrated initial ensemble, which means all particles are introduced to the system at t = 0. On the other hand, if the system reaches a certain equilibrium before we start measuring it, namely, the initial velocity v 0 has reached at the equilibrium distribution at t = 0, then it is named equilibrated initial ensemble [74]. Different initial ensemble usually leads to different diffusion behavior, which has been discussed a lot [45,63,64,65,66,67,68]. It is noteworthy that the way of initial ensemble affecting the diffusion behavior is quite special for the underdamped Langevin equation (3) with random relaxation timescale τ .

Summary
The rich dynamics of particle's motion in biological systems have been attracting the attention of many researchers in the field of statistical physics. The complexity of the cell environment brings many difficulties to the research of building models to describe the particle's motion and further theoretical analyses. In this paper, the particle's motion is characterized by an underdamped Langevin equation with a time-dependent diffusivity of velocity, in which the random relaxation timescale τ describes the complexity of the medium, and the random initial velocity v 0 represents a special initial ensemble. The aim of this paper is to detect how the two parameters, relaxation timescale τ and initial velocity v 0 , influence the diffusion behavior and ergodic property of the Langevin system.
For the relaxation timescale τ , we find its randomness slows down the decay rate of the velocity correlation function, and enhances the diffusion behavior of the Langevin system (3), which is opposite to the effect of the time-dependent diffusivity t β−1 , 0 < β < 1. Therefore, the exponents α and β play the competition roles and the stochastic process can present different diffusion behavior (from subdiffusion to superdiffusion) with different values of α and β. The smaller one between α and β plays the dominating role in the diffusion process. As to the ergodic property, it is the timedependent diffusivity t β−1 , rather than the random relaxation timescale τ , who leads to the aging behavior. Therefore, the Langevin system is nonergodic whether τ is fixed or random.
For the initial velocity v 0 , its randomness brings a very novel phenomenon. We mainly study the special case v 2 0 = ντ , which implies that the two random variable v 0 and τ are closely related. The random v 0 enhances the diffusion behavior in (21) (from t 1−α+β to t 2−α ) with respect to EAMSD. It also contributes to a ballistic term ∆ 2 T −α (22) to the ensemble-averaged TAMSD which cannot be omitted even for long time. Further, we investigate the special case β = 1, i.e., the diffusivity of velocity becomes a constant. Surprisingly, the random initial velocity v 0 does not change the diffusion any more when β = 1.
The random parameters in physical models are sometimes underestimated both from the aspects of the calculation procedure and the dynamic behavior. For the former aspect, we point that the quantity cannot be evaluated by directly performing ensemble average over the random parameters on the original quantity, and analyse the reasons around (27). For the latter aspect, although the relaxation timescale τ of the Langevin system (3) has finite first moment and the initial velocity v 0 has finite second moment, both of them change the diffusion behavior of particles and bring us the novel phenomena. The detailed derivations and meticulous analyses in this paper enormously enrich the theoretical part of physical models with random parameters.
Thus, the mean value of τ is where we have used the variable substitution in the second equality and the normalization of Lévy distribution in the last equality. Then we turn to the normalization of the random variable τ , where some technical procedures will be performed. The Laplace transform of g(τ ) iŝ where we have used the formula e −λτ /τ = The key of simulating the trajectories of Langevin equation (3) is to generate the random variable τ obeying distribution g(τ ). The first step of generating τ is to obtain the function L α (x) numerically (i.e., the values L α (x i ) at every point x i ) by performing the inverse Laplace transform on e −λ α through MATLAB. Then we obtain the function g(τ i ) numerically through linear transform on L α (x i ) with x i = Cτ i . Further, we sum g(τ i ) to obtain the numerical cumulative distribution function F (τ i ), together with its inverse function F −1 (z), z ∈ [0, 1]. Finally, by use of the well-known inverse transform sampling method, we can generate the random variable τ obeying the distribution g(τ ) in the way that where U obeys the uniform distribution in [0, 1].
Appendix B. Derivation of the aging EAMSD (9) The integrand of (9) can be rewritten as where the first two terms are known in (8), and the last term can be evaluated as We first calculate the integral with respect to t 2 on interval (t, t + ∆), and obtain It holds that for large t, since the Laplace transform (t → λ) of the left-hand side in (B.4) is Γ(β)λ −β λ+1/τ , which asymptotically equals τ Γ(β)λ −β for small λ.
for large t and ∆ ≪ t.
Appendix C. Derivation of the velocity correlation function (12) In the expression of the conditional correlation function v(t 1 )v(t 2 )|v 0 , τ (6), one finds that the terms containing parameter τ have the exponential form e − f τ . By use of the notation C = the ensemble average on the exponential term e − f τ over parameter τ is Denote the above integration as R(f ). It can be solved by using the residue theorem with poles s/α + 1 = −n or 1 − s = −n, where n = 0, 1, 2, ....

Appendix D. Derivation of the EAMSD (15)
We only present the derivation of the second term of the EAMSD (15), which is denoted as x 2 (t) 2 here. For long times, using the integral representation of the second term in velocity correlation function (C.4), i.e., we have Appendix E. The evaluations of two terms in (27) Denote the two terms in (27) as I 1 and I 2 , respectively. For I 1 , we have Substituting the expression (4) of g(τ ) into (E.1), and performing the variable substitution yield the large ∆ asymptotics of which can be evaluated by using Laplace transform. The Laplace transform (C∆ → λ) of the integral is for small λ, and its inverse Laplace transform is α(C∆) 1−α Γ(2−α) . Therefore, Similarly, I 2 can be calculated as where we evaluate the integral in second equality by using Laplace transform (C∆ → λ), i.e., 1 λ − 1 λ e −λ α ≃ λ α−1 , the inverse Laplace transform of which is (C∆) −α Γ(1−α) . Therefore, based on the approximation in (27), the ensemble-averaged TAMSD can be obtained as the coefficient of which is approximately equal to the one of (25).