Random diffusivity from stochastic equations: comparison of two models for Brownian yet non-Gaussian diffusion

The systematic study of the diffusive motion of tracer particles in fluids dates back to the 19th century, particularly referring to Robert Brown's experiments observing the erratic motion of granules extracted from pollen grains which were suspended in water [1]. Since then numerous scientists contributed by improving the experiments [2–4] as well as in defining the basis of the theory of diffusion [5–9]. In brief, Brownian or standard diffusion processes are mainly characterised by two central features: (i) the linear growth in time of the mean-squared displacement (MSD),

$$\begin{eqnarray}&&\langle {x}^{2}(t)\rangle =2{Dt},\end{eqnarray} \tag{ 1 }$$

where D is the diffusion coefficient, and (ii) the Gaussian probability density function (PDF) for the particle displacement,

$$\begin{eqnarray}&&G(x,t| D)=\displaystyle \frac{1}{\sqrt{4\pi {Dt}}}\exp \left(-\displaystyle \frac{{x}^{2}}{4{Dt}}\right).\end{eqnarray} \tag{ 2 }$$

Here and in the following we focus on a one-dimensional formulation of the model, a generalisation to higher dimensions can be achieved component-wise.

Discoveries of deviations from the linear time dependence (1) have a long history. Thus, Richardson already in 1926 reported his famed t-cubed law for the relative particle diffusion in turbulence [10]. Scher and Montroll uncovered anomalous diffusion of the power-law form

$$\begin{eqnarray}&&\langle {x}^{2}(t)\rangle \simeq {D}_{\alpha }{t}^{\alpha },\end{eqnarray} \tag{ 3 }$$

with the anomalous diffusion exponent 0 < α < 1 and the generalised diffusion coefficient D_α [11], for the motion of charge carriers in amorphous semiconductors [12]. With the advance of modern microscopy techniques, in particular, superresolution microscopy, as well as massive progress in supercomputing, anomalous diffusion of the type (3) has been reported in numerous complex and biological systems [13, 14]. Thus, subdiffusion with 0 < α < 1 was observed for submicron tracers in the crowded cytoplasm of biological cells [15–19] as well as in artificially crowded environments [20–23]. Further reports of subdiffusion come from the motion of proteins embedded in the membranes of living cells [24–26]. Subdiffusion is also seen in extensive simulations studies, for instance, of lipid bilayer membranes [27–30] and relative diffusion in proteins [31]. Superdiffusion, due to active motion of molecular motors, was observed in various biological cell types for both introduced and endogenous tracers [16, 17, 32, 33].

Most of the anomalous diffusion phenomena mentioned here belong to two main classes of anomalous diffusion: (i) the class of continuous time random walk processes, in which scale-free power-law waiting times in between motion events give rise to the law (3) [12, 34], along with a stretched Gaussian displacement probability density G(x, t) [11, 12, 34] as well as weak ergodicity breaking and ageing [35, 36]. We note that similar effects of non-Gaussianity, weak non-ergodicity, and ageing also occur in spatially heterogeneous diffusion processes [37–40]. (ii) The second one is the class of viscoelastic diffusion described by the generalised Langevin equation with power-law friction kernel [41, 42] and of fractional Brownian motion (FBM) [43]. These processes are both fuelled by long-range, power-law correlated noise. Its distribution is Gaussian, so that the displacement probability density G(x, t) is Gaussian, as well. Moreover, these are ergodic processes [23, 42, 44–46].

Over the last few years a new class of diffusive processes has been reported, namely, so-called Brownian yet non-Gaussian diffusion [47, 48]. This class identifies a dynamics characterised by a linear growth (1) of the MSD combined with a non-Gaussian PDF for the particle displacement. The emergence of a non-Gaussian distribution, despite the Brownian MSD scaling, suggests the presence of an inhomogeneity that can be located both on the single tracer particle and on the ensemble levels. The study of these processes is becoming increasingly relevant with the growing number of complex systems discovered to exhibit such statistical features. For instance, we mention soft matter and biological systems, in which the motion of biological macromolecules, proteins and viruses along lipid tubes and through actin networks [47, 48], as well as along membranes and inside colloidal suspension [49] and colloidal nanoparticles adsorbed at fluid interfaces [50–52] are studied. We also mention ecological processes, involving the characterisation of organism movement and dispersal [53, 54], as well as processes, that are Brownian but non-Gaussian in certain time windows of their dynamics. These concern the dynamics of disordered solids, such as glasses and supercooled liquids [55–57] as well as interfacial dynamics [58, 59]. Also anomalous diffusion processes of the viscoelastic class that typically are expected to exhibit Gaussian statistic of displacements, were reported to have non-Gaussian displacements along with distinct distributions of diffusivity values. These concern the motion of tracer particles in the cellular cytoplasm [60–62] and the motion of lipids and proteins in protein-crowded model membranes [29].

Here we study two alternative stochastic approaches to non-Gaussian diffusion due to random diffusivity parameters, namely, generalised grey Brownian motion (ggBM) and diffusing diffusivities (DD). We analyse their exact behaviour and relate these approaches to the idea of superstatistics. To prepare the discussion, section 2 presents a primer on the approach of superstatistics and what has been done in the context of ggBM and DD models. In section 3 we then study the ggBM model with a random diffusivity distributed according to the generalised Gamma distribution. In particular, ggBM will be shown to represent a stochastic description of the superstatistics approach and is equivalent to the short time (ST) limit of the DD model. In section 4 we formulate a set of stochastic equations for the dynamics within the DD framework, in which the diffusivity statistic is governed by the generalised Gamma distribution. This is then incorporated in the framework of the minimal model of DD in section 5. In section 5.4 we describe the behaviour of the kurtosis of the two models, an important quantity for data analysis. Section 6 introduces an analysis for an initial non-equilibrium setting for the random diffusivity, relevant, for instance, for the description of single particle trajectories. To transfer this concept to the ggBM approach we propose a non-equilibrium version of ggBM. Finally our conclusions are reported in section 7. In the appendices some mathematical details are collected.

When we talk about an ensemble of particles, we could imagine that non-Gaussian statistic in this ensemble sense emerges due to the fact that different particles are located in different environments with different transport characteristics, such as the diffusion coefficient. If during the observation time each particle remains in its own environment characterised by a given value D of the diffusivity, the ensemble of particles shows a mixture of individual Gaussians, weighted by some distribution p(D) of local diffusivities. This is the idea behind superstatistics, an approach promoted by Beck and Cohen [63], see also [64]. As a result, the ensemble dynamics is still Brownian yet the PDF of particle displacements will correspond to a sum or integral of single Gaussians with specific value of D, weighted by the distribution p(D). For instance, an exponential form for p(D) will produce an exponential shape of the ensemble displacement PDF, sometimes called a Laplace distribution. We note that there also exist superstatistical formulations on the basis of the stochastic Langevin equation, leading to Brownian yet non-Gaussian behaviour [65]. A quite general superstatistical formulation in terms of the gamma distribution was put forward by Hapca et al [53].

More recently, similar concepts have been sought to describe non-Gaussian viscoelastic subdiffusion. Thus, Lampo et al [61] observed exponential distributions of the generalised diffusivity D_α for the motion of submicron tracers in living bacteria and eukaryotic cells. As a theoretical description they used a superstatistical formulation of the stochastic equation for FBM [61]. Following the observation of stretched Gaussian shapes of the displacement PDF in protein-crowded lipid bilayer membranes [29], more general forms for the distribution of the generalised diffusion coefficient were introduced, see, for instance, [66, 67]. Viscoelastic, non-Gaussian diffusion was also described in terms of the generalised Langevin equation with superstatistical distribution of the friction amplitude [68, 69].

Some other models instead introduce a fluctuating diffusivity, for instance to describe segregation in solids [70] or to analyse data from diffusion processes assessed by modern measurement techniques [71]. Brownian motion in fluctuating environments, or governed by temperature or friction fluctuations has been studied in [72–74] and models with intermittency between two values of the diffusivity are considered in [75, 76]. Anomalous diffusion in a disordered system was also described in terms of a superstatistical model based on a Langevin equation formulation, combining a Rayleigh-shaped diffusivity distribution with deterministic power-law growth or decay of the mean diffusivity [77].

A general framework for the description of diffusion in a complex environment is provided also by the class of stochastic processes identified as ggBM [78–82]. The basic idea of this approach is that the complexity or heterogeneity of the medium is completely described by the random nature of a specific parameter. Choosing this parameter to be the diffusivity leads to a stochastic interpretation of the system that may be viewed as complementary to the superstatistics concept and thus suitable for the description of the class of Brownian yet non-Gaussian processes. We will define ggBM with a random diffusivity in more detail in the next section 3, and in the following demonstrate that ggBM is equivalent to the ST limit of the DD model.

Recently the idea of DD has received considerable attention. According to this approach, in addition to the introduction of a population of diffusivities, each particle during its motion is affected by a continuously changing diffusivity. Chubynsky and Slater first introduced this model describing the dynamics of the diffusion coefficient by a biased, stationary random walk with reflecting boundary conditions [83]. With this assumption the diffusivity changes slowly step by step, in the ST limit giving rise to normal diffusion with exponential displacement PDF⁶ . In the long time (LT) regime simulations showed a crossover to Gaussian diffusion with a single, effective diffusion coefficient [83]. In a more recent work a direct test of the DD mechanism for diffusion in inhomogeneous media is reported [86].

