Brownian yet non-Gaussian diffusion: from superstatistics to subordination of diffusing diffusivities

# Brownian yet non-Gaussian diffusion: from superstatistics to subordination of diffusing diffusivities

Aleksei V. Chechkin Institute for Physics & Astronomy, University of Potsdam, 14476 Potsdam-Golm, Germany INFN, Padova Section, and Department of Physics and Astronomy ”G. Galilei”, University of Padova, Via Marzolo 8 35131, Padova, Italy Akhiezer Institute for Theoretical Physics, Kharkov 61108, Ukraine    Flavio Seno INFN, Padova Section, and Department of Physics and Astronomy ”G. Galilei”, University of Padova, Via Marzolo 8 35131, Padova, Italy    Ralf Metzler Institute for Physics & Astronomy, University of Potsdam, 14476 Potsdam-Golm, Germany    Igor M. Sokolov Institute of Physics, Humboldt University Berlin, Newtonstrasse 15, D-12489 Berlin, Germany
July 3, 2019
###### Abstract

A growing number of biological, soft, and active matter systems are observed to exhibit normal diffusive dynamics with a linear growth of the mean squared displacement, yet with a non-Gaussian distribution of increments. Based on the Chubinsky-Slater idea of a diffusing diffusivity we here establish and analyze a minimal model framework of diffusion processes with fluctuating diffusivity. In particular, we demonstrate the equivalence of the diffusing diffusivity process with a superstatistical approach with a distribution of diffusivities, at times shorter than the diffusivity correlation time. At longer times a crossover to a Gaussian distribution with an effective diffusivity emerges. Specifically, we establish a subordination picture of Brownian but non-Gaussian diffusion processes, that can be used for a wide class of diffusivity fluctuation statistics. Our results are shown to be in excellent agreement with simulations and numerical evaluations.

## I Introduction

Thermally driven diffusive motion belongs to the fundamental physical processes. To a big extent inspired by the groundbreaking experiments of Robert Brown in the 1820ies brown () the theoretical foundations of the theory of diffusion were then laid by Einstein, Sutherland, Smoluchowski, and Langevin between 1905 and 1908 einstein (); sutherland (); smoluchowski (); langevin (). On their basis novel experiments, such as the seminal works by Perrin and Nordlund perrin (); nordlund (), in turn delivered ever better quantitative information on molecular diffusion as well as the atomistic nature of matter. Typically, we now identify two fundamental properties with Brownian diffusive processes: (i) the linear growth in time of the mean squared displacement (MSD)

 (1)

typically termed normal (Fickian) diffusion. Here denotes the spatial dimension and is called the diffusion coefficient. (ii) The second property is the Gaussian shape

 P(r,t)=1(4πDt)d/2exp(−r24Dt) (2)

of the probability density function to find the diffusing particle at position at some time vankampen (). From a more mathematical viewpoint the Gaussian emerges as limit distribution of independent, identically distributed random variables (the steps of the random walk) with finite variance and in that sense assumes a universal character gnedenko ().

Deviations from the linear time dependence (1) are routinely observed. Thus, modern microscopic techniques reveal anomalous diffusion with the power-law dependence of the MSD, where according to the value of the anomalous diffusion exponent we distinguish subdiffusion for and superdiffusion with bouchaud (); report (); pccp (); bba (); hoefling (). Examples for subdiffusion of passive molecular and submicron tracers abound in the cytoplasm of living biological cells golding (); goldinga (); goldingb () and in artificially crowded fluids banks (), as well as in quasi two-dimensional systems such as lipid bilayer membranes membranes (); membranesa (); membranesb (); membranesc (). Superdiffusion is typically associated with active processes and also observed in living cells elbaum (). Anomalous diffusion processes arise due to the loss of independence of the random variables, divergence of the variance of the step length or the mean of the step time distribution, as well as due to the tortuosity of the embedding space. The associated probability density function of anomalous diffusion processes may have both Gaussian and non-Gaussian shapes bouchaud (); report (); pccp ().

A new class of diffusive dynamics has recently been reported in a number of soft matter, biological and other complex systems: in these processes the MSD is normal of the form (1), however, the probability density function is non-Gaussian, typically characterized by a distinct exponential shape

 P(r,t)≃exp(−|r|λ(t)), (3)

