Unexpected crossovers in correlated random-diffusivity processes

The passive and active motion of micron-sized tracer particles in crowded liquids and inside living biological cells is ubiquitously characterised by"viscoelastic"anomalous diffusion, in which the increments of the motion feature long-ranged negative and positive correlations. While viscoelastic anomalous diffusion is typically modelled by a Gaussian process with correlated increments, so-called fractional Gaussian noise, an increasing number of systems are reported, in which viscoelastic anomalous diffusion is paired with non-Gaussian displacement distributions. Following recent advances in Brownian yet non-Gaussian diffusion we here introduce and discuss several possible versions of random-diffusivity models with long-ranged correlations. While all these models show a crossover from non-Gaussian to Gaussian distributions beyond some correlation time, their mean squared displacements exhibit strikingly different behaviours: depending on the model crossovers from anomalous to normal diffusion are observed, as well as unexpected dependencies of the effective diffusion coefficient on the correlation exponent. Our observations of the strong non-universality of random-diffusivity viscoelastic anomalous diffusion are important for the analysis of experiments and a better understanding of the physical origins of"viscoelastic yet non-Gaussian"diffusion.

correlated Gaussian noise, we demonstrate that the similarity of these models in the Brownian case disappears in the anomalous diffusion case. We present detailed results for this non-universality in the viscoelastic anomalous diffusion case in terms of the time evolution of the MSDs, the effective diffusivities, and the PDFs of these processes. Specifically, we show that in some cases anomalous diffusion persists beyond the correlation time while in others normal diffusion emerges. Comparing our theoretical predictions with experiments will allow us to pinpoint more precisely the exact mechanisms of viscoelastic yet non-Gaussian diffusion with its high relevance to crowded liquids and live cells.