The DD concept was further studied by Jain and Sebastian [87, 88] and Chechkin et al [67]. While Jain and Sebastian use a path integral approach, Chechkin et al invoke the concept of subordination and an explicit exact solution for the PDF in Fourier space. Despite the different mathematical approach, both models recover the linear trend of the MSD and a distribution of displacements that at ST is exponential, while, at LT, it crosses over to a Gaussian with effective diffusivity, in agreement with the results in [83]. Tyagi and Cherayil [89] present a hybrid procedure between the two approaches, finding that the modulation of white noise by any stochastic process, whose time correlation function decays exponentially, is likely to have features similar to the ones obtained in [67, 83, 87, 88]. As a recent result we also report the work by Lanoiselée and Grebenkov in which the concept of DD is further investigated, for instance, with respect to time averages and ergodicity breaking properties [90].

In this paper we present a detailed comparison between the concept of ggBM with random diffusivity and the DD model. The main difference between the DD and ggBM model is represented by the interaction between environment and particles. On the one hand, in the DD model two different statistical levels are taken into account, one for the motion of the environment and one for the motion of the particles. The relation between these two gives rise to specific characteristics. Thus, at ST the slow variability of the environment guarantees the superstatistical limit. In the LT regime the diffusivity reaches a stationary average value leading the particles to develop a Gaussian statistic. On the other hand, the ggBM model does not directly involve an environment dynamics but only implies a dynamics in which the statistical features of the environment continuously drives the particles in their motion, see below for more details.

Concretely, for both ggBM and DD models a set of stochastic equations is introduced to generate a random diffusivity with a well defined stationary distribution. Until now mainly exponential or Gamma distributions have been considered for the random diffusivity. We here base the discussion on the generalised Gamma distribution, which represents an even broader class of distributions including the ones mentioned above, as particular cases. We define the generalised Gamma distribution by

$$\begin{eqnarray}&&{\gamma }_{\nu ,\eta }^{\mathrm{gen}}(D)=\displaystyle \frac{\eta }{{D}_{\star }^{\nu }\,{\rm{\Gamma }}(\nu /\eta )}{D}^{\nu -1}\exp \left(-{\left[\displaystyle \frac{D}{{D}_{\star }}\right]}^{\eta }\right),\end{eqnarray} \tag{ 4 }$$

where D_⋆ is a positive and dimensional constant and ν and η are positive constants. This distribution encodes the nth order stationary moments

$$\begin{eqnarray}&&\langle {D}^{n}{\rangle }_{\mathrm{stat}}={D}_{\star }^{n}\displaystyle \frac{{\rm{\Gamma }}([\nu +n]/\eta )}{{\rm{\Gamma }}(\nu /\eta )}.\end{eqnarray} \tag{ 5 }$$

The choice of the generalised Gamma distribution is based on experimental evidence demonstrating its role as a versatile description for generalised distributions in various complex systems. Indeed, in the context of superstatistics the generalised Gamma distribution was studied by Beck in [91]. Importantly, the generalised Gamma distribution includes those cases labelled as Gamma or exponential distribution that have already shown good agreement with several systems [53, 55–57]. Moreover it comprises the cases of stretched and compressed exponential distributions which may be useful for the interpretation of various systems [26, 53, 92, 93].

In the following we generalise the ggBM model from [78–82] to incorporate the generalised Gamma function (4). We then demonstrate how to reformulate the Ornstein–Uhlenbeck picture of the DD minimal model [67] and the closely related DD models [83, 87, 88] to include the distribution (4). With this extension both models are considerably more flexible for the description of measured data. Moreover, we will show that the ggBM model is a powerful stochastic representation of the superstatistics approach, and that the ggBM model equals the ST limit of the DD model. Finally, we consider non-equilibrium conditions in the DD model and propose a non-equilibrium extension of the ggBM model to consider similar effects in the stochastic setting of superstatistics. Such non-equilibrium initial conditions represent an important extension of the random diffusivity models, especially for experimentally relevant cases of single particle trajectory measurements.

GgBM is defined through the stochastic equation [78–82]

$$\begin{eqnarray}&&{X}_{\mathrm{ggBM}}(t)=\sqrt{2D}\times W(t),\end{eqnarray} \tag{ 6 }$$

for the particle trajectory ${X}_{\mathrm{ggBM}}(t)$, in which $W(t)={\int }_{0}^{t}\xi (t^{\prime} ){\rm{d}}{t}^{\prime}$ is standard Brownian motion, the Wiener process defined as the integral over the white Gaussian noise ξ(t) with zero mean. Moreover, D is a random diffusivity, here taken to be distributed according to the generalised Gamma distribution (4). The idea is that different, but physically identical particles move in disjointed environments, in which they experience different diffusivities, the essential view of the superstatistics approach. Alternatively, we could also think of physically different particles, with different diffusion coefficients, moving in an identical environment. The latter could, for instance, correspond to an ensemble of tracer beads with varying radius or different surface properties.

More mathematically speaking, ggBM is defined through the explicit construction of the underlying probability space based on self-similar increments, and it can be represented by the stochastic equation ${X}_{\mathrm{ggBM}}=\sqrt{{\rm{\Lambda }}}{X}_{g}$, where Λ is an independent, non-negative random variable, and X_g is a Gaussian process [78–82]. The characterisation of this class has also been studied for the case when X_g is a standard FBM and Λ is distributed according to the quite general class of M–Wright functions [81, 94]. We note that the definition (6) is similar to the superstatistical Langevin equation models in [65, 77].

Figure 1 shows trajectories obtained from direct simulations of the scheme (6), for which the diffusivity values D are chosen from the generalised Gamma distribution (4). As a result we obtain a Brownian motion characterised by a random amplitude, as demonstrated explicitly by the MSD plots for the same trajectories shown in the bottom panel of figure 1. For the value $\nu =1.5$ (right panels) larger D values are observed, in accordance with the shape of the distribution (4). The ggBM description is indeed close to the superstatistical concept and fundamentally different from the time evolution of the sample paths for the DD model, compare figure 7. However, at very ST both processes look much alike, as the DD model at ST will be shown to reduce to the ggBM model.

Zoom In Zoom Out Reset image size

The particle displacement distribution can be recovered following Pagnini and Paradisi [94]. If we define with Z₁ and Z₂ two real independent random variables whose PDFs are p₁(z₁) and p₂(z₂) with $-\infty \leqslant {z}_{1}\leqslant +\infty$ and $0\leqslant {z}_{2}\leqslant +\infty$, respectively, and with the random variable Z obtained by the product of Z₁ and Z₂^γ, that is, $Z={Z}_{1}{Z}_{2}^{\gamma }$, then, if we denote the PDF of Z with p(z), we find that

$$\begin{eqnarray}&&p(z)={\int }_{0}^{\infty }{p}_{1}\left(\displaystyle \frac{z}{{\lambda }^{\gamma }}\right){p}_{2}(\lambda )\displaystyle \frac{{\rm{d}}\lambda }{{\lambda }^{\gamma }}.\end{eqnarray} \tag{ 7 }$$

In the present case we identify X_ggBM(t), W(t), and the random diffusivity D with Z, Z₁, and Z₂, respectively. The PDF for the particle displacement encoded by equations (6) and (7) is given by