with the decay length granick (). This form of the probability density function is also sometimes called a Laplace distribution. The Brownian yet non-Gaussian feature appears quite robustly in a large range of systems, including beads diffusing on lipid tubes wang () or in networks wang (); gels (), tracer motion in colloidal, polymeric, or active suspensions susp (), in biological cells bio (), as well as the motion of individuals in heterogeneous populations such as nematodes hapca (). For additional examples see bng (); bnga (); bngb (); bngc () and the references in granick (); gary (); klsebastian ().

How can this combination of normal, Brownian scaling of the mean squared displacement be reconciled with the existence of a non-Gaussian probability density function? One argument not brought forth in the discussion of anomalous diffusion above is the possibility that the random variables making up the observed dynamics are indeed not identically distributed. This fact can be introduced in different ways. First, Granick and co-workers wang () as well as Hapca et al. hapca () employed distributions of the diffusivity of individual tracer particles to explain this remarkable behavior: indeed, averaging the Gaussian probability density function (2) for a single diffusivity over the exponential distribution with the mean diffusivity , the exponential form (3) of the probability density function emerges gary (); hapca (). In fact, this idea of creating an ensemble behavior in terms of distributions of diffusivities of individual tracer particles is analogous to the concept of superstatistical Brownian motion: based on two statistical levels describing, respectively, the fast jiggly dynamics of the Brownian particle and the slow environmental fluctuations with spatially local patches of given diffusivity this concept demonstrates how non-Gaussian probability densities arise physically beck (). In what follows we refer to averaging over a diffusivity distribution as superstatistical approach. An important additional observation from experiments that cannot be explained by the superstatistical approach is that “the distribution function will converge to a Gaussian at times greater than the correlation time of the fluctuations” granick (). This is impressively demonstrated, for instance, in Fig. 1C in wang (). This crossover cannot be explained by the superstatistical approach. At the same time the normal-diffusive behavior is not affected by the crossover between the shapes of the distribution.

Second, Chubinsky and Slater came up with the diffusing diffusivity model, in which the diffusion coefficient of the tracer particle evolves in time like the coordinate of a Brownian particle in a gravitational field gary (). For short times they indeed find an exponential form (3). At long times, they demonstrate from simulations that the probability density function crosses over to a Gaussian shape. Jan and Sebastian formalize the diffusing diffusivity model in an elegant path integral approach, which they explictly solve in two spatial dimensions klsebastian (). Their results are consistent with those of Ref. gary ().

Here we introduce a simple yet powerful minimal model for diffusing diffusivities, based on the concept of subordination. Based on a double Langevin equation approach our model is fully analytical, providing an explicit solution for the probability density function in Fourier space. The inversion is easily feasible numerically, and we demonstrate excellent agreement with simulations of the underlying stochastic equations. Moreover, we provide the analytical expressions for the asymptotic behavior at short and long times, including the crossover to Gaussian statistics, and derive explicit results for the kurtosis of the probability density function. The bivariate Fokker-Planck equation for this process and its connection to the subordination concept are established. Finally, we show that at times shorter than the diffusivity correlation time our analytical results are fully consistent with the superstatistical approach. Our approach has the distinct advantage that it is amenable to a large variety of different fluctuating diffusion scenarios.

In what follows we first formulate the coupled Langevin equations for the diffusing diffusivity model. Section 3 then introduces the subordination concept allowing us to derive the exact form of the subordinator as well as the Fourier image of the probability density function. The Brownian form of the MSD is demonstrated and the short and long time limits derived. Moreover, the connection to the superstatistical approach is made. The kurtosis quantifying the non-Gaussian shape of the probability density function is derived. In section 4 the bivariate Fokker-Planck equation for the joint probability density function is analyzed, before drawing our conclusions in section 5. Several Appendices provide additional details.

## Ii Superstatistical approach to Brownian yet non-Gaussian diffusion

As mentioned above, it was suggested by Granick and coworkers granick () as well as by Hapca et al. hapca () that the Laplace distribution

 P(x,t)=1√4⟨D⟩texp(−|x|(⟨D⟩t)1/2) (4)

with effective diffusivity emerging from a standard Gaussian distribution

 G(x,t|D)=1√4πDtexp(−x24Dt) (5)

with diffusivity , through the averaging procedure

 P(x,t)=∫∞0pD(D)G(x,t|D)dD (6)