FBM-generalisation of the minimal diffusing-diffusivity model
We first analyse the FBM-generalisation of our minimal DD model [7], whose Langevin equation for the particle position reads dx/dt = 2D(t)ξ H (t) ( 1 ) in dimensionless form (see appendix A). The dynamics of D(t) is assumed to follow the square of an auxiliary Ornstein-Uhlenbeck process Y(t) [7], In the above ξ H (t) represents fractional Gaussian noise, understood as the derivative of smoothed FBM with zero mean and autocovariance ξ 2 H τ ≡ ξ H (t)ξ H (t + τ ) [35,36] decaying as ξ 2 H τ ∼ H(2H − 1)τ 2H−2 for τ longer than the physically infinitesimal (smoothening) time scale δ [35]. η(t) is a zero-mean white Gaussian noise of unit variance. We assume equilibrium initial conditions for Y(t), i.e., Y(0) is taken randomly from the equilibrium distribution f eq (Y) = π −1/2 exp(−Y 2 ) [7,17]. Thus the process Y(t) is stationary with variance Y 2 = D = 1/2. The autocorrelation is Y(t)Y(t + τ ) = exp(−|τ |)/2 with unit correlation time in our dimensionless units. From equation (1) we obtain the MSD (see appendix B) with kernel K(τ ) = √ D(t 1 )D(t 2 ) = (1/π)[b(τ ) + a(τ ) arctan(a(τ )/b(τ ))], where τ = |t 1 − t 2 |, a(τ ) = e −τ , and b(τ ) = √ 1 − a 2 (τ ). We first demonstrate how to get the main results for the MSD from simple estimates at short and long times compared to the correlation time of the D(t) dynamics. As the diffusion coefficient does not change considerably at times shorter than the correlation time, K(0) ≈ √ D(t)D(t) = D = 1/2, equation (4) yields For long times t 1, more care is needed: as we will see, the long-time limit is different for the persistent and anti-persistent cases. For the persistent case H > 1/2 we assume that the main contribution to the integral in equation (4) at long times comes from large τ , since the noise autocorrelation decays very slowly. We thus approximate K(τ ) ≈ |Y(t)| |Y(t + τ )| = |Y(t)| 2 = 1/π. Then, In the anti-persistent case H < 1/2 we split equation (4) into two integrals, 4t In the first integral at long times it is eligible to replace the upper limit of the integral by infinity, since it converges 9 . The second integral produces a subleading term, since it is bounded from above by Ct 2H , C being a constant. We therefore have the following asymptotic result for the MSD in the anti-persistent case at long times, with D eff = lim δ→0 2 +∞ 0 K(τ ) ξ 2 H τ dτ . Thus, the FBM-DD model demonstrates surprising crossovers in the behaviour of the MSD. In the persistent case the MSD scales as t 2H at both short and long times, but 9 If the diffusivity is constant, then K(τ ) is constant as well, and this approximation cannot be used, since necessarily   (1) and the exact MSD (C.13) for H = 1 as well as numerical integration of (4) for different H. The MSD approaches the limits (dashed lines) t 2H at short times and, at long times, anomalous [(2/π)t x2H ] or normal [2D eff t] scaling for super-and subdiffusion, respectively. Middle: effective diffusion coefficient as function of H. The theoretical curve [equation (D.10) for H < 1/2 and 1/π for H > 1/2] shows a distinct discontinuity at the Brownian value H = 1/2. Results from numerical evaluation of equations (D.1), (4), and simulations are shown to gradually approach the theoretical values (see text and appendix D). Right: crossover of the PDF from short-time non-Gaussian form with exponential tails to long-time Gaussian. The crossover is described in terms of the kurtosis (see figure E1).
with different diffusion coefficients. This is in a sharp contrast with the Brownian yet non-Gaussian diffusion characterised by the same, invariant diffusivity at all times. In the antipersistent case the situation is even more counterintuitive: the subdiffusive scaling of the MSD at short times crosses over to normal diffusion at long times.
The behaviour of the MSD is shown in figure 1. For superdiffusion, the change of the diffusivity between the short and long time superdiffusive scaling t 2H is distinct. Excellent agreement is observed between the exact and numerical evaluation for H = 1 and H = 0.7, 0.8, respectively. The exact analytical expression for H = 1 is derived in appendix C. In the subdiffusive case simulations and numerical evaluation nicely coincide and show the crossover from subdiffusion to normal diffusion. Figure 1 also shows the effective long time diffusivity. For superdiffusion the constant value 2/π ≈ 0.63 [see equation (6)] is distinct from the H-dependency for subdiffusion [H < 1/2, see equation (D. 10)]. For the Brownian case, D eff = 1/2, leading to a distinct discontinuity at H = 1/2. Note the slow convergence to the theory of simulations results and numerical evaluation of the respective integrals (see appendix D for details).
Given the above arguments that at short times (t 1) the diffusivity is approximately constant, we expect that in this regime the PDF corresponds to the superstatistical average of a single Gaussian over the stationary diffusivity distribution of the Ornstein-Uhlenbeck process, where x 2 (t) ST = t 2H and K 0 is the modified Bessel function of the second kind [7]. In the relevant large value limit of the scaling variable z = xt −H the Bessel function has the expansion K 0 (z) ∼ π/(2z) exp(−z) and thus represents the desired exponential tails, with a power-law correction [7]. 10 For long times (t 1) the diffusivity correlations decay and the Gaussian limit P(x, t) = G( x 2 (t) LT ) is recovered, where we introduce the general definition For H > 1/2, the long-time MSD is x 2 (t) LT = (2/π)t 2H while for H < 1/2, x 2 (t) LT = 2D eff t. The crossover behaviour of P(x, t) is indeed corroborated in figure 1 for different values of H. How do these observations compare to generalisations of other established random-diffusivity models? While in the normal-diffusive regime these models encode very similar behaviour, we show now that striking differences in the dynamics emerge when the motion is governed by long-range correlations.