$$\begin{eqnarray}{f}_{\mathrm{ggBM}}(x,t) & = & {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{1}{\sqrt{2\pi t}}{\rm{\exp }}\left(-\displaystyle \frac{{(x/\sqrt{2D})}^{2}}{2t}\right){p}_{D}(D)\displaystyle \frac{{\rm{d}}{D}}{\sqrt{2D}}\\ & = & {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{1}{\sqrt{4\pi {Dt}}}{\rm{\exp }}\left(-\displaystyle \frac{{x}^{2}}{4{Dt}}\right){p}_{D}(D){\rm{d}}{D}\\ & = & {\displaystyle \int }_{0}^{\infty }G(x,t| D){p}_{D}(D){\rm{d}}{D},\end{eqnarray} \tag{ 8 }$$

where $G(x,t| D)$ is the Gaussian distribution (2) for given D. Such a representation of the PDF corresponds to the one of the superstatistical approach, proving the similarity of the two methods. The distribution p_D(D) is defined in (4) and the integral in (8), which can be solved exactly through different methods (appendix), provides the result (A.6) in terms of a Fox H-function (see appendix, where also the series expansion is given). The asymptotic behaviour of this result acquires the generalised exponential shape

$$\begin{eqnarray}&&{f}_{\mathrm{ggBM}}(x,t)\sim \displaystyle \frac{{({x}^{2}/[4{D}_{\star }t])}^{(2\nu -\eta -1)/(2[\eta +1])}}{{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi {D}_{\star }t}}\exp \left(-\displaystyle \frac{\eta +1}{\eta }{\eta }^{\tfrac{1}{\eta +1}}{\left[\displaystyle \frac{{x}^{2}}{4{D}_{\star }t}\right]}^{\eta /(\eta +1)}\right).\end{eqnarray} \tag{ 9 }$$

In particular, the choice η = 1 leads us back to exponential distributions, with power-law prefactor. Figure 2 demonstrates the agreement between the analytical result (9) for the PDF and the result of stochastic simulations of the underlying ggBM process, for different times and a fixed set of the parameters ν and η. In particular, we see that the shape of the distribution remains invariant—as for the superstatistical approach—and in contrast to the DD model analysed below.

Zoom In Zoom Out Reset image size

The MSD follows immediately from the following transformations,

$$\begin{eqnarray}\langle {X}_{{\rm{ggBM}}}^{2}(t)\rangle & = & {\displaystyle \int }_{-\infty }^{\infty }{x}^{2}{f}_{{\rm{ggBM}}}(x,t){\rm{d}}{x}\\ & = & {\displaystyle \int }_{0}^{\infty }{p}_{D}(D){\displaystyle \int }_{-\infty }^{\infty }{x}^{2}G(x,t| D){\rm{d}}{x}{\rm{d}}{D}\\ & = & {\displaystyle \int }_{0}^{\infty }{p}_{D}(D)2{D}{t}{\rm{d}}{D}=2t{\displaystyle \int }_{0}^{\infty }{{Dp}}_{D}(D){\rm{d}}{D}\\ & = & 2\langle D{\rangle }_{{\rm{stat}}}t,\end{eqnarray} \tag{ 10 }$$

where, according to (5), the effective diffusivity becomes

$$\begin{eqnarray}&&\langle D{\rangle }_{\mathrm{stat}}={D}_{\star }{\rm{\Gamma }}([\nu +1]/\eta )/{\rm{\Gamma }}(\nu /\eta ).\end{eqnarray} \tag{ 11 }$$

Figure 3 demonstrates the linearity of the variance. The fitted parameters are consistent with the model prediction, $\langle D{\rangle }_{\mathrm{stat}}=0.20$ comparing to the values chosen in the simulations.

Zoom In Zoom Out Reset image size

By means of the ggBM approach and with the introduction of a generalised Gamma distribution for the diffusivity we are able to reproduce a diffusive motion with a linear scaling of the MSD and a PDF characterised by a stretched or compressed Gaussian with a power-law prefactor. This is our first main result.

We now consider the diffusion coefficient D(t) to be a random function of time, defined by means of the auxiliary variable Y(t) through D(t) = Y²(t), similarly to the DD minimal model introduced earlier [67]. Our goal is to construct a stochastic equation for the additional variable Y(t) such that the stationary PDF for its square is the generalised Gamma distribution in (4). Thus, our present model is represented by the following set of stochastic equations

$$\begin{eqnarray}&&{\rm{d}}{Y}=a(Y){\rm{d}}{t}+\sigma \times {\rm{d}}{W}(t)\end{eqnarray} \tag{ 12a }$$

$$\begin{eqnarray}&&D(t)={Y}^{2}(t),\end{eqnarray} \tag{ 12b }$$

where $a(Y)$ is a nonlinear function whose explicit shape is obtained below, $\sigma$ is a constant and W(t) is a Wiener process with variance $\langle {W}^{2}(t)\rangle =t$. The physical dimension of the auxiliary variable is $[Y]=\mathrm{cm}\,{{\rm{s}}}^{-1/2}$ and for the constant $\sigma$ we have $[\sigma ]=\mathrm{cm}\,{{\rm{s}}}^{-1}$.

Our approach is based on the central idea that it is possible to establish a direct relation between the PDFs of the two variables Y(t) and D(t). This allows us to introduce a completely new dynamics for the auxiliary variable. Such a dynamics, even though more complex, allows to reproduce a more general class of PDFs for the random diffusivity and thus provides a significant extension of the DD model, which will be our second main result.

To proceed we set p(Y, t) to represent the PDF of the process Y(t) described in (12a). It fulfils the Fokker–Plank equation [9]

$$\begin{eqnarray}&&\displaystyle \frac{\partial p(Y,t)}{\partial t}+\displaystyle \frac{\partial a(Y)p(Y,t)}{\partial Y}=\displaystyle \frac{{\sigma }^{2}}{2}\displaystyle \frac{{\partial }^{2}p(Y,t)}{\partial {Y}^{2}}.\end{eqnarray} \tag{ 13 }$$

Considering the stationary situation the corresponding time independent PDF p_Y(Y) fulfils the equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial a(Y){p}_{Y}(Y)}{\partial Y}=\displaystyle \frac{{\sigma }^{2}}{2}\displaystyle \frac{{\partial }^{2}{p}_{Y}(Y)}{\partial {Y}^{2}},\end{eqnarray} \tag{ 14 }$$

from which we infer the relation

$$\begin{eqnarray}&&a(Y)=\displaystyle \frac{{\sigma }^{2}}{2{p}_{Y}(Y)}\displaystyle \frac{\partial {p}_{Y}(Y)}{\partial Y},\end{eqnarray} \tag{ 15 }$$

directly relating the drift coefficient a(Y) with the stationary distribution of Y(t) [95].

We then recall that, given two random variables Z₁ and Z₂ related by Z₂ = g(Z₁), for appropriate functions g(z) we have [96]

$$\begin{eqnarray}&&{p}_{{Z}_{2}}({z}_{2})={p}_{{Z}_{1}}({g}^{-1}({z}_{2}))\left|\displaystyle \frac{{\rm{d}}}{{{\rm{d}}{z}}_{2}}{g}^{-1}({z}_{2})\right|.\end{eqnarray} \tag{ 16 }$$

This implies that the distributions of the variables Y(t) and D(t) are related via

$$\begin{eqnarray}&&{p}_{Y}(Y,t)=| Y| {p}_{D}({Y}^{2},t).\end{eqnarray} \tag{ 17 }$$

Based on this we construct a set of stochastic equations for the desired quantity D(t). Starting from the chosen stationary distribution p_D(D) of the random diffusivity we define the stationary distribution p_Y(Y) for the auxiliary variable Y(t) by means of equation (17). Finally relation (15) allows us to recover the suitable coefficient a(Y) in equation (12a). Following the described scheme for the generalised Gamma distribution (4) we obtain

$$\begin{eqnarray}&&{p}_{Y}(Y)=| Y| \displaystyle \frac{\eta }{{D}_{\star }^{\nu }{\rm{\Gamma }}(\nu /\eta )}{Y}^{2(\nu -1)}\exp \left(-{\left[\displaystyle \frac{Y}{\sqrt{{D}_{\star }}}\right]}^{2\eta }\right),\end{eqnarray} \tag{ 18 }$$

and thus

$$\begin{eqnarray}&&\displaystyle \frac{\partial {p}_{Y}(Y)}{\partial Y}(Y)=\displaystyle \frac{\eta \mathrm{sgn}(Y)}{{D}_{\star }^{\nu }{\rm{\Gamma }}(\nu /\eta )}{Y}^{2(\nu -1)}\exp \left(-{\left[\displaystyle \frac{Y}{\sqrt{{D}_{\star }}}\right]}^{2\eta }\right)\left(2\nu -1-2\eta {\left[\displaystyle \frac{Y}{\sqrt{{D}_{\star }}}\right]}^{2\eta }\right).\end{eqnarray} \tag{ 19 }$$

This finally leads us to the desired drift coefficient

$$\begin{eqnarray}&&a(Y)=\displaystyle \frac{{\sigma }^{2}}{2Y}\left(2\nu -1-2\eta {\left[\displaystyle \frac{Y}{\sqrt{{D}_{\star }}}\right]}^{2\eta }\right).\end{eqnarray} \tag{ 20 }$$

The stochastic equations (12a) together with the explicit form (20) of the drift coefficient for the diffusivity fluctuations provide a complete and generalised analogue of the DD model, which is extremely flexible for the modelling of experimental data.

We notice that in the particular case of ν = 0.5 and η = 1 we recover the Ornstein–Uhlenbeck model (diffusion in an harmonic potential) considered in the original minimal DD model [67]. As already remarked in [67] in this setting the resulting stochastic equation for D(t) is nothing else than the Heston model, that is widely used in financial mathematics and specifies the time evolution of the stochastic volatility of a given asset [90, 97, 98].

Equation (12a) can be readily solved numerically with initial conditions taken randomly from the equilibrium distribution (18). Figures 4 and 5 show sample time evolutions of the auxiliary variable Y and the diffusivity D = Y² for the DD process based on the steady state generalised Gamma distribution, as obtained below. We note that for the case ν = 0.5 in figure 4 the sample paths of the variable Y(t) frequently cross the zero line, while for the case ν = 1.5 in figure 5 the zero line is avoided, corresponding to the uni- and bimodal shapes of the PDFs of the variable Y(t) evaluated in figure 6. The existence of a pole in the generalised Gamma distribution (4) at D = 0 for the case ν = 0.5 thus creates a very different behaviour than for the case ν = 1.5 without singularity. For the diffusivity variable D(t) in figures 4 and 5 the regions of Y(t) close to the zero line lead to smaller D(t) values in the same regions. Finally, figures 4 and 5 demonstrate the exponential shape of the autocorrelation functions for both Y(t) and D(t),

$$\begin{eqnarray}&&{\mathrm{ACF}}_{Y}(t,t^{\prime} )=\langle (Y(t)-\langle Y\rangle )(Y(t+t^{\prime} )-\langle Y\rangle )\rangle \end{eqnarray} \tag{ 21 }$$

and an analogous expression for $D(t)$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We know from previous studies of DD models that the correlation time of the random diffusivity represents a key factor in the study of the particle dynamics. The correlation time τ_c is evaluated both by means of a two-parametric numerical fit to the exponential function and through the integral

$$\begin{eqnarray}&&{\tau }_{c}\sim \displaystyle \frac{1}{{\rm{ACF}}(0)}{\int }_{0}^{\infty }{\rm{ACF}}(\tau ){\rm{d}}\tau ,\end{eqnarray} \tag{ 22 }$$

which is exact for pure exponential autocorrelation functions. The results obtained by the two methods are reported in figure 4 and 5 and they are in excellent agreement, from which we conclude that the diffusivity autocorrelation is exponential to leading order and thus the correlation time τ_c well defined.

It is interesting to notice that the auxiliary function $Y(t)$ in the case of a bimodal distribution possesses a non-zero correlation function in the stationary state. This is due to the fact that despite a vanishing global mean of the PDF, depending on the initial setting each trajectory is representative of only one side of the bimodal PDF.

With the set of equations defined in section 4 we can consider the generalisation of the DD minimal model described in [67], and obtain the process in position space, ${X}_{\mathrm{DD}}(t)$. Recalling the idea of introducing an analytic description for the dynamics of the random diffusivity, we take that the motion of the particle is defined by the integral version of the overdamped Langevin equation,

$$\begin{eqnarray}&&{X}_{\mathrm{DD}}(t)={\int }_{0}^{t}\sqrt{2D(t^{\prime} )}\times \xi (t^{\prime} ){\rm{d}}{t}^{\prime} ,\end{eqnarray} \tag{ 23 }$$

where $\xi (t)$ is white Gaussian noise and D(t) is the random time-dependent diffusivity obtained in section 4. This dynamics based on the above results for the diffusivity dynamics generalises the idea introduced in [67], where an Ornstein–Uhlenbeck process was selected for the auxiliary variable. Figure 7 shows trajectories obtained from the stochastic equation (23) where the diffusivity was generated from (12a) with initial conditions taken randomly from the stationary distribution. In ggBM each trajectory has the same D value, while in the DD model the value of D changes as function of time. In turn, individual trajectories of the DD model are quite similar.

Zoom In Zoom Out Reset image size

Since the DD model is a direct generalisation of the minimal DD model we expect a crossover to a Gaussian displacement PDF for times longer than the correlation time τ_c. We thus carry on our analysis for the ST and LT regimes separately, before analysing the MSD and kurtosis of this DD process.

5.1. Short time regime

Since the dynamics of the environment is determined by the correlation time τ_c we expect that on ST scales with $t\ll {\tau }_{c}$ the diffusion coefficient is approximately fixed for each particle and we thus suppose the validity of a superstatistical description at ST,

$$\begin{eqnarray}&&{X}_{\mathrm{DD}}^{\mathrm{ST}}\sim \sqrt{2D}{\int }_{0}^{t}\xi (t^{\prime} ){\rm{d}}{t}^{\prime} =\sqrt{2D}\times W(t).\end{eqnarray} \tag{ 24 }$$

The existence of the superstatistical regime at $t\ll {\tau }_{c}$ is consistent with the model considered in [67] and with the results reported in [89] concerning the modulation of white noise by any stochastic process whose time correlation function decays exponentially. The superstatistical approach allows us to estimate the ST distribution of the particle displacement by means of

$$\begin{eqnarray}&&{f}_{\mathrm{DD}}^{\mathrm{ST}}(x,t)\sim {\int }_{0}^{\infty }{p}_{D}(D)G(x,t| D){\rm{d}}{D}.\end{eqnarray} \tag{ 25 }$$

This representation corresponds to the ggBM scenario established above, which means that we can borrow its results in equations (A.6) and (9), considering that ${f}_{\mathrm{DD}}^{\mathrm{ST}}(x,t)\sim {f}_{\mathrm{ggBM}}(x,t)$.

The expected behaviour (9) is confirmed by extensive numerical simulations. Figures 8(a) and 9(a) show the ST PDFs for two different sets of the parameters ν and η, and in both cases we observe excellent agreement with the asymptotic behaviour (9).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Comparing figure 2 with figure 8(a) we notice that the ggBM model allows one to describe a process that preserves the exact non-Gaussian PDF, which is exactly the same PDF we obtain in the DD model in the ST regime. Both approaches describe the same superstatistical frame but the DD model then crosses over to a Gaussian beyond the correlation time τ_c, see below the discussion of the kurtosis. The establishment of the relation between the DD model and the previously devised ggBM at ST is our third main result.

5.2. Long time regime

At LT, again taking our clue from [67] and from the general results in [89], we expect that eventually a crossover to a Gaussian distribution is observed (as already anticipated in figures 8 and 9). Above the correlation time, that is, for times $t\gg {\tau }_{c}$ we thus look for a PDF given by

$$\begin{eqnarray}&&{f}_{\mathrm{DD}}^{\mathrm{LT}}(x,t)\sim \displaystyle \frac{1}{\sqrt{4\pi \langle D{\rangle }_{\mathrm{stat}}t}}\exp \left(-\displaystyle \frac{{x}^{2}}{4\langle D{\rangle }_{\mathrm{stat}}t}\right),\end{eqnarray} \tag{ 26 }$$

with the effective diffusivity (11). The numerical results reported in figures 8(b) and 9(b) prove the validity of this behaviour. At sufficient LT the particles have explored all the diffusivity space and a Gaussian behaviour with an effective diffusivity emerges. This leads to a standard Brownian diffusive behaviour. We stress again that the transition from a non-Gaussian to a Gaussian profile depends on the value of the correlation time τ_c of the diffusivity process.

5.3. Mean squared displacement

For the DD model we found a crossover of the PDF of the spreading particles. An initial non-Gaussian behaviour is slowly replaced by a Gaussian one. The superstatistical behaviour of the DD approach at ST is equivalent to the ggBM model and is characterised by the non-Gaussianity. Nevertheless, as expected from previous studies [67], the MSD does not change in the course of time and is the same at ST and LT regimes. Direct calculation indeed produces the invariant form

$$\begin{eqnarray}&&\langle {X}_{\mathrm{DD}}^{2}(t)\rangle =2\langle D{\rangle }_{\mathrm{stat}}t.\end{eqnarray} \tag{ 27 }$$

This continuity of the MSD is demonstrated in figure 10, together with a linear fit proving the validity of the linear trend.

Zoom In Zoom Out Reset image size

5.4. Kurtosis

In what follows the second and fourth moments of the non-Gaussian PDF identified in equations (8) and (25) are studied in terms of the kurtosis that represents one of the first checks for non-Gaussianity. We recall the second order moment calculated in (10) and in a similar way we obtain the fourth order moment

$$\begin{eqnarray}&&\langle {X}_{\mathrm{ggBM}}^{4}(t)\rangle =\langle {X}_{\mathrm{DD}}^{4}(t){\rangle }_{\mathrm{ST}}=12\langle {D}^{2}{\rangle }_{{\rm{stat}}}{t}^{2},\end{eqnarray} \tag{ 28 }$$

where $\langle {D}^{2}{\rangle }_{{\rm{stat}}}$ is the second moment of the diffusivity in the stationary state. By means of results (10) and (28) and recalling the definition of the diffusivity moments in equation (5), the kurtosis $K=\langle {x}^{4}(t)\rangle /\langle {x}^{2}(t){\rangle }^{2}$ is given by

$$\begin{eqnarray}{K}_{\mathrm{ggBM}} & = & {K}_{\mathrm{DD}}^{\mathrm{ST}}=4{t}^{2}{D}_{* }^{2}\displaystyle \frac{3{\rm{\Gamma }}([\nu +2]/\eta )}{{\rm{\Gamma }}(\nu /\eta )}\times \displaystyle \frac{1}{4{t}^{2}{D}_{* }^{2}}{\left(\displaystyle \frac{{\rm{\Gamma }}(\nu /\eta )}{{\rm{\Gamma }}([\nu +1]/\eta )}\right)}^{2}\\ & = & 3\displaystyle \frac{{\rm{\Gamma }}([\nu +2]/\eta ){\rm{\Gamma }}(\nu /\eta )}{{\rm{\Gamma }}{([\nu +1]/\eta )}^{2}}\gt 3,\end{eqnarray} \tag{ 29 }$$

for ggBM and the short-time DD process. The non-Gaussian PDF represents a leptokurtic behaviour as can be observed in figure 11, showing the kurtosis of the DD and ggBM models. The value for the kurtosis at ST is in agreement with the value reported in (29). At LT the DD kurtosis approaches the value 3 characteristic of the Gaussian distribution, while the ggBM one keeps fluctuating around the same initial value.

Zoom In Zoom Out Reset image size

The results discussed above consider equilibrium initial conditions for the diffusivity fluctuations. In particular, results (10) and (27) for the particle MSD exhibit the invariant form $\langle {X}^{2}(t)\rangle =2\langle D{\rangle }_{\mathrm{stat}}t$ in both cases. Such equilibrium initial conditions will in general not be fulfilled for particles that are initially seeded in a non-equilibrium environment. For instance, in single particle tracking a tracer bead can be introduced into the system at t = 0, or similar in computer simulations. After this disturbance the environment equilibrates again. To accommodate for such a case we here study a minimal model for the case of non-equilibrium initial conditions, which leads to another main result of this work. As we are going to see, this non-equilibrium scenario gives rise to differences in the characteristics of the two studied models. In particular, we observe an initial ballistic behaviour. The LT behaviour, of course, does not show differences since in this range the diffusivity reaches its stationary state and we can again consider the results obtained in the previous sections for the LT limit.

We illustrate the role of non-equilibrium conditions by taking a specific, and in fact the simplest, set of parameters, ν = 0.5 and η = 1. This defines the stochastic dynamical equation in (12a) as

$$\begin{eqnarray}&&D(t)={Y}^{2}(t)\end{eqnarray} \tag{ 30a }$$

$$\begin{eqnarray}&&{\rm{d}}{Y}=-\displaystyle \frac{{\sigma }^{2}}{{D}_{\star }}{Y}{\rm{d}}{t}+\sigma {\rm{d}}{W}(t),\end{eqnarray} \tag{ 30b }$$

that corresponds to the well known dynamics of the Ornstein–Uhlenbeck process for the study in [67] with the correlation time ${\tau }_{c}={D}_{\star }/{\sigma }^{2}$. We start considering the related Fokker–Planck equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial }{\partial t}p(Y,t)=\displaystyle \frac{{\sigma }^{2}}{{D}_{\star }}\displaystyle \frac{\partial }{\partial Y}{Yp}(Y,t)+\displaystyle \frac{{\sigma }^{2}}{2}\displaystyle \frac{{\partial }^{2}}{\partial {Y}^{2}}p(Y,t).\end{eqnarray} \tag{ 31 }$$