over . This approach corresponds to the idea of superstatistics beck (): accordingly the overall distribution function of a system of tracer particles, individually moving in sufficiently large, disjunct patches with local diffusivity , becomes the weighted average, where is the stationary state probability density for the particle diffusivities . While in this Section we restrict the discussion to the one-dimensional case, we will also provide results for higher dimensions below.

Fourier transforming Eq. (6) we obtain

 P(k,t)=∫∞0pD(D)e−Dk2tdD=~pD(s=k2t), (7)

where we used the fact that . On the right hand side we identified the integral of over as the Laplace transform to be taken at . Concurrently, the Fourier transform of expression (4) is

 P(k,t)=11+⟨D⟩k2t. (8)

Combining these results and recalling the Laplace transform , we uniquely find that indeed

 pD(D)=1⟨D⟩exp(−D⟨D⟩). (9)

To obtain the Laplace distribution (4) as superstatistical average of elementary Gaussians (5), the necessary distribution of the diffusivities is the exponential (9). This is exactly the result of Granick and coworkers granick () and Hapca et al. hapca (). (We note that Hapca and coworkers also report results for the case of a gamma distribution .)

Now, let us take the Fourier inversion of Eq. (7) and invoke the substitution ,

 P(x,t) = 12π∫∞−∞e−ikx~pD(k2t)dk (10) = 12πt1/2∫∞−∞e−iκx/t1/2~pD(κ2)dκ.

The right hand side defines a scaling function of the form

 P(x,t)=1t1/2F(ζ), (11)

where . Thus the form as function of the similarity variable is an invariant. In particular, no transition of from a Laplace distribution to a different shape is possible in this superstatistical framework. To account for the experimental observation, however, we are seeking a model to explain the crossover from an initial Laplace distribution to a Gaussian shape at long(er) times.

### Anomalous diffusion with exponential, stretched Gaussian, and power law shapes of the probability density function

We briefly digress to mention that for the case of anomalous diffusion with a mean squared displacement of the form a similar phenomena was observed. Namely, for the motion of particles in a viscoelastic environment with a fixed generalized diffusivity of dimension the motion is characterized by the Gaussian goychuk (); pccp ()

 (12)

with scaling variable. Instead, in a recent experimental study observing the motion of labeled messenger RNA molecules in living E.coli and S.cerevisiae cells an exponential distribution of the diffusivity was found spako (),

 pD(Dα)=1D⋆αexp(−DαD⋆α) (13)

on the single trajectory level, pointing at a higher inhomogeneity of the motion than previously assumed. The distribution (13) combined with the Gaussian (12) gives rise to the Laplace distribution spako ()

 Pα(x,t)=1√4D⋆αtαexp(−|x|√D⋆αtα). (14)

Similarly one can show that the stretched Gaussian observed for the lipid motion in protein-crowded lipid bilayer membranes membranesa () emerges from the Gaussian (12) in terms of a modified diffusivity distribution of the form

 pD(Dα)=1Γ(1+1/κ)D⋆αexp(−[DαD⋆α]κ). (15)

In that case the resulting distribution assumes the form

 Pα(x,t)≃exp⎛⎝−c[|x|(4D⋆αtα)1/2]2κ/(1+κ)⎞⎠ (16)

with an additional power law term in , see Appendix A. Depending on the value of one can then obtain stretched Gaussian shapes for or even broader than exponential forms (superstretched Gaussians).

We finally note that for a power law distribution

 pD(D)≃D−1−α (17)

with the resulting superstatistical distribution acquires long tails of the form

 Pα(x,t)≃1|x|2α+1, (18)

as demonstrated in Appendix A. This brief discussion shows the need for a more general model for the diffusing diffusivity, the basis of which is established here.

## Iii Langevin model for diffusing diffusivities

