Superstatistics and non-Gaussian diffusion

Brownian motion and viscoelastic anomalous diffusion in homogeneous environments are intrinsically Gaussian processes. In a growing number of systems, however, non-Gaussian displacement distributions of these processes are being reported. The physical cause of the non-Gaussianity is typically seen in different forms of disorder. These include, for instance, imperfect “ensembles” of tracer particles, the presence of local variations of the tracer mobility in heteroegenous environments, or cases in which the speed or persistence of moving nematodes or cells are distributed. From a theoretical point of view stochastic descriptions based on distributed (“superstatistical”) transport coefficients as well as time-dependent generalisations based on stochastic transport parameters with built-in finite correlation time are invoked. After a brief review of the history of Brownian motion and the famed Gaussian displacement distribution, we here provide a brief introduction to the phenomenon of non-Gaussianity and the stochastic modelling in terms of superstatistical and diffusing-diffusivity approaches.


Introduction
In the anni mirabiles from 1905 to 1908 Albert Einstein, William Sutherland, Marian Smoluchowski, and Paul Langevin introduced their theories of Brownian motion [1][2][3][4] and Jean Perrin published his seminal single particle tracking experiments [5,6]. Concurrently, Karl Pearson introduced the concept of the random walk [7,8]. Based on this early work, we now understand Brownian motion as the continuum limit of a Pearson walk, in which individual steps are independent and identically distributed. As long as the lengths of individual steps of a tracer particle undergoing such a process have a finite variance, following the central limit theorem the probability density function (PDF) is Gaussian a forteriori, a e-mail: rmetzler@uni-potsdam.de

712
The European Physical Journal Special Topics Fig. 1. Three microscopic trajectories measured by Perrin [5]. The turning points correspond to the particle position measured in 30 s time intervals. and the corresponding mean squared displacement (MSD) has the linear time dependence with diffusivity D 1 in a d-dimensional system. Properties (1) and (2) are the hallmarks of Brownian motion [9,10]. Figures 1 and 2 show exemplary results of the experiments by Perrin. The "random walk" traces shown in Figure 1 show the test particle positions as taken in 30 sec intervals. These positions are connected by straight lines. 1 The trajectories contain relatively few points, as the particles left the focus of the microscope quite quickly. For a quantitative analysis, Perrin therefore used the independence of individual steps and shifted each displacement vector to a common origin, Figure 2. The resulting radial histogram was then used to determine the parameters of the Gaussian distribution from which Perrin ultimately obtained Avogadro's number N A , making use of the famed Einstein-Sutherland-Smoluchowski relation with the mass m of the test particle, the viscosity η, as well as thermal energy k B T and the gas constant R [1][2][3][4]9]. The Gaussian distribution was directly mapped out by Eugen Kappler in his torsional Brownian motion setup using a small mirror suspended on a long, thin quartz thread [11]. Given that the elongations induced by bombarding ambient air  Figure 1 are shifted to a common origin and mapped onto a Gaussian profile [5].
molecules are fairly small, the restoring force by the thread can be assumed to be Hookean. The resulting diffusion in an harmonic confining potential is actually an Ornstein-Uhlenbeck process, whose probability density function remains Gaussian at all times [9,12,13]. Kappler's result displayed in Figure 3 impressively corroborates the expected Gaussian shape.
Normal Brownian diffusion is not the only process with a native Gaussian PDF. We mention two classes of stochastic processes, both describing anomalous diffusion of the power-law form [14][15][16] Here, the anomalous diffusion coefficient D α has dimension cm 2 /sec α , and, depending on the value of the anomalous diffusion exponent α, we distinguish subdiffusion (0 < α < 1) and superdiffusion (α > 1).
A closely related process is described by the fractional Langevin equation (FLE) [37][38][39][40] with 0 < α < 1 and γ = γ 0 t α−2 . After an initial ballistic regime x 2 (t) t 2 the motion crosses over to the subdiffusive MSD (4) at long times. In contrast to FBM, the fractional Langevin equation (6) fulfils detailed balance and thus in a confining external potential relaxes to thermal equilibrium. This requires that the noise covariance function is coupled to the power-law friction kernel, ξ [37,38,41]. Subdiffusion here emerges despite the persistent (positively correlated) noise covariance, as large noise spikes are compensated by large friction values, which is the physical explanation that equation (6) describes viscoelastic systems measured in [42][43][44][45][46][47]49,50]. We note that within the framework of FBM or FLE motion crossovers from anomalous (for FLE motion at time scales beyond the initial ballistic motion) to normal diffusion or between anomalous diffusion characterised by different scaling exponents can be included by tempering of the fractional Gaussian noise [51].
Here we address cases in which the MSD and displacement correlation function point at either normal (Fickian) or anomalous diffusion driven by fractional Gaussian noise, while the displacement PDF is clearly non-Gaussian. These "Brownian yet non-Gaussian" and "viscoelastic yet non-Gaussian" processes have been reported in numerous systems. We provide a physical scenario for these cases and present different theoretical models describing these phenomena, including the superstatistical formulation and the diffusing-diffusivity model.