We can solve this equation with a non-equilibrium condition, for instance, p(Y, 0) = δ(Y − Y₀), using the method of characteristics in Fourier space. We readily derive the general solution

$$\begin{eqnarray}&&p(Y,t| {Y}_{0})={(\pi {D}_{\star }[1-\exp (-2t{\sigma }^{2}/{D}_{\star })])}^{-1/2}\exp \left(-\displaystyle \frac{{(Y-{Y}_{0}\exp (-t{\sigma }^{2}/{D}_{\star }))}^{2}}{{D}_{\star }(1-\exp (-2t{\sigma }^{2}/{D}_{\star }))}\right).\end{eqnarray} \tag{ 32 }$$

Recalling relation (16) for the diffusivity PDF we then obtain

$$\begin{eqnarray}{p}_{D}(D,t| {D}_{0}={Y}_{0}^{2}) & = & \displaystyle \frac{1}{2\sqrt{D}}[p(\sqrt{D},t)+p(-\sqrt{D},t)]\\ & = & {(4\pi {D}_{\star }D[1-\exp (-2t{\sigma }^{2}/{D}_{\star })])}^{-1/2}\\ & & \times \left\{\exp \left(-\displaystyle \frac{{(\sqrt{D}-\sqrt{{D}_{0}}\exp (-t{\sigma }^{2}/{D}_{\star }))}^{2}}{{D}_{\star }(1-\exp (-2t{\sigma }^{2}/{D}_{\star }))}\right)\right.\\ & & \left.+\exp \left(-\displaystyle \frac{{(-\sqrt{D}-\sqrt{{D}_{0}}\exp (-t{\sigma }^{2}/{D}_{\star }))}^{2}}{{D}_{\star }(1-\exp (-2t{\sigma }^{2}/{D}_{\star }))}\right)\right\}.\end{eqnarray} \tag{ 33 }$$