To describe Brownian but non-Gaussian diffusion we start with the combined set of stochastic equations

 ddtr(t) = √2D(t)ξ(t), (19a) D(t) = Y2(t), (19b) ddtY(t) = −1τY+ση(t). (19c) The independent noise terms ξ(t) and η(t) are white and Gaussian, and both are specified by their first two moments ⟨ξ(t)⟩=0,⟨ξi(t1)ξj(t2)⟩=δijδ(t1−t2) (19d) ⟨η(t)⟩=0,⟨ηl(t1)ηm(t2)⟩=δlmδ(t1−t2), (19e) for i,j=x,y,z and l,m=1,…,n. As explained below, the dimension n of the process Y(t) may differ from the value of d of the process r(t) in real space. In the above set of coupled stochastic equations (19), expression (19a) designates the well-known overdamped Langevin equation driven by the white Gaussian noise ξ(t) risken (). However, we consider the diffusion coefficient D(t) to be a random function of time, and we express it in terms of the square of the Ornstein-Uhlenbeck process Y(t) (see below for the reasoning). The physical dimension of the latter is [Y]=cm/sec1/2. In Eq. (19c) the correlation time of the Ornstein-Uhlenbeck process is τ, and σ of units [σ]=cm/sec characterizes the amplitude of the fluctuations of Y. We complete the set of stochastic equations with the initial conditions, chosen as r(0)=0,Y(0)=Y0. (19f)

Physically, the choice of the above set of dynamic equations corresponds to the following reasonings. In the diffusing diffusivity picture we model the particle motion, on the single trajectory level, by the random diffusivity . Taking as the square of the auxiliary variable guarantees the non-negativity of . This way we avoid the need to impose reflecting boundary condition on at , which is more difficult to handle analytically gary (). The reason to choose the Ornstein-Uhlenbeck process (19c) for is two-fold. First, it makes sure that the diffusivity dynamics is stationary, with a given correlation time. Second, the ensuing distribution has exponential tails, thus guaranteeing the emergence of the Laplace-like distribution for at short times, as we will show. At long times, the above choice corresponds to a particle moving with an effective diffusivity , and thus leads to the crossover to the long time Gaussian behavior of . The above set of Langevin equations not only fulfill these requirements but also allows for an analytical solution, as shown below.

For simplicity, we introduce dimensionless units via the transformations and (and similar for and ). The process is renormalized according to . As detailed in Appendix B we then obtain the set of stochastic equations

 ddtr(t) = √2D(t)ξ(t), (20a) D(t) = Y2(t), (20b) ddtY(t) = −Y+η(t). (20c)

for our minimal diffusing diffusivity model.

We note that the above minimal model for the diffusing diffusivity allows different choices for the number of components of . The number is thus essentially a free parameter of the model. It defines the number of ‘modes’ necessary to describe the random process . This is actually another advantage of the present approach, since it provides additional flexibility.

In the Discussion section we will show that the above compound process is analogous to the Heston model heston () and thus a special case of the Cox-Ingersoll-Ross (CIR) model cir (), which are widely used for return dynamics in financial mathematics. Our approach therefore has a wider appeal beyond stochastic particle dynamics.

### Properties of the Ornstein-Uhlenbeck process

The stochastic equation (20c) contains a linear restoring term, corresponding to the motion of the process in a centered harmonic potential. The formal solution of this Ornstein-Uhlenbeck process reads

 Y(t)=Y0e−t+∫t0η(t′)e−(t−t′)dt′. (21)

The associated autocorrelation function is

 ⟨Y(t1)Y(t2)⟩ = Y20e−(t1+t2)+e−(t1+t2) (22) ×∫t10dt′1∫t20dt′2⟨η(t′1)η(t′2)⟩et′1+t′2 =Y20e−(t1+t2)+n2(e−|t2−t1|−e−(t1+t2)).

Thus, for long times (), we find the exponential decay

 ⟨Y(t1)Y(t2)⟩∼n2e−|t2−t1| (23)

of the autocorrelation, and thus the stationary variance

 ⟨Y2(t)⟩=⟨D⟩st=n2. (24)

We note that the Fokker-Planck equation for this Ornstein-Uhlenbeck process reads

 ∂∂tf(Y,t)=∂∂Y(Yf(Y,t))+12∂2∂Y2f(Y,t). (25)

The distribution converges to the normalized equilibrium Boltzmann form

 fst(Y)=1πn/2e−Y2. (26)

In what follows and in our simulations we assume that the initial condition is taken randomly from the equilibrium distribution (26). Then, the process becomes stationary starting from , and Eq. (23) is exact at all times.

The stationary diffusivity distribution encoded in Eq. (26) in terms of the variable can then be obtained as follows.

(i) In dimension , the variance of in the stationary state is and the mapping to reads

 pstD(D)=∫∞−∞fst(Y)δ(D−Y2)dY=1√πDe−D. (27)

In dimensional form we have

 pstD(D)=1√πD⋆De−D/D⋆, (28)