Brownian yet non-Gaussian diffusion and superstatistics
The combination of the linear time dependence (2) of the MSD with a non-Gaussian PDF is in a way counterintuitive, as it contradicts our naive expectation that such "normal" diffusion should be Gaussian. The case for this "Brownian [or Fickian] yet non-Gaussian" diffusion phenomenon was championed by Granick [52], whose group reported non-Gaussian diffusion for the Fickian motion of submicron tracers along linear tubes and in entangled actin networks [53], as well as for tracer dynamics in hard sphere colloidal suspensions [54]. Other experimental evidence for non-Gaussian behaviours comes from the diffusion of nanoparticles in nanopost arrays [55], colloidal nanoparticles adsorbed at fluid interfaces [56][57][58] and moving along membranes and inside colloidal suspension [59], and the motion of nematodes [60]. We also mention the non-Gaussian dynamics in disordered solids such as glasses and supercooled liquids [61][62][63] as well as interfacial dynamics [64][65][66] and dynamics in actively remodelling semiflexible networks [67,68]. Figure 4 shows, along with a cartoon of the tracer bead in the F-actin network, the original data from [54]. We see both the linear time dependence of the MSD and the non-Gaussian shape of the displacement PDF: while at short distances the shape is Gaussian the tails of the PDF are exponential ("Laplace distribution") [54], Supplementing this information, Figure 5 demonstrates that the exponential tails are present at different lag times and collapse to a master curve with exponential tail. Concurrently, the width λ is shown to scale with time as λ(t) t 1/2 , such that the position-time scaling is diffusive in the sense that r 2 t [54].
A way to understand this non-Gaussianity was already proposed in the paper by Granick [52], namely, the concept of superstatistics as formulated by Beck and Cohen. Accordingly, the measured PDF P sup (r, t) of an ensemble of particles corresponds to the mean [69][70][71]