FBM-generalisation of the Tyagi-Cherayil (TC) model
The generalisation of the Tyagi-Cherayil (TC) model [16] in dimensionless units reads 10 Such sub-dominant power-law corrections may indeed account for the deviations from the pure exponential shape of the PDF reported in [3]. However, many experimental data sets may not have sufficient resolution for smaller z values to pin down sub-dominant corrections. The exact result (E.14) gradually converges to the theoretical curve for different δ and t. Right: crossover of the PDF from short-time non-Gaussian shape with exponential tails to a long-time Gaussian. The crossover is described in terms of the kurtosis in appendix E.
This expression is obtained from the original equations (appendix E) via the transformations t → t/τ c and ). Using the same notation as before, η represents zero-mean white Gaussian noise and ξ H (t) is fractional Gaussian noise with Hurst exponent H.
The TC model looks quite similar to the minimal DD model as stochastically modulated Brownian motion, however, there exists a decisive difference: in equations (10) the OU-process Z(t) enters without the absolute value used in the minimal DD model (1). In expression (10) the prefactor Z(t) is therefore not a diffusion coefficient (by definition, a non-negative quantity). In the case H = 1/2, the analysis in [16] shows that on the level of the diffusion equation the quantity Z 2 (t) (in our notation here) takes on the role of the diffusion coefficient, and in this sense is thus well defined. The extension to fractional Gaussian noise therefore appears justified, yet we stress that the process (10) is intrinsically different from the FBM-DD model (1). As our discussion shows, the close similarity between the TS and DD models in the case H = 1/2 is replaced by a distinct dissimilarity in the emerging dynamics for H = 1/2.
The MSD of the FBM-TC model reads where the kernel K(τ ) is defined as It is shown in figure B1 along with the corresponding Langevin simulations. Before presenting the exact solution, let us apply an analogous reasoning for the behaviour of the MSD as developed for the FBM-DD model above. Namely, at short times we approximate K(τ ) ≈ Z 2 = 1/2. Then equation (11) At long times the MSD can be composed of the two parts x 2 (t) = 4t The upper limit of the first integral can be replaced by infinity because the first integral converges in both persistent and anti-persistent cases at long times [K(τ ) decays to 0 exponentially, different from the FBM-DD model]. The second term is subleading in comparison to the first term. As a result the MSD at long times scales linearly in time, H τ dτ . Indeed, from the exact form of the MSD in appendix E we obtain the limiting behaviours Thus for both sub-and superdiffusion this model shows a crossover from anomalous to normal diffusion, as demonstrated in figure 2. The effective long-time diffusion coefficient in this model varies continuously as D eff = Γ(2H + 1)/2 for all H. In particular, this means that for H = 1/2, D eff = 1/2. Figure 2 shows the exact match of the simulations results and the numerical evaluation at finite integration step. The PDF at short times coincides with the superstatistical limit in expression (8) above, as shown explicitly in equation (E.15). At long times we recover the Gaussian P(x, t) = G(Γ(2H + 1)t). Note that for H = 1 the noise is equal to unity at all times and the dynamics of x(t) is completely determined by the superstatistic encoded by the OU-process Z(t). The tails of the PDF are thus always exponential, reflected by the fact that the kurtosis has the invariant value 9 (see appendix E).
Despite the strong similarity between the DD and TC models in the Brownian case, for correlated driving noise their detailed behaviour is strikingly dissimilar, due to the different asymptotic forms of the kernel K(τ ) (figure B1).