with

 D⋆=x20t0=σ2τ. (29)

From comparison with the direct superstatistical approach in Section II we see that the pure exponential form (9) in our diffusing diffusivity model is being modified by the additional prefactor . From numerical comparison, however, the exponential dependence is dominating and thus the result (28) practically indistinguishable from (9) for sufficiently large values.

(ii) For , the stationary state variance of is , the mapping from to reads

 pstD(D)=2π∫∞0Yfst(Y)δ(D−Y2)dY=e−D, (30)

where and . Moreover, we made use of the property

 δ(Y2−D)=12√D[δ(Y+√D)+δ(Y−√D)]. (31)

of the -function. In dimensional units, we have

 pstD(D)=1D⋆e−D/D⋆, (32)

in conjunction with relation (29).

(iii) Finally, for we have and

 pstD(D)=4π∫∞0Y2fst(Y)δ(D−Y2)=2√D√πe−D, (33)

where . In dimensional form,

 pstD(D)=2√D√πD3⋆e−D/D⋆. (34)

## Iv Subordination concept for diffusing diffusivities

Subordination, introduced by Bochner bochner (), is an important concept in probability theory feller (). Simply put, a subordinator associates a random time increment with the number of steps of the subordinated stochastic process. For instance, continuous time random walks with power-law distributions of waiting times can be described as a Brownian motion in terms of the number of steps of the process, while the random waiting times are introduced in terms of a Lévy stable subordinator, as originally formulated by Fogedby fogedby () and developed as a stochastic representation of the fractional Fokker-Planck equation report () and generalized master equations for continuous time random walk models igor (); friedrich (); friedricha (); friedrichb ().

Here we apply and extend the subordination concept to a new class of random diffusivity based stochastic processes. Our results for our minimal model of diffusive diffusivities demonstrates that the subordination approach leads to a superstatistical solution at times shorter than typical diffusivity correlation times.

To start with, we note that the stochastic probability density function fulfills the diffusion equation

 ∂∂t¯¯¯¯P(r,t)=D(t)∇2¯¯¯¯P(r,t). (35)

With this in mind we can rewrite the Langevin equation (20a) in the subordinated form

 ddτr(τ) = √2ξ(τ) (36a) ddtτ(t) = D(t). (36b)

After this change of variables the Green function of the diffusion equation has the form

 G(r,τ)=1√(4πτ)dexp(−r24τ) (37)

with . The path variable , for any given instant of time , according to Eq. (36b) is a random quantity. In order to calculate the probability density of the variable at time we need to eliminate the path variable . This is achieved by averaging the Green function (37) over the distribution of in the form

 P(r,t)=∫∞0Tn(τ,t)G(r,τ)dτ. (38)

Here is the probability density function of the process

 τ(t)=∫t0D(t′)dt′=∫t0Y2(t′)dt′. (39)

Equation (38) is but the well known subordination formula, implying the following: the probability for the walker to arrive at position at time equals the probability of being at on the path at time , multiplied by the probability of being at position for this path length , summed over all path lengths fogedby ().

By help of relation (38) we write the Fourier transform

 ^P(k,t)=∫∞−∞eik⋅rP(r,t)dr (40)

in the subordinated form

 ^P(k,t) = ∫∞0Tn(τ,t)^G(k,τ)dτ (41) = ∫∞0Tn(τ,t)e−k2τdτ=~Tn(k2,t),

with . Thus, the Fourier transform of is expressed in terms of the Laplace transform of the density function with respect to ,

 ~Tn(s,t)=∫∞0e−sτTn(τ,t)dτ, (42)

with argument .

The subordination approach established here introduces a superior flexibility into the diffusing diffusivity model. By specific choice of the subordinator density we may study a broad class of normal and anomalous diffusion processes caused by diffusivities, that are randomly varying in time and/or space. In turn, the advantage of our minimal model for diffusing diffusivities introduced here is, that the process is the integrated square of the Ornstein-Uhlenbeck process, for which in the one-dimensional case the Laplace transform of the probability density function is known dankel (),

 ~T1(s,t) = exp(t/2)/[12(√1+2s+1√1+2s) (43) ×sinh(t√1+2s)+cosh(t√1+2s)]1/2.