716
The European Physical Journal Special Topics showing a λ(t) √ t scaling, for radii and mesh sizes a = 50 nm and ξ = 300 nm (crosses), a = 100 nm and ξ = 450 nm (triangles), and a = 100 nm and ξ = 300 nm (circles). Reproduced from [52]. (This figure is subject to copyright protection and is not covered by a Creative Commons license).
where P (r, t|D 1 ) is the Gaussian (1) for a specific value D 1 of the diffusion coefficient, and p(D 1 ) is the PDF of D 1 values. The MSD of this process becomes Independently of the diffusivity distribution p(D 1 ), that is, the particle ensemble is characterised by Fickian diffusion, with an effective diffusion coefficient D eff = D 1 . Within the superstatistical approach it is straightforward to show that the Laplace distribution (7) uniquely emerges from an exponential diffusivity distribution, [72]. The case of a gamma distribution for p(D 1 ) was considered in [60], see also [73]. Finally, it was shown that power-law forms for p(D 1 ) effect superstatistical distributions P sup (r, t) with power-law tails [72,74].
Physically, superstatistics naturally emerges when we consider an imperfect "ensemble" of diffusing particles with a distribution of mobilities. This was in fact the case in Perrin's original measurements and even occurs in contemporary experiments with tracer beads that can be ordered from specialist providers: when we measure ensembles of such particles, the formulations (8) and (9) naturally emerge. In particular, in this case the value D 1 -specific for each particle -is constant in time. The picture envisaged by Beck and coworkers was in fact that of an heterogeneous environment. Imagine that all particles are identical, but each particle is moving on its own patch in space, characterised by a specific D 1 value. Again, measuring over an ensemble of particles in an array of patches with different D 1 produces the behaviour encoded in (8) and (9), and p(D 1 ) is then given by the statistic of the patches. In this scenario, of course, once a given particle reaches the border of its patch, it will move into a patch with a different D 1 value. Once all particles explore many different local environments, the overall ensemble behaviour will cross over to an effectively Gaussian statistic: beyond some correlation time the system is Brownian and Gaussian, with an effective diffusion constant. A specific scenario for such a crossover dynamic is discussed in more detail in Section 3.
Beck and coworkers applied superstatistical concepts to a range of dynamic phenomena, including turbulence [75,76], high energy physics [77], power grid fluctuations [78], and delay time statistics of British trains [79]. Naturally, concepts similar to superstatistics were previously discussed, for instance, by Shraiman and Siggia [80] in the context of stretched exponential distributions in turbulence [81]. However, Beck introduced superstatistics as a physical concept and made the connection to statistical mechanics [82,83]. We note that the superstatistical formulation was also achieved starting from a stochastic Langevin equation [84]. We also note that a similar, random-parameter formulation of diffusion processes is given by the concept of (generalised) grey Brownian motion [73,[85][86][87][88].
That the superstatistical ensemble is characterised by the effective MSD (9) is in fact not surprising, as such a behaviour necessarily follows for any shape P (r, t) = t −d/2 g(r 2 /t) of the PDF in terms of some scaling function g(·). To show this, we start from expression (8). A Fourier transform then takes us to with the Fourier transform exp(−D 1 k 2 t) of the Gaussian (1). This expression defines the Laplace transformp(s) of p(D 1 ), to be taken at k 2 t. Fourier inversion of result (10) and substituting κ = kt 1/2 , we arrive at

718
The European Physical Journal Special Topics This shows that we can write P sup = t −d/2 g(ζ 2 ), in terms of the scaling function g(ζ 2 ) that solely depends on the similarity variable ζ = r/t 1/2 . Thus, a given shape of the function g(ζ 2 ) is an invariant, and no transition to a different shape is possible. Crossovers to other shapes, for instance, an effective Gaussian, at long times can be explained in the diffusing-diffusivity framework below.

Anomalous non-Gaussian diffusion
Non-Gaussianity of the PDF P (r, t) is a common feature for anomalous diffusion processes of the continuous time random walk type, in which jumps are interrupted by random waiting times with a scale-free distribution of the form ψ(τ ) τ −1−α such that no characteristic waiting time τ exists [14,36,89]. Instead, stretched Gaussian forms are obtained. Similarly, non-Gaussian shapes of the PDF are known from other anomalous diffusion processes, most prominently for diffusion on fractal supports [90] or diffusion in heterogeneous diffusion processes, in which the diffusion coefficient is explictly position-dependent [91,92]. Experimentally and in simulations, non-Gaussian patterns were observed in membrane dynamics [50,[93][94][95][96], confined diffusion of water molecules in soft environments [97], polymer diffusion along surfaces decorated with nano-pillars [98], intermittent hopping on solid-liquid interfaces [99] similar to bulk-mediated diffusion [100][101][102][103], diffusion of colloids in dense crowded suspensions [104], tracer diffusion in glassy systems [105,106], as well as in mucin hydrogels [107,108], fibrin gels [109], and in disordered micropillar matrices [110]. In simulations, non-Gaussianity was found in crowded and interactive environments [111]. Additionally, non-Gaussian displacement distributions were studied in static disordered media [113] and colloidal liquid crystals [114].
The motion of individual lipid molecules in a bilayer membrane at sufficiently short times was shown to be anomalous-diffusive with displacements that are Gaussian distributed and whose correlations are consistent with the fractional Langevin equation motion defined in Section 1 [49]. 3 However, once the membrane is crowded with embedded proteins (Fig. 6A) the displacement PDF of the lipids can be adequately described by a stretched Gaussian of the form [50] where r is the radial co-ordinate. This functional behaviour is demonstrated in Figure 6B, where the cumulative distribution Π(r, ∆) is plotted in the form − log[1 − Π(r 2 , ∆)] ∼ (r/c∆ α/2 ) δ such that the power-law in the exponent becomes a straight line (Fig. 6, middle). Plotting this behaviour for various measurement times T and lag times ∆ shows that the stretching exponent δ shows only minor variations around values of δ ∈ (1.35, 1.66) (Fig. 6C) [50]. However, it is expected that for times beyond the reach of the simulations normal diffusion with an effective diffusivity will be restored, similar to the observations in the protein-free bilayer membranes [49,51]. If we interpret this non-Gaussian behaviour in the superstatistical language, the stretched Gaussian (12) emerges from a diffusivity distribution of the form p(D α ) ∝ exp(−[cD α ] κ ) such that δ = 2κ/(1 + κ) [72].  shows homogeneous-in-time fluctuations for the dilute bilayer, in the protein-crowded case of Figure 6A, lipids may get associated with the less mobile proteins or become trapped inside a cordon of proteins. The traces in Figure 7A, in particular, the blue trace, show a clear intermittent behaviour. Concurrently, the relatively narrow distribution of diffusivity values K α in the dilute case significantly broadens in the crowded case, with a slightly bimodal shape (Fig. 7B). Remarkably, several of the features of the full protein-crowded membrane system can already be observed in two-dimensional excluded volume systems with narrowly-placed obstacles [50].
In single particle tracking experiments in heterogeneous membranes anomalous diffusion with an almost-exponential displacement PDF and diffusivity distribution were observed [95]. Similarly, in [115,116] the diffusion of submicron tracers in bacteria and yeast cells were demonstrated to show antipersistent motion consistent with the fractional Langevin equation model, however, the displacements showed a Laplace distribution -including an impressive scaling behaviour. The associated diffusivity distribution is, as expected, of exponential form.
To phrase anomalous diffusion in a superstatistical language a generalised Langevin equation model with distributed diffusivity was studied by Beck and van der Straeten [117], while a more general approach for a superstatistical generalised  Langevin equation was introduced byŚlȩzak et al. [118] in which it was shown that the distribution of the position variable is characterised by a relaxation from a Gaussian to a non-Gaussian distribution. Random parameter diffusion models for normal and anomalous diffusion are very actively studied, and we can here only give a limited overview. Apart from the developments sketched above we mention the study by Cherstvy et al. [119] in which scaled Brownian motion for massive and massless particles was analysed for a Rayleigh distribution of the diffusion coefficient. Stylianidou et al. [120] show that in a random barrier model anomalous diffusion with exponential-like step size distribution and anticorrelations emerge, similar to the behaviour measured by Lampo et al. [115], with a crossover to Brownian and Gaussian behaviour at sufficiently long times. Sokolov et al. compare the diffusing-diffusivity model with the emerging dynamics when the quenched nature of a disordered environment is explicitly taken into account [121]. A model for Brownian yet non-Gaussian diffusion based on perpetual multimerisation and dissociation is discussed by Hidalgo-Soria and Barkai [122]. Moreover, we mention a study by Barkai and Burov [123], in which the authors use extreme value statistic arguments to derive a robust exponential shape of the displacement PDF. Finally, in a recent workŚlȩzak et al. [124] show that random coefficient autoregressive processes of the ARMA type can be used to describe Brownian yet non-Gaussian processes, and thus connect the world of physics of such dynamics with the world of time series analysis. From the data analysis side, apart from measuring diffusivity distributions and displacements PDFs, the codifference [125] is a well suited measure to detect non-Gaussianity [126].

Diffusing-diffusivity models
In its formulation, as shown above, superstatistics incorporating the time independent diffusivity distribution p(D 1 ) or p(D α ) cannot account for a crossover to an effective Gaussian PDF at times longer than some correlation time as observed in some experiments [52,53]. An example is shown in Figure 8. In this experiment, colloidal beads diffuse along lipid tubes, the associated displacement PDF clearly shows a crossover from a Laplace-like distribution at earlier times to a Gaussian at later times.
A theoretical model describing this observed crossover behaviour was proposed by Chubinsky and Slater in their "diffusing-diffusivity" model [127]. This approach was further developed by Jain and Sebastian [128,129], Chechkin et al [72], Tyagi and Cherayil [130], Lanoiselée and Grebenkov [131], as well as Sposini et al. [73]. The implementation of the diffusing-diffusivity model in a Bayesian analysis scheme of single trajectory data was investigated in [132]. The basic idea of the diffusing-diffusivity The displacement distribution (C) has exponential tails at earlier times and crosses over to a Gaussian shape at longer times. Reproduced from [52].
picture by Chubinsky and Slater is that the diffusion coefficient in a single trajectory is a stochastic quantity itself, changing its value perpetually along the trajectory of the tracer particle. Physically, this is a simplified picture for a particle moving in a heterogeneous environment, imposing continuous changes in the particle mobility along its path [72,133]. The diffusing-diffusivity dynamics is characterised by an intrinsic correlation time, beyond which the diffusion becomes effectively Gaussian.
In a minimal formulation of the diffusing-diffusivity model, the crossover dynamics can be captured by the set of coupled stochastic equations [72] d dt Here expression (13a) is the Langevin equation for the position r(t) of a particle, driven by the white Gaussian noise ξ(t). Instead of the regular Langevin equation, however, the associated noise strength amplitude contains the explicitly time-dependent diffusion coefficient. This property is specified by equations (13b), that maps D onto the squared auxiliary quantity Y thus guaranteeing positivity of the diffusivity, and (13c). The latter, stochastic equation describes the time evolution of the auxiliary variable Y driven by another white Gaussian noise η(t). However, in contrast to equation (13a), the motion of Y is confined and thus will relax to equilibrium above the crossover time τ -in fact, equation (13c) is the famed Ornstein-Uhlenbeck process [9]. In the analysis of [72] it was shown that this formulation of the diffusing-diffusivity model at short times reproduces the superstatistical approach with exponential tails of the PDF, while at times longer than the correlation time τ of the auxiliary Y process a crossover occurs to a Gaussian PDF characterised by a single, effective diffusion coefficient. This crossover can be conveniently characterised by the kurtosis K = r 4 (t) / r 2 (t) 2 , which reaches the value for a Gaussian distribution at times longer than the Ornstein-Uhlenbeck correlation time τ [72]. More technically, the formulation in terms of the minimal model (13a) to (13c) corresponds to a subordination approach, which is helpful in obtaining exact analytical results and in formulating a two-variable Fokker-Planck equation for the diffusing-diffusivity process [72]. In Figure 9 we show the behaviour encoded in the minimal diffusing-diffusivity model (13a) to (13c). The three panels respectively show the crossover from an initial  Laplace-like distribution with exponential tails to a Gaussian (A) and the fact that the MSD of the process always is linear in time with a constant coefficient (B). Particularly, the crossover behaviour measured by the kurtosis is shown panel (C), indicating the crossover time from exponential to Gaussian shapes, equivalent to the characteristic time scale τ of the Ornstein-Uhlenbeck process (13c). This behaviour is characteristic for the equilibrium nature of the auxiliary variable Y. The more general situation for a non-equilibrium initial condition with crossovers in the associated MSD is analysed in [73]. In particular, for the initial condition D 0 = 0 of the diffusivity D 0 = Y(0) 2 the MSD shows a crossover from initial ballistic scaling proportional to t 2 , to a time-linear scaling at times longer than the correlation time of the Ornstein-Uhlenbeck process Y(t) [73]. We remark that the diffusing-diffusivity model developed here is closely related to the Cox-Ingersoll-Ross (CIR) and Heston models for monetary returns widely used in financial mathematics [134][135][136][137], see the discussion in [72,131].

Conclusions
A growing number of processes from a wide range of systems is being reported in which the measured stochastic motion deviates from the expected Gaussian shape of the displacement PDF. The two most prominent cases are Brownian yet non-Gaussian diffusion, in which the MSD is (a) linear in time (i.e., Fickian) but the PDF is not given by the classical Gaussian form and (b) non-Gaussian stochastic motion driven by fractional Gaussian noise. While in some experiments the non-Gaussianity is seen over the entire time range of observation, in others a crossover to an effective Gaussian behaviour is observed beyond some correlation time. This non-Gaussian behaviour is accompanied by distributions of associated diffusion coefficients.
On a more general level, stochastic processes with random parameters with both stationary-distributed values as in the superstatistical approach or as time-stochastic processes themselves, are important tools in the description of heterogeneous particle "ensembles" or heterogeneous environments. Given such stochastic formulations, follow-up processes such as chemical reactions can be described quantitatively in terms of first-passage formalisms. For the diffusing-diffusivity model the distribution of first-passage times was calculated [133,138]. These, or more general models, may be important for the description of passive diffusion processes in biological cells, that are highly heterogeneous. Current models describing the diffusive search of proteins for specific binding sites on the cell's DNA are mainly based on Brownian and Gaussian diffusion [139][140][141] as well as anomalous diffusion [142] in homogeneous environments. Concurrently, single-trajectory power spectra statistics are being derived for diffusingdiffusivity models [143], extending the theories for single-trajectory power spectra in Brownian motion, FBM, and SBM [144][145][146].
In fact, also completely different behaviours of non-Gaussianity have been observed. As an example, Figure 10 shows the displacement PDF of the twodimensional amoeboid motion of dictyostelium dicoideum cells, projected onto the x-axis [147]. 4 For growing lag time the PDF does not converge to a Gaussian. In contrast, the stretching exponent δ evolves from δ ≈ 1.2 for lag times of few seconds to δ ≈ 1.06 (δ = 1.03 for the y motion) for longer lag times. In other words, that is, the motion becomes more exponential over time. For living cells it appears perfectly reasonable to have a distribution of absolute speeds or "persistence", both effecting a non-Gaussian shape of the displacement PDF. We may speculate whether cells have constant, cell-specific speeds with non-Gaussian distribution, translating into an exponential distribution of their "diffusivity". Further studies are necessary to clarify this point.
We note that there exist similar models with time-varying diffusion parameters, in particular, dichotomous diffusivity models [148,149]. Moreover, random walk models with correlated waiting times [150,151] have been discussed, a variant of which are diffusing waiting times [152][153][154]. For the description of heterogeneous systems other random walk models have also been developed, such as the annealed transit time model [155]. We finally mention a recent result in active random walker systems, in which even non-monotonic displacement distributions were studied [156]. Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.