We point out that in the limit of LT this result provides exactly the stationary distribution described in (4) with the specific set of parameters defined above. This is also verified by the trend of the average value

$$\begin{eqnarray}&&\langle D(t)\rangle =\displaystyle \frac{1}{2}({D}_{\star }(1-{{\rm{e}}}^{-2t{\sigma }^{2}/{D}_{\star }})+2{D}_{0}{{\rm{e}}}^{-2t{\sigma }^{2}/{D}_{\star }}),\end{eqnarray} \tag{ 34 }$$

in agreement with result (4).

In contrast to the previous analysis, we observe an explicit dependence on time of p_D(D, t), which makes the calculations more involved. Thus, we select an initial condition for the diffusivity, D₀ = 0, which is convenient for the study of the particles displacement distribution. This leads to a reduction in (33), namely,

$$\begin{eqnarray}&&p(D,t| {D}_{0}=0)={(\pi {D}_{\star }D(1-\exp (-2t{\sigma }^{2}/{D}_{\star })))}^{-1/2}\exp \left(\displaystyle \frac{D}{{D}_{\star }(1-\exp (-2t{\sigma }^{2}/{D}_{\star }))}\right).\end{eqnarray} \tag{ 35 }$$

We now study the two models in this particular case of a non-equilibrium initial condition for the diffusivity.

6.1. Diffusing diffusivities with non-equilibrium initial diffusivity condition

The dynamics for the diffusivity encoded in equations (30a) and (30b) when choosing the specific set of parameters ν = 0.5 and η = 1 is the same as described in [67] when d = n = 1. Thus, in this paragraph, we extend the description of the minimal DD model studied in [67] to the case of non-equilibrium initial conditions for the diffusivity. In order to proceed with the same notation we introduce dimensionless units for relations (30a) and (30b) as well as for the overdamped Langevin equation describing the particle motion [67], such that the full set of stochastic equations reads

$$\begin{eqnarray}{X}_{\mathrm{DD}} & = & {\displaystyle \int }_{0}^{t}\sqrt{2D(t^{\prime} )}\xi (t^{\prime} ){\rm{d}}{t}^{\prime} \\ D(t) & = & {Y}^{2}(t)\\ {\rm{d}}{Y} & = & -{Y}{\rm{d}}{t}+{\rm{d}}{W}(t).\end{eqnarray} \tag{ 36 }$$

A subordination approach can then be used to obtain the distribution of the particle displacement [67], namely,

$$\begin{eqnarray}&&{f}_{{\rm{DD}}}(x,t)={\int }_{0}^{\infty }T(\tau ,t)G(x,\tau ){\rm{d}}\tau ,\end{eqnarray} \tag{ 37 }$$

where G(x, τ) is the Gaussian (2) and T(τ, t) represents the PDF of the process $\tau (t)={\int }_{0}^{t}{Y}^{2}(t^{\prime} ){\rm{d}}{t}^{\prime}$. Starting from the subordination formula (37) we obtain the relation

$$\begin{eqnarray}&&{\hat{f}}_{{\rm{DD}}}(k,t)=\tilde{T}(s={k}^{2},t)\end{eqnarray} \tag{ 38 }$$

where with the symbols $\hat{\cdot }$ and $\tilde{\cdot }$ we indicate the Fourier and Laplace transforms, respectively. For the particular initial condition D₀ = 0, which is equivalent to y₀ = 0, the solution is known [99, 100],

$$\begin{eqnarray}&&\tilde{T}(s,t)=\exp \left(\displaystyle \frac{t}{2}\right)/{\left(\displaystyle \frac{1}{\sqrt{1+2s}}\sinh (t\sqrt{1+2s})+\cosh (t\sqrt{1+2s})\right)}^{1/2}.\end{eqnarray} \tag{ 39 }$$

This latter quantity is directly related to the MSD of the particles through [67]

$$\begin{eqnarray}&&\langle {X}_{\mathrm{DD}}^{2}(t)\rangle =-2{\left.\displaystyle \frac{\partial \tilde{T}(s,t)}{\partial s}\right|}_{s=0}.\end{eqnarray} \tag{ 40 }$$

We readily obtain the closed form result