We thus directly obtain the exact analytical result for the Fourier transform

 ^P(k,t) = exp(t/2)/[12(√1+2k2+1√1+2k2) (44) ×sinh(t√1+2k2)+cosh(t√1+2k2)]1/2.

of the probability density function . The inverse Fourier transform can be performed numerically. Figure 1 demonstrates excellent agreement between this result and simulations of the stochastic starting equations (20a) to (20c). Below we provide analytical estimates of for short and long times and establish a connection of the subordination approach with the superstatistical framework.

Our approach can be easily generalized to the case of the -dimensional Ornstein-Uhlenbeck process considered by Jain and Sebastian klsebastian (). Namely, let us consider

 D(t)=Y2(t), (45)

where is a -dimensional Ornstein-Uhlenbeck process. Since the components of are independent and

 τ(t) = ∫t0Y2(t)dt′ (46) = ∫t0(Y21(t′)+Y22(t′)+…+Y2n(t′))dt′,

the Laplace transform and the characteristic function are simply -fold products of identical, one-dimensional functions (43) and (44), respectively:

 ~Tn(s,t) = exp(nt/2)/[12(√1+2s+1√1+2s) (47) ×sinh(t√1+2s)+cosh(t√1+2s)]n/2.

We thus directly obtain the exact analytical result for the Fourier transform

 ^P(k,t) = exp(nt/2)/[12(√1+2k2+1√1+2k2) (48) ×sinh(t√1+2k2)+cosh(t√1+2k2)]n/2.

This result is consistent with that of Eq. (25) in Ref. klsebastian (), up to numerical factors, which appear due to the difference in numerical coefficients entering the Ornstein-Uhlenbeck process.

### iv.1 Brownian mean squared displacement and leptokurtic behavior

The mean squared displacement and the fourth moment encoded in our minimal model can be directly obtained from the Fourier transform (44) through differentiation,

 ⟨r2(t)⟩ = −∇2k^P(k,t)∣∣k=0, ⟨r4(t)⟩ = ∇4k^P(k,t)∣∣k=0. (49)

For the isotropic case considered here the Laplace operator is defined as

 ∇2k=1kd−1∂∂k(kd−1∂∂k). (50)

Expanding relation (41) for small , we obtain up to the fourth order

 ^P(k,t) = ∫∞0e−k2τTn(τ,t)dτ (51) = 1−k2∫∞0τTn(τ,t)dτ +k42∫∞0τ2Tn(τ,t)dτ+…

From this we directly obtain the mean squared displacement

 ⟨r2(t)⟩ = 2d∫∞0τTn(τ,t)dτ=2d⟨τ⟩ (52) = −2d∂~Tn(s,t)∂s∣∣ ∣∣s=0

and the fourth order moment

 ⟨r4(t)⟩ = 4d(2+d)∫∞0τ2Tn(τ,t)dτ=4d(2+d)⟨τ2⟩ (53) = 4d(d+2)∂2~Tn(s,t)∂s2∣∣ ∣∣s=0

With the results of Appendix D we find

 ⟨r2(t)⟩=dnt=2d⟨D⟩stt, (54)

where is given by Eq. (24), and

 ⟨r4(t)⟩=4d(2+d)⟨D⟩st[−1−e−2t2+t+⟨D⟩stt2]. (55)

Eq. (54) is the famed result of the normal Brownian, linear dispersion of the mean squared displacement with time. In dimensional units the result (54) reads

 (56)

Fig. 2 demonstrates excellent agreement of the analytical result (54) with direct simulations of the set of Langevin equations (20a) to (20c) with respect to both slope and amplitude.

The deviation of the shape of a distribution function from a Gaussian can be conveniently quantified in terms of the kurtosis

 K=⟨r4(t)⟩⟨r2(t)⟩2. (57)

We note that the kurtosis is closely related to the (first) non-Gaussian parameter, introduced in the classical text by Rahman rahman (). Inserting results (54) and (55) we obtain in the short time limit that

 K∼(1+2d)(1+1⟨D⟩st)=⎧⎪⎨⎪⎩9,d=14,d=225/9,d=3 (58)

for the choice . At long times,

 K∼(1+2d)=⎧⎪⎨⎪⎩3,d=12,d=25/3,d=3. (59)