FBM-generalisation of the switching (S) model
The third case model we consider here is the S-model with generalised noise [19], where n(t) is a two-state Markov chain switching between the values {0, 1} and ξ H (t) represents again fractional Gaussian noise. The constants D i are the diffusivities in the two states. The switching rates are k 12 and k 21 , such that the correlation time is τ c = 1/(k 12 + k 21 ). Note that the S-model (14) for white Gaussian noise with H = 1/2 is well known in the theory of stochastic processes [1,2]. In nuclear magnetic resonance literature it is known as the Kärger model [20,54]. From the first and second moments of the process θ(t), equations (F.5) and (F.6), we calculate the MSD of the process. In the Brownian limit H = 1/2 the MSD has a linear dependence at all times, This result was also obtained in [20]. For the general case with the correlation function based on fractional Gaussian noise, we have where At short times t τ c we find the scaling behaviour At long times t τ c the same scaling law is obtained, but with a different prefactor for the persistent case (H > 1/2), In contrast, for the anti-persistent case (H < 1/2), we derive a crossover to normal diffusion, From equations (15), (18), and (19), the long-time effective diffusivity can be obtained as The crossover behaviours of the MSD in the persistent and anti-persistent cases, analogous to the difference in the long-time scalings of the FBM-DD model, are displayed in figure 3. We also see some similarities between the FBM-S and FBM-DD models for the effective diffusivity. For the FBM-S model an H-dependent behaviour for H < 1/2 is followed by a discontinuity at H = 1/2 and then a constant value for H > 1/2. The results of the MSD for finite values δ and t are given in appendix F.
Next we discuss the PDF and kurtosis. At short times the continuous superstatistic of the previous cases is reduced to the discrete case of two superimposed Gaussians, producing the non-exponential form At long times a single Gaussian dominates, where x 2 (t) LT is given by equations (18) and (19) for the super-and subdiffusive cases, respectively. Figure 3 shows the superimposed two Gaussians at short times and the single Gaussian at long times.

Conclusions
Viscoelastic anomalous diffusion with long-ranged correlations is a non-Markovian, natively Gaussian process widely observed in complex liquids and the cytoplasm of biological cells. Most data analyses have concentrated on the MSD and the displacement autocorrelation function. Yet, once probed, the PDF in many of these systems turns out to be non-Gaussian, a phenomenon ascribed to the heterogeneity of the systems. Building on recent results for Brownian yet non-Gaussian diffusion, in which the non-Gaussian ensemble behaviour is understood as a consequence of a heterogeneous diffusivity coefficient, we here analysed three different random-diffusivity models driven by correlated Gaussian noise. Despite the simplicity of these models we observed surprising behaviours. Thus, while in the Brownian case all models display a linear MSD with invariant diffusion coefficient, in the correlated case a crossover occurs from short to long-time behaviours, with respect to the intrinsic correlation times. In particular, whether the long-time scaling of the MSD is anomalous or normal, depends on the specific model. Moreover, the effective diffusivity exhibits unexpectedly complex behaviours with discontinuities in the FBM-DD and FBM-S models. We note here that an additional crossover may come into play when a cutoff time scale of the fractional Gaussian noise becomes relevant [55].
In all cases a crossover from an initial non-Gaussian to a Gaussian PDF occurs. We showed that the FBM-S model is different from the other models in that it encodes an initial superposition of two Gaussians, turning into a single Gaussian at long times. We note that while the short-time exponential shape may point towards a universal, extreme-value jump-dominated dynamics [21], data also show stretched-Gaussian shapes [43], as well as long(er)-time convergence towards an exponential [8]. Clearly, the phenomenology of heterogeneous environments is rich and needs further investigation.
From an experimental point of view, the behaviours unveiled here may be used to explore further the relevance of the different possible stochastic formulations of random-diffusivity processes. For instance, in artificially crowded media one may vary the Hurst exponent by changing the volume fraction of crowders or the tracer sizes, or add drugs to change the system from super-to subdiffusive [30]. Comparison of the resulting scaling behaviours of MSD and associated effective diffusivity may then yield decisive clues.
The results found here will also be of interest in mathematical finance. In fact, the original DD model is equivalent to the Heston model [56] used to describe return dynamics of financial markets. Fractional Gaussian noise in mathematical finance is used to include an increased 'roughness' to the emerging dynamics [57]. The different models studied here could thus enrich market models.
The CLT is a central dogma in statistical physics, based on the fact that the entry variables are identically distributed. For inhomogeneous environments, ubiquitous in many complex systems, new concepts generalising the CLT will have to be developed. While random-diffusivity models are a start in this direction and provide relevant strategies for data analyses [58], ultimately more fundamental models including the quenched nature of the disordered environment [23,59] and extensions of models for non-equilibrium situations [60] need to be conceived.