$$\begin{eqnarray} \langle {X}_{{\rm{DD}}}^{2}(t) \rangle =t-\displaystyle \frac{1}{2}(1-{{\rm{e}}}^{-2t})\sim \left\{\begin{array}{ll}{t}^{2}, & t\ll 1\\ 2t, & t\gg 1\end{array}\right..\end{eqnarray} \tag{ 41 }$$

The resulting dynamics is thus no longer Brownian at all times. In contrast, at times shorter than the correlation time (in the dimensionless units used here ${\tau }_{c}=1$) we obtain a ballistic scaling of the MSD. This behaviour reflects the fact that the diffusivity equilibration in this case with D₀ = 0 leads to an initial acceleration.

Starting from equations (38) and (39) we consider approximations of the PDF for ST and LT which, since we are in dimensionless units, correspond to $t\ll 1$ and $t\gg 1$ respectively. In the ST limit, the Fourier transform of the PDF becomes

$$\begin{eqnarray}&&{\hat{f}}_{{\rm{DD}}}^{{\rm{ST}}}(k,t)\sim \displaystyle \frac{{(1+t+{t}^{2}/2)}^{1/2}}{{(1+t+{t}^{2}/2+{k}^{2}{t}^{2})}^{1/2}}\sim \displaystyle \frac{1}{t}{\left({k}^{2}+\displaystyle \frac{1}{{t}^{2}}\right)}^{-1/2}.\end{eqnarray} \tag{ 42 }$$

Note that this expression is normalised, ${\hat{f}}_{{\rm{DD}}}(k=0,t)=1$. After taking the inverse Fourier transform we find

$$\begin{eqnarray}{f}_{{\rm{DD}}}^{{\rm{ST}}}(x,t) & \sim & \displaystyle \frac{1}{\pi t}{\displaystyle \int }_{0}^{\infty }\displaystyle \frac{\cos ({kx})}{{({k}^{2}+1/{t}^{2})}^{1/2}}{\rm{d}}{k}\\ & = & \displaystyle \frac{1}{\pi t}{K}_{0}\left(\displaystyle \frac{| x| }{t}\right).\end{eqnarray} \tag{ 43 }$$

Re-establishing dimensional units, this result becomes

$$\begin{eqnarray}&&{f}_{{\rm{DD}}}^{{\rm{ST}}}(x,t)\sim \displaystyle \frac{1}{\pi \sigma t}{K}_{0}\left(\displaystyle \frac{| x| }{\sigma t}\right).\end{eqnarray} \tag{ 44 }$$

Here K_ν(z) is the modified Bessel function of second type. The asymptotic behaviour of this distribution for $| x| \to \infty$ is the Laplace distribution

$$\begin{eqnarray}&&{f}_{{\rm{DD}}}^{{\rm{ST}}}(x,t)\sim \displaystyle \frac{1}{\sqrt{2\pi \sigma t| x| }}\exp \left(-\displaystyle \frac{| x| }{\sigma t}\right).\end{eqnarray} \tag{ 45 }$$

In the LT limit equations (38) and (39) yield

$$\begin{eqnarray}&&{\hat{f}}_{{\rm{DD}}}^{{\rm{LT}}}(k,t)\sim \displaystyle \frac{{2}^{1/2}\exp (t[1-\sqrt{1+2{k}^{2}}]/2)}{{(1+1/\sqrt{1+2{k}^{2}})}^{1/2}},\end{eqnarray} \tag{ 46 }$$

that again is normalised. If we focus on the tails of the distribution in the limit $k\ll 1$ we obtain the Gaussian

$$\begin{eqnarray}&&{\hat{f}}_{{\rm{DD}}}^{{\rm{LT}}}(k,t)\sim \exp (-{k}^{2}t/2)\end{eqnarray} \tag{ 47 }$$

in Fourier space, corresponding to the Gaussian

$$\begin{eqnarray}&&{f}_{{\rm{DD}}}^{{\rm{LT}}}(x,t)\sim \displaystyle \frac{1}{\sqrt{2\pi t}}\exp \left(-\displaystyle \frac{{x}^{2}}{2t}\right)\end{eqnarray} \tag{ 48 }$$

in direct space. Restoring dimensional units and recalling that $\langle D{\rangle }_{\mathrm{stat}}={D}_{\star }/2$, eventually provides

$$\begin{eqnarray}{f}_{{\rm{DD}}}^{{\rm{LT}}}(x,t) & \sim & \displaystyle \frac{1}{\sqrt{2\pi {D}_{\star }t}}{\rm{\exp }}\left(-\displaystyle \frac{{x}^{2}}{2{D}_{\star }t}\right)\\ & = & \displaystyle \frac{1}{\sqrt{4\pi \langle D{\rangle }_{{\rm{stat}}}t}}{\rm{\exp }}\left(-\displaystyle \frac{{x}^{2}}{4\langle D{\rangle }_{{\rm{stat}}}t}\right),\end{eqnarray} \tag{ 49 }$$

where in the last step we identified the equilibrium value $\langle D{\rangle }_{\mathrm{stat}}$ of the diffusivity. From the approximations (45) and (49) we readily recover the two limiting scaling laws for the variance in equation (41).

Figure 12 nicely corroborates these findings, comparing the non-equilibrium DD model results for the PDF obtained above with results from stochastic simulations. The crossover behaviour of the associated MSD is displayed in figure 13, again showing very good agreement with the theory.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

6.2. Non-equilibrium ggBM

The ggBM model discussed in section 3 is based on the static distribution p_D(D) of the diffusivity. In order to explore non-equilibrium effects as discussed above for the DD model also within the superstatistical approach, we here propose a non-equilibrium generalisation of the ggBM model. Thus, we generalise the standard ggBM definition (6) and introduce a variability of D in time, according to the stochastic equation

$$\begin{eqnarray}&&{X}_{\mathrm{ggBM}}(t)=\sqrt{2D(t)}\times W(t).\end{eqnarray} \tag{ 50 }$$

Physically, this new concept may be interpreted as fluctuations of the disjointed environments experienced by the different particles or to temporal changes of the particle size, for instance, due to agglomeration-separation dynamics.

Based on the definition (50) it is then straightforward to take the dynamics of $D(t)$ to be the same as the one considered for the DD model. This guarantees that the ensemble properties of this generalised process (50) are exactly the same as the ones of the standard ggBM model studied in section 3. In particular, the dependence on time of the diffusivity does not affect the validity of equation (7), so in order to estimate the PDF of the particle displacement of the ggBM model, we consider the distribution (35) in the calculation of the integral

$$\begin{eqnarray}&&{f}_{\mathrm{ggBM}}(x,t| {D}_{0})={\int }_{0}^{\infty }{p}_{D}(D,t| {D}_{0})G(x,t| D){\rm{d}}{D},\end{eqnarray} \tag{ 51 }$$

which may be defined in general as a dynamic superstatistics because of the dependence of p_D(D, t) on t. We obtain an explicit solution by means of the Mellin transform following the same procedure as described in appendix A.2,

$$\begin{eqnarray}{f}_{\mathrm{ggBM}}(x,t| {D}_{0}=0) & = & {({\pi }^{2}{D}_{\star }t(1-\exp (-2t{\sigma }^{2}/{D}_{\star })))}^{-1/2}\\ & & \times {K}_{0}\left(\displaystyle \frac{| x| }{\sqrt{{D}_{\star }t(1-\exp (-2t{\sigma }^{2}/{D}_{\star }))}}\right),\end{eqnarray} \tag{ 52 }$$

where K_ν(z) is the modified Bessel function of second type. The asymptotic behaviour for $| x| \to \infty$ is given by the exponential

$$\begin{eqnarray}&&{f}_{\mathrm{ggBM}}(x,t| {D}_{0}=0)\sim \displaystyle \frac{1}{\sqrt{2\pi | x| \sqrt{{D}_{\star }t(1-{{\rm{e}}}^{-2t{\sigma }^{2}/{D}_{\star }})}}}\exp \left(-\displaystyle \frac{| x| }{\sqrt{{D}_{\star }t(1-{{\rm{e}}}^{-2t{\sigma }^{2}/{D}_{\star }})}}\right).\end{eqnarray} \tag{ 53 }$$

However, in comparison with the result (9) in the equilibrium situation we now observe a different time scaling in the exponent. For ST we see that

$$\begin{eqnarray}&&{f}_{\mathrm{ggBM}}^{{\rm{ST}}}(x,t| {D}_{0}=0)\sim \displaystyle \frac{1}{\sqrt{\sqrt{2}\pi \sigma t| x| }}\exp \left(-\displaystyle \frac{| x| }{\sqrt{2}\sigma t}\right),\end{eqnarray} \tag{ 54 }$$

while at LT

$$\begin{eqnarray}&&{f}_{\mathrm{ggBM}}^{{\rm{LT}}}(x,t| {D}_{0}=0)\sim \displaystyle \frac{1}{\sqrt{2\pi | x| \sqrt{{D}_{\star }t}}}\exp \left(-\displaystyle \frac{| x| }{\sqrt{{D}_{\star }t}}\right).\end{eqnarray} \tag{ 55 }$$

Comparing the ST PDF in (54) with the DD model obtained in (45) we notice that they show a difference in the time scaling of a factor $\sqrt{2}$ which is exactly what we observe in figure 12.

Starting from equation (51) the MSD can be written as

$$\begin{eqnarray}\langle {X}_{\mathrm{ggBM}}^{2}(t)\rangle & = & \langle {(\sqrt{2D(t)}W(t))}^{2}\rangle =2\langle D(t)\rangle t\\ & = & [{D}_{\star }(1-{{\rm{e}}}^{-2t{\sigma }^{2}/{D}_{\star }})+2{D}_{0}{{\rm{e}}}^{-2t{\sigma }^{2}/{D}_{\star }}]t.\end{eqnarray} \tag{ 56 }$$

Note that this result is valid for any initial conditions D₀, not only for the case D₀ = 0. As already suggested above, the scaling of the variance is no longer linear at all times. According to the relation between the parameters it is possible to observe the different scaling behaviours

$$\begin{eqnarray}\langle {X}_{\mathrm{ggBM}}^{2}(t)\rangle =\left\{\begin{array}{ll}2{D}_{0}t, & {\sigma }^{2}t\ll {D}_{0}\\ 2{D}_{0}t+2{\sigma }^{2}{t}^{2}, & {D}_{0}\ll {\sigma }^{2}t\ll {D}_{\star }\\ {D}_{\star }t & {\sigma }^{2}t\gg {D}_{\star }\end{array}\right.,\ 0\leqslant {D}_{0}\ll {D}_{\star },\end{eqnarray} \tag{ 57 }$$

$$\begin{eqnarray}\langle {X}_{\mathrm{ggBM}}^{2}(t)\rangle =\left\{\begin{array}{ll}2{D}_{0}t-4{\sigma }^{2}({D}_{0}/{D}_{\star }){t}^{2} & {\sigma }^{2}t\ll {D}_{\star }\\ {D}_{\star }t & {\sigma }^{2}t\gg {D}_{\star }\end{array}\right.\quad \mathrm{elsewhere}.\end{eqnarray} \tag{ 58 }$$

Thus, when ${D}_{0}\ll {D}_{\star }$ we observe three regimes for the MSD. When D₀ = 0 or when the relation D₀ ≪ D_⋆does not hold we directly observe an initial ballistic behaviour followed by the stationary linear trend. This behaviour is nicely corroborated in figure 13.

A growing range of systems are being revealed which exhibit Brownian yet non-Gaussian diffusion dynamics. Often, an exponential (Laplace) shape of the displacement PDF is observed, however, also stretched Gaussian shapes have been reported. The comparison of diffusion processes recorded by new experimental techniques suggests that the complexity and inhomogeneity of the medium, interpreted as the cause of non-Gaussian behaviour, may influence the spreading of particles in specific fashion and at different levels. In particular, experiments have demonstrated that a non-Gaussian dynamic may persist throughout the observation window and that there are systems that, instead, at LT, exhibit a crossover to Gaussian diffusion. In this article we introduced an analytic approach to generate a random and time-dependent diffusivity with specific features and we proposed two possible models for the spreading dynamics of particles in complex systems: one belonging to the class of ggBM and the other supporting the idea of DD.

We saw that the two models have in common the idea that the non-Gaussianity of the PDF is a direct consequence of an inhomogeneity of the environment, represented by a population of diffusivities. The same PDF for the random diffusivity was introduced for both models. We defined an operative set of dynamic stochastic equations to study random diffusivity effects within the broad class of generalised Gamma distributions. This includes the Gamma distribution, or the exponential PDF which produces the Laplace distribution for the particle displacements.

We observed that the main difference between the ggBM and the DD model is the description of the particle dynamics in the LT regime, corresponding to different physical scenarios for the environment. GgBM does not consider an active dynamics of the environment, and the characteristic that mainly influences the particle motion is the randomness of the medium. This means that the statistical features of the medium completely drive the particles in their entire motion. In contrast the DD model supports the idea of randomly evolving diffusivity corresponding to a dynamics also for the environment. In this way the particles evolve experiencing both a continuous variability in time and a stochasticity in the ensemble. The first model delineates a specific non-Gaussian dynamics for the entire diffusion process, while the second allows to describe a transition from a non-Gaussian to a Gaussian diffusion. In fact, it was shown that the ST non-Gaussian dynamics is the same in the two models, whereas at longer times the ggBM model retains the diffusivity distribution and the DD model leads to an effective value for the diffusivity.

We here also studied the influence of non-equilibrium initial conditions for the diffusivity dynamics and found two main effects. First, the non-equilibrium case breaks the equivalence of the DD and the dynamic generalisation of the ggBM models at ST and, second, it causes changes in the temporal evolution of the MSD. In this case the ggBM model, which in the static case we showed to represent a stochastic interpretation of superstatistical Brownian motion, describes what we may call a dynamical superstatistics that leads to the presence of different time scaling regimes in the process. The DD model, which we investigated in this case via a subordination approach, at ST can no longer be described through a superstatistic approximation, since the subordination results in that regime diverge from the behaviour of ggBM. Furthermore we observed different time scaling regimes for the DD model, as well. Nevertheless, we note that for both models we never obtained an anomalous time scaling for the MSD, only a crossover between ballistic and linear (Brownian or Fickean) behaviour. In the LT regime we obtained a description of the two models which is in agreement with the one for the equilibrium case, as it should be.

It will be interesting to generalise the present findings to anomalous dynamics with stochastic diffusivity by implementing different types of noise. Maintaining the same population of diffusivities the results obtained for the PDF of the particle displacement will not be affected, yet the MSD scaling will become anomalous.

AVC and RM acknowledge funding from the DFG within project ME 1535/6-1. This research is supported by the Basque Government through the BERC 2014–2017 programme and by the Spanish Ministry of Economy and Competitiveness MINECO through BCAM Severo Ochoa excellence accreditation SEV-2013-0323 and through project MTM2016-76016-R MIP.

In this appendix we provide different methods to solve the integral representing the non-Gaussian PDF of the two models discussed in this work,

$$\begin{eqnarray}&&\bar{P}(x,t)={\int }_{0}^{\infty }{p}_{D}(D)G(x,t| D){\rm{d}}{D},\end{eqnarray} \tag{ A.1 }$$

where $G(x,t| D)$ represents a Gaussian distribution and p_D(D) is the generalised Gamma distribution (4).

A.1. Computation via Fox H-function

Recalling equation (A.1) we have

$$\begin{eqnarray}\bar{P}(x,t) & = & {\displaystyle \int }_{0}^{\infty }\displaystyle \frac{\eta }{{D}_{\star }^{\nu }{\rm{\Gamma }}(\nu /\eta )}{D}^{\nu -1}{{\rm{e}}}^{-{(D/{D}_{\star })}^{\eta }}\displaystyle \frac{1}{\sqrt{4\pi {Dt}}}{{\rm{e}}}^{-\tfrac{{x}^{2}}{4{Dt}}}{\rm{d}}{D}\\ & = & \displaystyle \frac{\eta }{{D}_{\star }^{\nu }{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi t}}{\displaystyle \int }_{0}^{\infty }{D}^{\nu -3/2}{{\rm{e}}}^{-{D}^{\eta }}{{\rm{e}}}^{-\lambda {D}_{\star }/D}{\rm{d}}{D},\end{eqnarray} \tag{ A.2 }$$

where we set $\lambda ={x}^{2}/4{D}_{\star }t$. Changing the variable of integration to $y={(D/{D}_{\star })}^{\eta }$ we get

$$\begin{eqnarray}\bar{P}(x,t) & = & \displaystyle \frac{\eta }{{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi {D}_{\star }t}}{\displaystyle \int }_{0}^{\infty }{y}^{\tfrac{\nu -3/2}{\eta }}{{\rm{e}}}^{-y}{{\rm{e}}}^{-\lambda {y}^{-1/\eta }}\displaystyle \frac{1}{\eta }{y}^{\tfrac{1}{\eta }-1}{\rm{d}}{y}\\ & = & \displaystyle \frac{1}{{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi {D}_{\star }t}}{\displaystyle \int }_{0}^{\infty }{y}^{-1-(1/2-\nu )/\eta }{{\rm{e}}}^{-y-\lambda {y}^{-1/\eta }}{\rm{d}}{y}.\end{eqnarray} \tag{ A.3 }$$

With the identification

with the Fox H-function and exploiting some (very convenient) properties of the H-function [101] we then obtain

The Fox function is defined as a generalised Mellin-Barnes integral and has very convenient properties under integral transformations. The Fox function comprises a large range of special functions, including Mejer's G-function, hypergeometric functions, or Bessel functions [102]. In the notation used here the vertical line separates the argument from the function's parameters, and the horizontal line denotes the lack of upper parameters [102].

Recalling that $\lambda ={x}^{2}/4{D}_{\star }t$, we finally obtain

The series expansion of this function reads [102]

The asymptotic behaviour is then obtained in the form [102]

$$\begin{eqnarray}&&{f}_{\mathrm{DD}}^{\mathrm{ST}}(x,t)\sim \displaystyle \frac{1}{{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi {D}_{\star }t}}{\left(\displaystyle \frac{{x}^{2}}{4{D}_{\star }t}\right)}^{\tfrac{2\nu -\eta -1}{2(\eta +1)}}\exp \left(-\displaystyle \frac{\eta +1}{\eta }{\eta }^{\tfrac{1}{\eta +1}}{\left(\displaystyle \frac{{x}^{2}}{4{D}_{\star }t}\right)}^{\eta /(\eta +1)}\right),\end{eqnarray} \tag{ A.8 }$$

for $| x| \to \infty$.

A.2. Computation via Mellin transform

It is possible to rearrange the integral in equation (A.1) as a convolution integral,

$$\begin{eqnarray}\bar{P}(x,t) & = & {\displaystyle \int }_{0}^{\infty }{p}_{D}(D)\displaystyle \frac{{\rm{\exp }}(-\tfrac{{(x/\sqrt{t})}^{2}}{4D})}{\sqrt{4\pi t}\sqrt{D}}{\rm{d}}{D}\\ & = & \displaystyle \frac{2}{\sqrt{4t}}{\displaystyle \int }_{0}^{\infty }\sqrt{D}{p}_{D}({(\sqrt{D})}^{2})\displaystyle \frac{{\rm{\exp }}(-\tfrac{{(x/\sqrt{D})}^{2}}{4})}{\sqrt{\pi }}\displaystyle \frac{{\rm{d}}\sqrt{D}}{\sqrt{D}}\\ & = & \displaystyle \frac{1}{\sqrt{4t}}{\displaystyle \int }_{0}^{\infty }2\xi {p}_{D}({\xi }^{2}){M}_{1/2}\left(\displaystyle \frac{\bar{x}}{\xi }\right)\displaystyle \frac{{\rm{d}}\xi }{\xi },\end{eqnarray} \tag{ A.9 }$$

where we defined $\bar{x}=x/{t}^{1/2}$ and $\xi ={D}^{1/2}$, and M_1/2 denotes the M-Wright function with parameter β = 1/2 [78]. Considering the convolution formula for the Mellin transform

$$\begin{eqnarray}&&{\int }_{0}^{\infty }f(\xi )g\left(\displaystyle \frac{x}{\xi }\right)\displaystyle \frac{{\rm{d}}\xi }{\xi }\ \mathop{\longleftrightarrow }\limits^{{ \mathcal M }}\ {f}^{M}(s){g}^{M}(s),\end{eqnarray} \tag{ A.10 }$$

and remembering the property

$$\begin{eqnarray}&&{x}^{\beta }f({{ax}}^{h})\ \mathop{\longleftrightarrow }\limits^{{ \mathcal M }}\ {h}^{-1}{a}^{-(s+\beta )/h}{f}^{M}\left(\displaystyle \frac{s+\beta }{h}\right),\end{eqnarray} \tag{ A.11 }$$

we compute the Mellin transform of the obtained integral in equation (A.9), recovering

$$\begin{eqnarray}&&\sqrt{4t}\bar{P}(\bar{x})={\displaystyle \int }_{0}^{\infty }2\xi {p}_{D}({\xi }^{2}){M}_{1/2}\left(\displaystyle \frac{\bar{x}}{\xi }\right)\displaystyle \frac{{\rm{d}}\xi }{\xi }\mathop{\longleftrightarrow }\limits^{{ \mathcal M }}\,{p}_{D}^{M}\left(\displaystyle \frac{s+1}{2}\right){M}_{1/2}^{M}(s).\end{eqnarray} \tag{ A.12 }$$

The Mellin transforms for the M-Wright function [78] and the generalised Gamma distribution [102] are known and given by

$$\begin{eqnarray}{M}_{\beta }(x) & & \mathop{\longleftrightarrow }\limits^{{ \mathcal M }}\,\displaystyle \frac{{\rm{\Gamma }}(s)}{{\rm{\Gamma }}(1+(s-1)/\beta )},\\ {\gamma }_{(\nu ,\eta )}^{{\rm{gen}}}(x) & & \mathop{\longleftrightarrow }\limits^{{ \mathcal M }}\,{D}_{\star }^{s-1}\displaystyle \frac{{\rm{\Gamma }}((\nu +s-1)/\eta )}{{\rm{\Gamma }}(\nu /\eta )}.\end{eqnarray} \tag{ A.13 }$$

We can thus rewrite equation (A.12) in the form

$$\begin{eqnarray}&&{p}_{D}^{M}\left(\displaystyle \frac{s+1}{2}\right){M}_{1/2}^{M}(s)={D}_{\star }^{\tfrac{s-1}{2}}\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{\nu +(s+1)/2-1}{\eta }\right)}{{\rm{\Gamma }}\left(\tfrac{\nu }{\eta }\right)}\displaystyle \frac{{\rm{\Gamma }}(s)}{{\rm{\Gamma }}\left(1+\tfrac{1}{2}(s-1)\right)}\end{eqnarray} \tag{ A.14 }$$

$$\begin{eqnarray} & = & \displaystyle \frac{2{D}_{\star }^{\tfrac{s-1}{2}}}{{\rm{\Gamma }}\left(\tfrac{\nu }{\eta }\right)}\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{\nu -1/2+s/2}{\eta }\right){\rm{\Gamma }}(s-1)}{{\rm{\Gamma }}\left(\tfrac{s-1}{2}\right)}\\ & = & \displaystyle \frac{2{D}_{\star }^{\tfrac{s-1}{2}}}{{\rm{\Gamma }}\left(\tfrac{\nu }{\eta }\right)}\displaystyle \frac{{\rm{\Gamma }}\left(\tfrac{\nu -1/2+s/2}{\eta }\right){\rm{\Gamma }}((s-1)/2+1/2)}{{2}^{1-(s-1)}\sqrt{\pi }}\\ & = & \displaystyle \frac{1}{2\sqrt{{D}_{\star }\pi }{\rm{\Gamma }}\left(\tfrac{\nu }{\eta }\right)}{\left(\displaystyle \frac{1}{4{D}_{\star }}\right)}^{-s/2}{\rm{\Gamma }}\left(\displaystyle \frac{\nu -1/2}{\eta }+\displaystyle \frac{s/2}{\eta }\right){\rm{\Gamma }}\left(\displaystyle \frac{s}{2}\right).\end{eqnarray} \tag{ A.15 }$$