The first relation characterizes exponential distributions according to Equations (66), (70), and (72) derived below, whereas the second result coincides exactly with the kurtosis of the multidimensional Gaussian distribution (note that this kurtosis does not depend on ). Taking along the next higher term in the long time expansion, we observe that the leptokurtosis vanishes as in the long time limit,

 K∼1+2d+2+dd⟨D⟩stt. (60)

In Fig. 3 we compare the analytical result (57) for the kurtosis for based on Equations (54) and (55) with simulations, showing excellent agreement from the short time behavior all the way to the saturation plateau at the Gaussian value . We note that the crossover time from strongly leptokurtic to Gaussian behavior occurs at , which in dimensional units corresponds to the correlation time of the diffusing diffusivity process. Thus in experiments the behavior of the kurtosis as function of time provides a direct means to extract the correlation time of , which also corresponds to the crossover time from the exponential to the Gaussian behavior of the probability density , as will be shown below. We also note that when the process is very highly dimensional and thus large, the exponential tails of do not exist, see Appendix C.

We now derive explicit analytical results for the probability density function in the short and long time limits, starting with the short time limit and its relation to the superstatistical formulation of the diffusing diffusivity. In our further exemplary calculations we take . However, in Appendix C we discuss the situations when and , as well as when goes to infinity while stays finite.

### iv.2 Short time limit

First we concentrate on the shape of the density in the short time limit . In dimensionless units this means that we consider the asymptotic behavior of the characteristic function (48) under the condition , for which

 sinh(t√1+2k2)∼t√1+2k2,cosh(t√1+2k2)∼1, (61)

Together with expression (48) we thus find

 ^P(k,t)∼(1+t)n/2(1+[1+k2]t)n/2∼t−n/2(k2+1t)−n/2. (62)

This expression is indeed normalized, . We can thus perform the inverse Fourier transform to obtain in the short time limit.

(i) For one dimension we find

 P(x,t)∼1πt1/2∫∞0cos(kx)(k2+1/t)1/2dk=1πt1/2K0(xt1/2), (63)

in terms of the Bessel function gradshteyn ()

 K0(aβ)=∫∞0cos(ax)√x2+β2dx. (64)

Apart from the normalization we observe that from Eq. (63) we also derive the Brownian behavior , in accordance with Eqs. (24) and (54).

Keeping in mind that here we are pursuing the large value limit of the scaling variable , we expand the Bessel function in the form abramowitz ()

 K0(z)∼√π2ze−z. (65)

We thus find the asymptotic result

 P(x,t)∼1√2π|x|t1/2exp(−|x|t1/2). (66)

This expression reproduces the exponential shape of the probability density function of the diffusing diffusivity model, with the power-law correction .

Figure 4 demonstrates excellent agreement of our short time result (63) and simulations. For the longest simulated time the wings of the distribution start to show some deviations, indicating that in this case the short time limit is no longer fully justified.

(ii) For we obtain

 P(r,t) = ∫eik⋅r^P(k,t)dk(2π)2 (67) = 12π∫∞0kJ0(kr)^P(k,t)dk,

and thus

 P(r,t)=12πtK0(r√t). (68)

Here we used the relation

 ∫π0cos(krcosφ)dφ=πJ0(kr) (69)

in terms of the modified Bessel function . The distribution is normalized and encodes the Brownian behavior . Expanding the Bessel function as above we find that

 P(r,t)∼12√2πrt3/2e−r/√t. (70)

(iii) Finally, in the asymptotic probability density becomes

 P(r,t) = ∫eik⋅r^P(k,t)dk(2π)3 (71) = 18π3∫∞0k2dk∫π0sinθdθ∫2π0dφeikrcosθ ×t−3/2(k2+1t)−3/2 = 14π2t3/2∫∞0k2dk∫1−1dξeikrξ(k2+1t)−3/2 = 12π2t3/2K0(r√t),

and we have . Expansion of the Bessel function produces the asymptotic exponential behavior

 P(r,t)∼1(2π)3/2r1/2t5/4e−r/√t. (72)

As we will show in the following Subsection, in all dimensions the short time limit reproduces the superstatistical behavior. The subordination formulation, the explicit result (63) in terms of the Bessel function, and the asymptotic exponential (Laplace) form (66) (and their multidimensional analogs (68) and (70) to (72)) constitute our first main result. We note that the asymptotic forms (66), (70), and (72) of the density function by itself leads to the Brownian scaling (54) of the corresponding mean squared displacement. When the exponential shape of the short time behavior is conserved while the subdominant prefactors change, as shown for the case and in Appendix C.