Acknowledgments
We acknowledge funding from DFG (ME 1535/7-1). RM acknowledges the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP) for an Alexander von Humboldt Polish Honorary Research Scholarship. FS acknowledges Davide Straziota for helpful discussions and financial support of the 191017 BIRD-PRD project of the Department of Physics and Astronomy of Padua University. We acknowledge the support of the German Research Foundation (DFG) and the Open Access Publication Fund of Potsdam University.

Appendix A. Dimensionless units for the FBM-DD model
In dimensional form the starting equations governing the evolution of the position x(t) of the diffusing particle in the fractional version of the minimal DD-model read Here D(t) is the diffusion coefficient of dimension [D] = cm 2 s −1 , ξ H represents fractional Gaussian noise with the Hurst index H ∈ (0, 1] whose dimension is [ξ H ] = s H−1 and whose correlation function reads [35] ξ Noting that for the Gaussian noise sources we have ξ Now, we choose the temporal and spatial scales such that τ c = σ 2 = 1, to find With this choice of units, the stochastic equations of our minimal FBM-DD model are then given by equations (1) and (2) of the main text.

Using the integral
where erfc(z) = 1 − erf(z) = 2π −1/2 +∞ z e −t 2 dt is the complementary error function, we rewrite B 1 and B 2 as and Plugging equations (B.2) and (B.3) into (B.1) and after some transformations, we get which is equation (4) in the main text. This result is verified by simulation of the Ornstein-Uhlenbeck process in figure B1. We immediately obtain the first-order and second-order derivatives of K(τ ) with respect to τ , and K(τ ), K (τ ) and K (τ ) are all monotonic and have the following limits

Appendix E. FBM-generalisation of the Tyagi-Cherayil model
We now consider the fractional TC model Here ξ H (t) represents fractional Gaussian noise, η(t) is a white Gaussian noise, and the respective correlation functions are the same as in equation Equation (11) can be solved analytically, (E.5) and Considering the leading term of the Taylor expansion in terms of δ we get and After plugging equation (E.7) into (E.3) we get At long times t satisfying δ 1 t we have Here D eff can be calculated as For both persistence and anti-persistence cases, a crossover from anomalous diffusion to normal diffusion emerges. . (E.14) The second term on the right hand side contributes to the discrepancies near H → 0 in figure 2(b) of the main text. We expect the same behaviour of the PDF as for the DD model of reference [7] but with the rules of FBM. In particular, at short times we expect the superstatistical behaviour to hold and the PDF should be given by the weighted average of a single Gaussian over the stationary diffusivity distribution of the OU process. Therefore the expected PDF reads where G(σ 2 ) = (2πσ 2 ) −1/2 exp(−x 2 /(2σ 2 )) is the Gaussian distribution, p Z (Z) is the PDF of the dimensionless OU-process, and K 0 is the modified Bessel function of the second kind. At longer times the Gaussian limit will be reached, P(x, t) = G(Γ(2H + 1)t).
(E. 16) In particular, for H = 1, the PDF is always exponential at both short and long times. This can be seen from examination of the kurtosis, namely, the fourth order moment of the displacement reads  This means that for H = 1, the crossover to the Gaussian will never emerge at any time. This is a fundamental distinction from the FBM-DD model. The behaviour of the kurtosis is shown in figure E1. 2 ) 2 τ 2 c . The correlation (shown in figure B1 in comparison to Langevin simulations) approaches a 2 + a 2 at short times and a 2 at long times.
The fourth order moment of the displacement reads