Quantifying non-ergodicity of anomalous diffusion with higher order moments

Anomalous diffusion is being discovered in a fast growing number of systems. The exact nature of this anomalous diffusion provides important information on the physical laws governing the studied system. One of the central properties analysed for finite particle motion time series is the intrinsic variability of the apparent diffusivity, typically quantified by the ergodicity breaking parameter EB. Here we demonstrate that frequently EB is insufficient to provide a meaningful measure for the observed variability of the data. Instead, important additional information is provided by the higher order moments entering by the skewness and kurtosis. We analyse these quantities for three popular anomalous diffusion models. In particular, we find that even for the Gaussian fractional Brownian motion a significant skewness in the results of physical measurements occurs and needs to be taken into account. Interestingly, the kurtosis and skewness may also provide sensitive estimates of the anomalous diffusion exponent underlying the data. We also derive a new result for the EB parameter of fractional Brownian motion valid for the whole range of the anomalous diffusion parameter. Our results are important for the analysis of anomalous diffusion but also provide new insights into the theory of anomalous stochastic processes.

Superresolution microscopy allows unprecedented insight into the motion of fluorescently labelled, single molecules in complex liquid environments and even inside living biological cells. The observed tracer dynamics also reveals new insights into the physical properties of the systems and thus provides a handle for the modelling of followup processes such as molecular reactions in the system. Concurrently, due to ever increasing computational power, large scale simulations uncover longer and longer time windows of the atomistic or coarse grained dynamics in molecular systems [1][2][3][4][5][6] .
The MSD (1) is obtained as the average of the squared particle position over an ensemble of particles at a fixed time t. In many single particle tracking studies and large scale computer simulations few but long time series x(t) of the tracer particle position are available. These are typically evaluated in terms of the time averaged MSD where Δ is called the lag time and t is the overall length of the time series 4,37 . In the Boltzmann-Khinchin sense we call a stochastic system ergodic when the time averaged MSD (2) converges to the MSD (1) in the long measurement time limit: . For certain anomalous diffusion processes this equality no longer holds, and we observe a so-called weak ergodicity breaking: 4, 37-44 . For finite measurement times even for ergodic processes the time averaged MSD (2) will exhibit more or less pronounced amplitude variations around the mean i N i 2 1 2 taken over N garnered trajectories 4,37,38 . Defining the dimensionless variable ξ δ δ = Δ Δ ( )/ ( ) 2 2 the variations of δ Δ ( ) 2 around δ Δ ( ) 2 are then typically characterised in terms of the variance of ξ, 2 the ergodicity breaking parameter 4,37,38,45,46 . The EB parameter has become a widely used tool to quantify the trajectory-to-trajectory fluctuations in single particle tracking. Note that in some papers EB is defined only in the limit Δ/t → 0, however, we here use it as a measure for the spread ξ also at finite Δ/t ratios. The canonical Brownian motion is characterised by the scaling EB = 4/3 × (Δ/t) of the EB parameter in the limit Δ/t → 0 4, 45 , see ref. 47 for the full expression.
Here we study in detail the full distribution φ(ξ) of the amplitude fluctuations of the time averaged MSD. For the most popular anomalous diffusion processes we demonstrate that while the EB parameter is an important measure for the amplitude of these fluctuations, in many cases it is not sufficient to adequately characterise the distribution φ(ξ). Namely, in many cases φ(ξ) is significantly skewed, that is, asymmetric around its mean. In this case higher order moments should be used to complement the EB parameter. Specifically, we here analyse the skewness where the denominator represents the standard deviation σ = (〈ξ 2 〉 − 〈ξ〉 2 ) 1/2 = EB 1/2 of the amplitude scatter of ξ and we used the fact that 〈ξ〉 = 1 by definition. When γ = 0 the distribution is symmetric, for negative/positive γ it is skewed to the left/right. When |γ| is larger than unity the distribution is considered significantly skewed. As we will see, this is frequently the case for the commonly used anomalous diffusion models, in particular for fractional Brownian motion (FBM), for which a symmetric Gaussian distribution for φ(ξ) was previously suggested. We also investigate in detail the kurtosis which is a measure for the outliers of the scatter distribution φ(ξ). Alternatively in literature also the non-Gaussianity measure is used which involves the fourth time averaged moment 4 . As a reference process the analysis of the TAMSD of normal Brownian motion is discussed in ref. 47 in detail. It can be recovered as the special case α = 1 of our anomalous diffusion models. We note that φ(ξ) provides meaningful information even for relatively sparse data, and it can be reliably used to infer the physical character of the stochastic process underlying the observed data 48,49 .
In the following section II we study the amplitude scatter for FBM at finite values of Δ/t and show, in particular, that the EB parameter varies smoothly as function of the anomalous diffusion exponent α. We also investigate the skewness and kurtosis for FBM and demonstrate the relevance to consider higher order moments. Section III is devoted to the subdiffusive continuous time random walk (CTRW), for which φ(ξ) is naturally skewed, as it has a finite value at ξ = 0. In addition to the standard CTRW we also quantify the shape parameters of φ(ξ) of ageing CTRW processes. Heterogeneous diffusion processes (HDP) with their systematic, quenched variation of the diffusion coefficient are then studied with respect to the skewness and kurtosis in section IV. Finally, in section V we draw our conclusions and argue, why this analysis of additional data inference techniques is important for a reliable quantitative and physical analysis of stochastic data. In the appendix some mathematical details are presented.