### iv.3 Relation to the superstatistical approximation

Above we formulated the concept of diffusing diffusivities in terms of coupled stochastic equations for the particle position and the random diffusivity . An alternative approach suggested in Refs. granick (); gary () is that of the superstatistical distribution of the diffusivity, as laid out in Section II. In this superstatistical sense the overall distribution function is given as the weighted average of a single Gaussian over the stationary diffusivity distribution,

 Ps(r,t)=∫∞0pstD(D)G(r,t|D)dD. (73)

(i) In dimension our minimal model produces with Eq. (28)

 Ps(x,t) = 12π√D⋆t∫∞01Dexp(−DD⋆−x24Dt)dD (74) = 1π√D⋆tK0(|x|√D⋆t).

(ii) In , we have with Eq. (32)

 Ps(r,t) = 14πD⋆t∫∞01Dexp(−DD⋆−r24Dt)dD (75) = 12πD⋆tK0(r√D⋆t).

(iii) Finally, in , we find with Eq. (34)

 Ps(r,t) = 2π2(4D⋆t)3/2∫∞01Dexp(−DD⋆−r24Dt)dD (76) = 12π2(D⋆t)3/2K0(r√D⋆t).

For all the mean squared displacement acquires the linear Brownian scaling in time , as it should.

This is but exactly the result of our subordination scheme in the short time limit, expressions (66), (68), and (71), written in dimensional form. Thus, in our approach to the diffusing diffusivity the short time regime of the subordination formalism leads directly to the superstatistical result. The reason is as follows: At times less than the diffusivity correlation time the diffusion coefficient does not change considerably, and the subordination scheme describes an ensemble of particles, each diffusing with its own diffusion coefficient. This mimics a spatially inhomogeneous situation, when the local diffusion coefficient is random, but stays constant within confined spatial domains. In this case the ensemble of particles moving in different domains exhibits a superstatistical behavior, as assumed in the original works beck (). However, in any system with finite patch sizes, we would not expect the particles to stay in their local patch of diffusivity forever, thus violating the assumption of the superstatistical approach. Our annealed approach in some sense delivers a mean field approximation to the spatially disordered situation, and adequately describes the transition from short time superstatistical behavior to the Gaussian probability law at long times, which will be shown in the subsequent section. The full consistency in the short time limit between the subordination approach and superstatistics is our second main result.

### iv.4 Long time limit

We now turn to the long time limit encoded in the Fourier transform (44) of the probability density , that is, the times larger than the diffusivity correlation time . In dimensionless units it corresponds to , and the hyperbolic functions assume the limiting behaviors

 sinh(t√1+2k2) ∼ cosh(t√1+2k2) (77) ∼ 12exp(t√1+2k2).

Combined with result (48) we find

 ^P(k,t)∼2n/2exp(nt2[1−√1+2k2])(1+12[√1+2k2+1√1+2k2])n/2. (78)

As in the short time limit above, this expression is normalized, .

Now let us focus on the tails of the probability density , corresponding to the limit , for which Eq. (78) gives , and thus

 P(r,t)∼1(4π⟨D⟩stt)n/2exp(−r24⟨D⟩stt). (79)

At long times the probability density function assumes a Gaussian form, with the effective diffusivity . This is a consequent result given the Ornstein-Uhlenbeck variation of the diffusivity encoded in the starting equations (19b) and (19c): at sufficiently long times the process samples the full diffusivity space and behaves like an effective Gaussian process with renormalized diffusivity. The explicit derivation of the crossover to the Gaussian behavior is our third main result.

Fig. 5 shows the crossover from the initial exponential to the long time Gaussian behavior of the probability density function for by comparison to the Gaussian distribution (79) for short time , the crossover time and the longer time .

## V Bivariate Fokker-Planck equation and relation to the subordination approach

In this section we derive the Fokker-Planck equation corresponding to the set of stochastic equations (20a) to (20c) of our diffusing diffusivity model. We will also establish the relation to the subordination approach of section IV. Note that we here restrict the discussion to the case , as higher dimensional cases are completely equivalent.

Following our notation we thus seek the Fokker-Planck equation for the bivariate probability density function , which has the structure

 ∂∂tf(x,y,t)=Lyf(x,y,t)+y2∂2∂x2f(x,y,t). (80)