Now we notice that the Mellin transform of the H-function is [102]

$$\begin{eqnarray}&&{H}_{p,q}^{m,n}\left[{ax}\left|\begin{array}{l}({a}_{p},{A}_{p})\\ ({b}_{q},{B}_{q})\end{array}\right.\right]\mathop{\,}\limits^{{ \mathcal M }}\longleftrightarrow \,{a}^{-s}\displaystyle \frac{\left(\displaystyle \prod _{j=1}^{m}{\rm{\Gamma }}({b}_{j}+{B}_{j}s)\right)\left(\displaystyle \prod _{j=1}^{n}{\rm{\Gamma }}(1-{a}_{j}-{A}_{j}s)\right)}{\left(\displaystyle \prod _{j=m+1}^{q}{\rm{\Gamma }}(1-{b}_{j}-{B}_{j}s)\right)\left(\displaystyle \prod _{j=n+1}^{p}{\rm{\Gamma }}({a}_{j}+{A}_{j}s)\right)}.\end{eqnarray} \tag{ A.16 }$$

Thus, recalling also the property of the Mellin transform in equation (A.11) we obtain that

and finally

The result here recovered is consistent with equation (A.6).

A.3. Asymptotic trend via Laplace method

Starting again from equation (A.1) it is also possible to calculate directly the asymptotic behaviour through the Laplace method. We introduce the new variable $y={D}_{\star }/D$ in equation (A.1),