Results
Fractional Brownian Motion. FBM is one of the most widely used anomalous diffusion processes. For instance, it is a standard model for stock market dynamics 50 , and it is used as a polymer model 51 . It also describes the effective dynamics of single file diffusion 52 . Moreover, it is commonly used to model single particle diffusion experiments 12,16,17,[53][54][55][56] in living cells. Its physical relevance stems from the fact that it corresponds to the overdamped limit of the fractional Langevin equation associated with particle motion in viscoelastic environments 4,57,58 .
FBM was formulated in 1968 by Mandelbrot and van Ness 59 as a family of Gaussian random functions. We note that a similar process was introduced by Kolmogorov in 1940 60 . FBM is conveniently defined in terms of the Langevin equation where ζ α (t) represents fractional Gaussian noise, defined in terms of the Gaussian yet power-law correlated random noise ζ α (t) of zero mean and two-time correlation for t 1 , t 2 > 0 and t 1 ≠ t 2 45, 61-63 The anomalous diffusion exponent α used here relates to the Hurst exponent, often encountered in the discussion of FBM, via α = 2H. Due to the sign of the factor (α − 1) subdiffusive FBM for 0 < α < 1 is anti-correlated (antipersistent), i.e., more erratic and with a smaller span than Brownian motion. Conversely, positively correlated (persistent) motion occurs in the range 1 < α < 2. In the limit α = 2 the steps are completely correlated and ballistic motion is recovered.
A sample path of FBM is generated analogously to Brownian motion in the form There are several ways to simulate fractional Gaussian noise, including the Hosking method 64 , which we used to generate our simulation paths. Accordingly the path is generated as the sum of a fractional Gaussian noise realisation in that every step is explicitly calculated recursively using the entire path history. The method is defined in detail in ref. 65. Ergodicity breaking parameter. FBM as well as ordinary Brownian motion is known to be ergodic in the 61 . Note, however, that for FBM and the related fractional Langevin equation motion transiently non-ergodicity and ageing effects may occur 18,61,66,67 . Moreover the mean time averaged MSD δ 2 equals the ensemble averaged MSD at all lag times Δ = t 45 . The EB parameter of FBM was proposed to split up in the tripartite scaling law 45 In this expression we see that the EB parameter converges to zero as Δ  t / proportionally to the Brownian case as long as α < 3/2. For α > 3/2 a slower decay Δ α −  t ( / ) 4 2 is observed. We also see that the result (10) includes a divergence at the critical point α = 3/2 45 , as shown in Fig. 1. Here we present the exact analytical result valid for any Δ/t based on a systematic, strictly converging series expansion (see Methods). In particular, we prove that in the limit Δ  t / 1 EB does not display any discontinuity at α = 3/2. To show this let us reconsider the full expression of the EB parameter without approximation 45 , As detailed in Methods, this expression can be modified such that the consistent approximation Δ  t / 1 to leading order yields the FBM EB parameter where the prefactors are defined by This result is clearly different from Eq. (10) and shows that the leading scaling behaviour suggested in ref. 45 is only assumed sufficiently far away from the point α = 3/2, with a quantitative dependence on Δ/t. As α approaches α = 3/2 the EB parameter according to expressions (12), (13), and (14) converges to the same value, and at α = 3/2 a linear term with a logarithmic correction in Δ/t emerges due to the exact limit Scientific RepoRts | 7: 3878 | DOI:10.1038/s41598-017-03712-x Hence, at α = 3/2 the EB parameter assumes the form 16 ln 539 240 (16) The behaviour according to expressions (12), (13), and (14) is demonstrated in Fig. 1 in comparison with the numerically calculated full expression (11) of the EB parameter. Moreover, Fig. 1 confirms that the result of ref. 45 is only valid for α values, that are sufficiently far away from the point α = 3/2. Different values of the ratio Δ/t are shown and exhibit excellent agreement with our results. Only when we consider small values of Δ/t the scaling around the point α = 3/2 is somewhat off the leading order expansion obtained above. More specifically, in the range 0 < α < 3/2 the scaling of EB is Δ  t / , hence the respective prefactor can be calculated by EB/(Δ/t). Analogously for 3/2 < α < 2 the prefactor is EB/(Δ/t) 4−2α and at exactly α = 3/2 it is EB/(Δ/t × lnt). We conclude that our analytical approximation (12) represents the behaviour of the EB parameter of FBM to numerically sufficient accuracy. Figure 2 shows the EB parameter as a function of the anomalous diffusion exponent for different ratios Δ/t of lag time to measurement time. We observe an excellent agreement between our analytical approximation (12) and the numerical result for the full integral (11). In particular, the value of EB is continuous across the value α = 3/2. Further one can see the transition in the scaling. Up to α = 3/2 the curves are equidistant due to the proportionality of EB to Δ/t. For larger values of α the curves progressively approach each other until the ballistic limit α = 2 where they reach the value EB = 2.
The continuity of EB as function of α and the approximative result for EB and its good agreement with the exact result are our first main finding. Figure 1. EB parameter of FBM divided by the respective scaling laws with respect to Δ/t for the regimes 0 < α < 3/2 and 3/2 < α < 2 according to Eq. (10). The behaviour according to Eq. (10) from 45 is represented by the dashed red lines and the red symbol at α = 3/2: it shows a distinct divergence at α = 3/2. Our analytical result (12) based on a less severe approximation is represented by the symbols and the full line corresponding to the numerical evaluation of the integral (11), calculated via the trapezoidal rule with integration step 10 −3 . No divergence remains, and in fact a continuous albeit non-smooth behaviour is revealed. The agreement between our approximation (12) and the numerical evaluation of the full expression (11) remains quite good even for relatively large values of Δ/t. Note the left and right ordinates referring to the cases 0 < α < 3/2 and 3/2 < α < 2.
Scientific RepoRts | 7: 3878 | DOI:10.1038/s41598-017-03712-x Skewness and kurtosis of the amplitude scatter distribution. As stated before, for finite measurement times the time averaged MSD of a single trajectory of any stochastic process is a random variable, whose distribution φ(ξ) provides information about the underlying stochastic process 4,49,68 . The result for φ(ξ) for FBM is depicted in Fig. 3. The top panels show the scatter distribution for different values of the anomalous diffusion exponent α, in the middle and at the bottom rows the scatter distribution is shown for different lag times Δ for the case of subdiffusion and superdiffusion with α = 0.6 and α = 1.6 respectively. The left panels show the data in linear scales whereas the right panels use a semi-logarithmic scale. The tails of the distributions are fitted by exponentials. The special case of Brownian motion with α = 1 is fitted by a Gaussian distribution in the top panel.
For FBM the amplitude scatter distribution φ(ξ) was shown to be approximately Gaussian at sufficiently short lag times for subdiffusive FBM-the underlying argument was based on the assumption that for small Δ any two displacements do not overlap for successive time intervals 49 . It was also shown that φ(ξ) becomes highly asymmetric for larger Δ 49 . As shown in Fig. 3 the Gaussian approximation indeed holds in the subdiffusive regime with 0 < α < 1 for Δ = 1. However, already for Δ = 5 the distribution φ(ξ) shows a pronounced asymmetry. In the superdiffusive case shown in the bottom panel of Fig. 3 this asymmetry is much more pronounced, and even in the case Δ = 1 the shape is obviously far beyond any Gaussian approximation. In all analysed cases the exponential tail appears to capture the behaviour well. We quantify the asymmetry of φ(ξ) by using the skewness parameter (5), which by definition is a function of the ratio of the lag time versus the measurement time as well as the anomalous diffusion exponent α. A symmetric distribution is characterised by a zero-valued skewness. We consider a distribution significantly skewed when the absolute value of the skewness parameter exceeds unity. Depending on the sign of the skewness a distribution can be skewed to the left or right. The skewness parameter evaluated from our simulations is shown in Fig. 4. For the special case of Brownian motion the explicit expression is known from 47 . The comparison to simulated data can be seen in Fig. 5, showing excellent agreement. Although the explicit formula for the skewness, apart from this special case, is unknown, the data can be tentatively fitted to a power-law in the ratio t/Δ of the measurement and lag times (the data are obtained for the fixed lag time Δ = 1) in the form The scaling exponent β and its derivative dβ/dα are shown in the right panel of Fig. 4. For subdiffusive FBM the value of β remains approximately constant at around 1/2 for the entire subdiffusive range 0 < α < 1, leading to a rapid decay of the skewness as function of t/Δ, and thus an approximate validity of the Gaussian distribution for φ(ξ) proposed in ref. 49. In particular, in this analysis the Brownian limit α = 1 does not appear distinguished. However, once we reach the superdiffusive domain 1 < α < 2 the exponent β exhibits a distinct decay to zero. This means that for more superdiffusive FBM the asymmetry of φ(ξ) is a long-ranging characteristic of the process. We also note that the variation of β with the anomalous diffusion exponent becomes most delicate at the value α = 3/2, at which the (Δ/t)-scaling of the EB parameter crosses over, see Eq. (12). Figure 6 further analyses these observations in terms of the kurtosis. For the subdiffusive domain 0 < α < 1 the kurtosis converges relatively quickly to the expected value  = 3 for a Gaussian amplitude scatter distribution. Even for α = 1.2 the deviation is moderate. In contrast to this, for larger values of α the kurtosis converges significantly more slowly and assumes larger values, thus pointing at increasingly large outliers in the superdiffusive domain.
The strong asymmetry of the amplitude scatter distribution for FBM quantified in terms of the skewness and the kurtosis is our second main result.   6 and Δ = 1. The case of Brownian motion with α = 1 is shown with a Gaussian fit. Middle: different lag times for α = 0.6 and t = 2 6 . A Gaussian fit for the short lag time Δ = 1 shows nice agreement with the data, whereas at longer Δ the asymmetry of φ(ξ) becomes obvious in the semi-log scale on the right. Bottom: different lag times for α = 1.6 and t = 2 6 . High asymmetry of φ(ξ) is observed. In the semi-log scale on the right the tails are fitted exponentially. Shlesinger and are jump processes characterised by random waiting times τ in between successive jumps 7,8,[69][70][71] . These τ are assumed to be distributed according to the probability density function ψ(τ). As long as the mean waiting time remains finite, the associated CTRW process renormalises to a discrete time random walk with absolutely sharp distribution ψ(τ) = δ(τ − 〈τ〉) of waiting times for sufficiently long times τ  t 72 . However, once a power-law form ψ τ τ τ τ with 0 < α < 1 is chosen the characteristic waiting time 〈τ〉 diverges and ψ(τ) is scale free. The resulting process describes anomalous diffusion of the form (1). Note that this is consistent with the form ψ τ − α  u u ( ) 1 ( ) 0 in Laplace space 73 which for 0 < α < 1 relates to the characteristic function of a completely one-sided Lévy stable law, and for α = 1 corresponds to ψ(τ) = δ(τ − τ 0 ) with τ 0 = 〈τ〉 72,73 . A power-law form for the waiting time density in terms of τ is sometimes chosen, and in that case the characteristic waiting diverges marginally for α = 1. However, in our formulation based on the Shlesinger-Hughes idea in Laplace space, the limit α = 1 with the finite mean waiting time results consistently.
Power-law waiting time densities are intimately connected with random energy landscapes with exponential distribution of depths 7 , comb models 74 , weakly chaotic and turbulent systems [75][76][77] , as well as dynamic maps [78][79][80] . In experiments they are widely observed 4 . Subdiffusive CTRWs are non-stationary and weakly non-ergodic in the sense that the inequivalence δ Δ ≠ Δ x ( ) ( ) 2 2 remains true at arbitrarily long times 4, 37, 38 . In the limit of many jumps subdiffusive CTRWs are characterised by the amplitude scatter distribution 38 φ ξ α αξ  Figure 7 shows the amplitude scatter distribution φ(ξ) for various α. In particular, the high degree of asymmetry of the shape of φ(ξ) becomes obvious: even for quite large α the maximum of φ(ξ) stays to the right of the ergodic value ξ = 1. In the Brownian limit the distribution (19) converges to a δ-peak, φ(ξ) = δ(ξ − 1) 38, 47 , reflecting the fact that all realisations have the same, fully reproducible amplitude.
Due to the relation between the Mittag-Leffler distribution ρ α (x) and the one-sided Lévy stable law ( with s ≥ 0 and 0 < α < 1, we can rewrite relation (19) as t Using the moment-generating function M ξ (s) = E α (−sΓ(1 + α)) the moments can be generated by Consequently the EB parameter is given by  Figure 8 shows the EB parameter, the skewness, and the kurtosis of the subdiffusive CTRW as a function of the anomalous diffusion exponent α. The EB parameter shows a relatively moderate variation over the interval α∈(0,1] and decays monotonically from EB = 1 at α → 0 to EB = 0 in the Brownian case α = 1. In the Brownian limit α = 1 corresponding to Brownian motion with finite moments 〈τ k 〉 of the waiting times the skewness and kurtosis reduce to γ = 0 and  = 3, respectively. The approach to these Gaussian values for finite t/Δ is detailed in Fig. 9. Note once more that according to our choice for the waiting time density starting from the expression in Laplace space the form ψ(τ) = δ(τ − τ 0 ) emerges in the Brownian limit, for which all moments are indeed finite.
Let us now turn to the subdiffusive case with 0 < α < 1. The skewness γ in Fig. 8 decays monotonically from γ = 2 in the limit α → 0 and generally has a more complex form than EB. γ changes its sign at around α = 0.78 meaning that the scatter distribution switches from a right-skewed to a left-skewed distribution. For most of the range of the anomalous diffusion exponent the absolute value of the skewness is greater than unity-the two arrows in Fig. 8 indicate where γ crosses unity-and thus significantly skewed. In this region the EB parameter is thus not a sufficient quantity to characterise the shape of the amplitude scatter distribution φ(ξ). When we gradually approach the Brownian case α = 1 we see that γ diverges to negative values and apparently does not converge to its zero Brownian value. The reason for this lies in the very shape of φ(ξ) shown in Fig. 7. Namely, the maximum of φ(ξ) is not centred at the ergodic value ξ = 1. This off-centre behaviour persists until we fully reach the Brownian case α = 1. When φ(ξ) is sufficiently sharp such that we can evaluate the moments as if they were generated by a δ-function positioned at ξ = 1 + ε, we obtain 〈ξ k 〉 = (1 + ε) k . In the limit ε → 0 the skewness (5) indeed diverges to minus infinity. This heuristic argument is consistent with the limiting behaviour of expression (26) when we approach the Brownian case: In contrast to the EB parameter, the skewness thus exhibits a discontinuous behaviour when we approach the Brownian limit. This is an important observation with respect to the analysis of subdiffusive CTRW processes. Note that in experiments values of α ≈ 0.9 are indeed observed and significant 26 . for Brownian motion 47 can be clearly anticipated.

Figure 9.
Skewness and kurtosis for a discrete-time random walk with Gaussian jump length distribution. Solid lines represent the analytic result for Brownian motion from 47 . The convergence to γ = 0 and = 3  is nicely displayed by the data. While the skewness exhibits a power-law decay to zero the kurtosis appears to converge to its value much earlier and then does not change much any more.
Similarly, the kurtosis also varies from an intermediate value for small α over a minimum, reaching quite large values when the anomalous diffusion exponent approaches unity. This divergent behaviour and the discontinuity at the point α = 1 can be rationalised analogously to the discussion of the skewness.
In comparison to FBM we note that for the subdiffusive CTRW we consider the shape parameters at t → ∞. Interestingly for both the EB parameter and the skewness the general trend is opposite: for FBM the smallest values occur for low α, while for the CTRW this corresponds to the largest values (not taken the divergencies for larger α into account). The higher order moments may therefore be good inference indicators for the CTRW mechanism. We also note that due to the larger variation of the skewness with α it may provide a useful tool to extract the anomalous diffusion exponent from sufficiently good data. The delicate variation of the kurtosis around α ≤ 1 makes it a particularly good indicator for α in this region.
The exact variation of γ and  with α and, in particular, their apparent divergence for α → 1 and discontinuous behaviour at α = 1 are our third main result.
Ageing CTRW. Often, the start of a measurement does not coincide with the initiation of the time evolution of the system. For instance, when we measure the motion of an endogenous lipid granule in a living cell, it is not clear what the starting point for this process is. Or, we could probe the charge current in amorphous semiconductors and on purpose introduce a time shift between the creation of the charge carriers and the moment when we switch on the driving electrical field 81,82 . If the dynamics is stationary, a delay between system initiation and start of the measurement does not make any difference, as the correlation functions solely depend on the time difference. In a non-stationary system, in contrast, we observe ageing phenomena 4 , similar to those in glassy systems 83,84 . The delay between initiation and start of the measurement is called the ageing time t a 4 . Subdiffusive CTRWs represent a prototype case to study ageing effects [85][86][87][88][89] . It was shown that ageing in subdiffusive CTRWs leads to a population splitting into a fraction of particles, that never move during an observation period of duration t starting after the system has been allowed to age for the period t a -and a fraction of variable mobility 88,89 . For the aged process the time averaged MSD is calculated via 88,89 The distribution φ(ξ) for the ageing CTRW then yields in the form 88,89 φ ξ δ ξ α where m α is the probability to observe at least one jump during the observation period t 88,89 , In analogy to the procedure outlined in Methods we determine the moments of ξ via the Mellin transform of the scatter distribution (29), yielding Figure 10. EB parameter, skewness, and kurtosis of ageing CTRW as a function of the anomalous diffusion exponent α in case of different ratios t a /t of the ageing time-to-measurement time. For comparison, we also include the behaviours of the parameters in the non-aged regime.
Scientific RepoRts | 7: 3878 | DOI:10.1038/s41598-017-03712-x Figure 10 shows the EB parameter, skewness, and kurtosis for CTRW for different ageing time-to-measurement time ratios. When the ageing time becomes more dominant we observe an increase of the degree of non-ergodicity, due to the increasing inhomogeneity of the jump statistic following the above-mentioned population splitting 88,89 . Except for a small parameter range of larger α the distribution is significantly asymmetric, as shown by the skewness parameter. Also outliers are more pronounced for an aged system, as can be seen for the kurtosis. These results for the shape parameters of the ageing process complete our study of CTRWs.

Heterogeneous diffusion process and gamma-distributed amplitude scatter distributions.
Another important class of anomalous stochastic processes are heterogeneous diffusion processes (HDPs), which may be defined in terms of the multiplicative Langevin equation 91 in which ζ(t) represents Gaussian, delta-correlated noise with zero mean. We here briefly discuss the properties of HDPs for completeness. One of the well studied cases for the (quenched) position dependent diffusion coefficient is the power-law form 91 where D 0 has the dimension cm 2−β sec −1 . The small shift x off = 10 −2 prevents the particle from being trapped at x = 0. In the Stratonovich interpretation the process was studied in detail in refs 91-93. It is known that the Stratonovich interpretation of the Langevin equation produces a so-called spurious drift, whereas a thermodynamically consistent interpretation follows from the 'isothermal' discretisation [94][95][96] . Note that rigorous results on the general symmetries of overdamped multiplicative white noise processes, clarifying some prevailing misconceptions on the Itô-Stratonovich dilemma and presenting an exact mapping of the different interpretations to a purely additive white noise process, were reported only recently 97 . Nevertheless, the notion of a 'correct' interpretation of multiplicative white noise Langevin equations is always context dependent 98-100 , as it was shown, for example, that general non-linear relaxation processes not satisfying a fluctuation-dissipation theorem in fact obey a Fokker-Planck equation following from the Itô interpretation of the underlying Langevin equation 101 . Moreover, as we here advocate the general importance of higher-order statistics in the quantification of stationary and/or ergodic properties of single-particle time series data, the specific interpretation of the multiplicative white noise Langevin equation is not critical. A particle in the diffusivity field D(x) is effectively pushed into regions of ever smaller diffusivity. The MSD of the process follows the power-law form (1) with α = 2/(2 − β) 91 . Figure 11 shows two representative probability density functions of the HDP amplitude scatter generated by the simulation scheme This probability density function can be fitted by the generalised Gamma distribution 91 with three fitting parameters a, b and ν and ξ > 0. Here K ν is the modified Bessel function of the second kind 90 . This distribution has exponential cutoffs for both small and large argument, as specified by the parameters a and b. In between, it exhibits power-law scaling with exponent ν − 1. This form is similar to generic first passage densities of stochastic processes 102 . The above density has the moments 103 From which the EB parameter, the skewness, and the kurtosis can be readily constructed. Investigating the respective quantities numerically shows that the skewness can reach values larger than unity and the kurtosis may attain very large values, pointing at significant outliers in φ(ξ).

Discussion
We presented an analysis of the amplitude scatter distribution of the time averaged MSD of various, widely used anomalous diffusion processes: FBM is an ergodic process in the sense that time and ensemble averages of physical observables converge to each other. Subdiffusive CTRW and HDPs are weakly non-ergodic and show a disparity between time and ensemble averages 4 . In all cases we identified relevant situations when the amplitude scatter distribution is not sufficiently characterised by the ergodicity breaking parameter EB due to the significant skewness of the distribution. A rough criterion for when the skewness is non-negligible is given when the absolute value of the skewness exceeds unity. As important additional parameters for the quantitative description of the amplitude scatter of the time averaged MSD, we analysed this skewness as well as the kurtosis of these processes. We believe that the results presented here are important additional characteristics for the analysis of anomalous stochastic time series measured in experiments and simulations in complex systems. Moreover they demonstrate the importance of the higher order moments of φ(ξ) in the analysis of data garnered from experiments or simulations. The amplitude scatter distribution φ(ξ) is a very useful quantity to classify anomalous stochastic processes. Higher order moments, as shown here, provide even more relevant information for the physical classification and interpretation of stochastic data. We note that, indeed, following better data from experiments, higher order moments provide useful information for various processes, for instance, shown in a recent study of enzyme kinetics 104 .
Specifically for FBM we revisited the EB parameter and showed that with a less severe approximation EB is continuous over the entire range of the anomalous diffusion exponent, 0 < α < 2, and no singularity at α = 3/2 occurs. We showed that the previous results are in good agreement with the ones found here sufficiently far away from α = 3/2. We also showed that while a previous assumption that the amplitude scatter distribution of FBM is Gaussian for short lag times is approximately correct for subdiffusive FBM, for anomalous diffusion exponents in the superdiffusive range strong asymmetries of the scatter distribution arise for which the Gaussian description becomes inadequate and the skewness therefore is an important parameter to quantify the distribution. Cognisance of this fact is particularly important for generalisations of first passage processes in chemical and molecular biological contexts [105][106][107] , especially when active modes are present 108,109 . Note that indeed in ref. 34 the superdiffusive transport inside amoeba cells was shown to be of FBM form.
As parameter inference from stochastic time series is being increasingly recognised as an important field in the theory of stochastic processes, analyses demonstrate that it is vital to compare several complementary measures for a given time series, in order to identify the very physical nature of the underlying stochastic mechanism (FBM, HDP, CTRW, etc.) and correctly infer the parameters in this process 4 . Apart from the ensemble and time averaged MSDs, the EB parameter as well as the amplitude scatter distribution φ(ξ) and its shape parameters γ and  analysed in this work we mention ratios of ensemble moments and the mean maximal excursion method 110 , fractional order moments 33 , the p-variation method [54][55][56] , Bayesian ranking 111 , or ageing analysis 13,14,[26][27][28][29] , see also the review 112 on analysis of diffusive motion. Generally, using several complementary methods improves the quality of the data inference. As we show here γ and  are relevant parameters to distinguish stochastic processes from each other but also to extract good estimates for the anomalous diffusion exponent α itself, due to their high α-sensitivity.
An alternative approach is to use the time averaged van Hove cross-correlation functional 22 Here 4πχ 2 G tt (χ,Δ) is the probability density to find the particle at time Δ + t′ at a relative displacement between χ and χ + dχ. This cross-correlation function essentially describes the frequency of occurrences of jump lengths χ along the trajectory during the lag time delta. Analysing this distribution for various values of Δ provides additional information about the process. Especially the physical aspects underlying the bimodality of the scatter distribution for aged CTRW processes might be visible in an analysis of the van Hove cross-correlation function, an aspect requiring substantial simulations work in the future. We expand the integrands using the generalised binomial series n n 0 which converges absolutely for |x| < 1. In order to assure an absolute convergence for all Δ/t as well as all α, we rewrite the integral as Explicitly, Eq. (11) becomes We interchange the integral and the sums and take the limit ε → 0. Note that we have not made any approximation up to this point. Focusing on Δ  t / 1 we only keep terms of leading order. As we show below, the special result for α = 3/2 can be obtained smoothly from this result by taking the limit α → 3/2 from below and from above. For α ≠ 3/2 after performing the integrals we obtain     For the reason of clarity we decided to use the approximation above since it shows excellent agreement with the numerical results.
We note in passing that a direct use of Mathematica to evaluate the integrals leads to problematic results, as the necessary analytic continuations for the involved special functions are neglected.
Continuous time random walk: calculation of the moments and skewness in the Brownian limit. We here consider an alternative derivation of the amplitude scatter distribution and its moments for subdiffusive CTRW in terms of a Fox H-function. Identifying the Laplace image of a one-sided Lévy stable density l α (t) with the corresponding H-function yields 113 1, 1 1,0 Given the argument in l α and the prefactor in Eq. (19), this means that