$$\begin{eqnarray}\bar{P}(x,t) & = & \displaystyle \frac{\eta }{{D}_{\star }^{\nu }{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi t}}{\displaystyle \int }_{0}^{\infty }{D}^{\nu -3/2}{{\rm{e}}}^{-{(D/{D}_{\star })}^{\eta }}{{\rm{e}}}^{-\lambda {D}_{\star }/D}{\rm{d}}{D}\\ & = & \displaystyle \frac{\eta }{{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi {D}_{\star }t}}{\displaystyle \int }_{0}^{\infty }{y}^{-\nu -1/2}{{\rm{e}}}^{-{y}^{-\eta }-\lambda y}{\rm{d}}{y}\end{eqnarray} \tag{ A.19 }$$

Now the integral looks like a Laplace integral of the form

$$\begin{eqnarray}&&I(\lambda )={\int }_{0}^{\infty }f(y){{\rm{e}}}^{-\lambda y}{\rm{d}}{y}.\end{eqnarray} \tag{ A.20 }$$

In order to apply the Laplace method we need f(0) $\ne$ 0 which is not our case since f(0) = 0 together with all its derivatives. Thus, to evaluate the asymptotics, we define the maximum of the function,

$$\begin{eqnarray}&&\phi (y)=-\lambda y-{y}^{-\eta },\end{eqnarray} \tag{ A.21 }$$

which is located at y_m = (η/λ)^1/η + 1. Introducing of the new variable $\bar{z}=z{\eta }^{\tfrac{1}{\eta +1}}$ the integral becomes

$$\begin{eqnarray}I(\lambda ) & = & {\lambda }^{\tfrac{2\nu -1}{2(\eta +1)}}{\displaystyle \int }_{0}^{\infty }{\bar{z}}^{\tfrac{2\nu -1}{2(\eta +1)}}\exp [-{\lambda }^{\tfrac{\eta }{\eta +1}}({\bar{z}}^{-\eta }+\bar{z})]{\rm{d}}\bar{z}\\ & = & {\lambda }^{\tfrac{2\nu -1}{2(\eta +1)}}{\displaystyle \int }_{0}^{\infty }{\bar{z}}^{\tfrac{2\nu -1}{2(\eta +1)}}{{\rm{e}}}^{\bar{\lambda }f(\bar{z})}{\rm{d}}\bar{z},\end{eqnarray} \tag{ A.22 }$$

where we defined $f(\bar{z})=-\bar{z}-{\bar{z}}^{-\eta }$ and $\bar{\lambda }={\lambda }^{\eta /(\eta +1)}$. Now the standard Laplace method can be applied considering that the function $f(\bar{z})$ reaches its maximum at ${\bar{z}}_{m}={\eta }^{1/(\eta +1)}$, such that

$$\begin{eqnarray}I(\lambda ) & = & {\lambda }^{\tfrac{2\nu -1}{2(\eta +1)}}{\bar{z}}_{m}^{\tfrac{2\nu +1}{2(\eta +1)}}{{\rm{e}}}^{\bar{\lambda }f({\bar{z}}_{m})}\sqrt{\displaystyle \frac{2\pi }{\bar{\lambda }| f^{\prime\prime} ({\bar{z}}_{m})| }}\\ & = & \sqrt{\displaystyle \frac{2\pi }{\eta +1}}{\eta }^{\tfrac{2\nu +\eta +1}{2{(\eta +1)}^{2}}}\,{\lambda }^{\tfrac{2\nu -\eta -1}{2(\eta +1)}}\exp \left[-\displaystyle \frac{\eta +1}{\eta }{\eta }^{\tfrac{1}{\eta +1}}{\lambda }^{\tfrac{\eta }{\eta +1}}\right].\end{eqnarray} \tag{ A.23 }$$

This finally leads to

$$\begin{eqnarray}&&\bar{P}(x,t)\simeq \displaystyle \frac{{\eta }^{\tfrac{2\nu +\eta +1}{2{(\eta +1)}^{2}}+1}}{{\rm{\Gamma }}(\nu /\eta )\sqrt{4\pi {D}_{\star }t}}\sqrt{\displaystyle \frac{2\pi }{\eta +1}}{\left(\displaystyle \frac{{x}^{2}}{4{D}_{\star }t}\right)}^{\tfrac{2\nu -\eta -1}{2(\eta +1)}}\exp \left[-\displaystyle \frac{\eta +1}{\eta }{\eta }^{\tfrac{1}{\eta +1}}{\left(\displaystyle \frac{{x}^{2}}{4{D}_{\star }t}\right)}^{\tfrac{\eta }{\eta +1}}\right],\end{eqnarray} \tag{ A.24 }$$

for $| x| \to \infty$. This result is, up to a numerical prefactor, identical to the asymptotic behaviour obtained in (A.1).