Anomalous diffusion and ergodicity breaking in heterogeneous diffusion processes

Since the early systematic studies of Perrin and Nordlund [1], single particle tracking has become a routine method for tracing individual particles in microscopic systems [2], but also for animal [3] and human [4] motion patterns. In a wide range of systems and over many time and length scales, these measurements reveal anomalous diffusion with mean squared displacement (MSD) 〈x²(t)〉 ≃ t^p [5]: subdiffusion (0 < p < 1) was found for the passive motion of tracers in living biological cells [6–11], a field in which single particle tracking has proved to be particularly useful, or for bacteria in biofilms [12]. A comprehensive review on anomalous diffusion in the crowded environment of living cells was recently published [13]. More generally, subdiffusion is observed in a large variety of systems, including charge carrier transport in amorphous semiconductors [14], tracer dispersion in groundwater [15] or in fractal systems [16, 17]. Superdiffusion (p > 1) is observed for active motion in the biological cells [11] of mussels [18], plant lice [19] or higher animals and humans [3, 4].

Given the broad use of single particle tracking methods to garner diffusion data in both micro- and macroscopic systems and thus the need to evaluate and physically interpret the recorded time series x(t) of the particle position, the ergodic properties of a stochastic process are of particular current interest: is the information obtained from time averages of a single or few trajectories representative for the ensemble? In other words, are sufficiently long time averages reproducible and equivalent to the corresponding ensemble averages [20, 21]? This question has motivated some vivid activity in the field of statistical physics recently. Thus, of the various stochastic models for anomalous diffusion, diffusion on random fractals has been shown to be ergodic [22]. Fractional Brownian and fractional Langevin equation motion driven by long-range correlated Gaussian noise are transiently non-ergodic [23–26]. In contrast, subdiffusive Scher–Montroll continuous time random walks (CTRW) [14, 27] with diverging characteristic waiting time exhibit weak ergodicity breaking [28], and time averages remain random even in the long time limit [29]. Weak ergodicity breaking was indeed observed for the blinking dynamics of quantum dots [30], protein motion in human cell walls [7], lipid granule diffusion in yeast cells [8] and of insulin granules in MIN6 cells [9].

While CTRW processes represent an elegant tool to describe anomalous diffusion, the occurrence of power-law distributed waiting times is far from universal. We therefore ask the question of whether we can create a simple, alternative stochastic model that still features this intricate, weakly non-ergodic behaviour. Here we show that the MSD of Markovian, heterogeneous diffusion processes (HDPs) with space dependent diffusion constant D(x) ≃ |x|^α scales like 〈x²(t)〉 ≃ t^p with p = 2/(2 − α), while the time averaged MSD $\overline {\delta ^2}$ scales linearly in both sub- (α < 0) and superdiffusive (α > 0) regimes. We quantify the non-ergodicity in terms of the ergodicity breaking parameters and the amplitude scatter distribution of $\overline {\delta ^2}$. We showcase the analogies and differences between HDPs and other non-ergodic diffusion processes.

Physically, a space dependent diffusivity appears a natural description for diffusion in heterogeneous systems. Examples include Richardson diffusion in turbulence [31] as well as mesoscopic approaches to transport in heterogenous porous media [32] and on random fractals [33]. Recently, maps of the local cytoplasmic diffusivities in bacterial and eukaryotic cells showed a highly heterogeneous landscape [34], recalling the strongly time-varying diffusion coefficients of tracers in cells, see e.g. [35].

We start with the Markovian Langevin equation for the displacement x(t) of a test particle in a heterogeneous medium with space dependent diffusivity D(x), namely,

$$\begin{eqnarray*} &&\frac{{\rm d}x}{{\rm d}t}=\sqrt {2D(x)}\zeta (t). \end{eqnarray*}$$

Here, ζ(t) is white Gaussian noise with δ-correlation 〈ζ(t)ζ(t')〉 = δ(t − t') and zero mean 〈ζ(t)〉 = 0. In the following we interpret the Langevin equation in the Stratonovich sense [36], to ensure the correct limiting transition for the noise with infinitely short correlation times [37]⁵. Specifically, for D(x) we employ the power-law form [38, 39]⁶

$$\begin{equation} D(x)=D_0|x|^{\alpha}. \end{equation} \tag{ 1 }$$

The dimension of D₀ is [D₀] = cm^2−α s⁻¹. In the Stratonovich interpretation, we introduce the substitution $y=\int ^x[2D(x')]^{ -1/2}\,\mathrm {d}x'$, where y(t) is the Wiener process [36]. For the initial condition P(x,0) = δ(x) a compressed (α < 0) or stretched (α > 0) Gaussian probability density function (PDF)

$$\begin{equation} P(x,t)=\frac{|x|^{1/p-1}}{\sqrt{4\pi D_0t}}\exp\!\left(-\frac{|x|^{2/p}}{ (2/p)^{2}D_0t}\right) \end{equation} \tag{ 2 }$$

emerges (see appendix A), with p = 2/(2 − α). The ensemble averaged MSD $\langle x^{2}(t) \rangle =\int x^2P(x,t)\,\mathrm {d}x$ becomes

$$\begin{equation} \langle x^{2}(t)\rangle=\frac{\Gamma(p+1/2)}{\pi^{1/2}}\left(\frac{2}{p}\right) ^{2p}(D_0t)^p. \end{equation} \tag{ 3 }$$

Thus, for α < 0 we find subdiffusion, and superdiffusion for α > 0. Brownian motion with 〈x²(t)〉 = 2D₀t emerges for α = 0, and α = 1 produces ballistic motion. The diffusion becomes increasingly fast when α tends to 2.

Figure 1 compares the PDF (2) to simulations results (see below). Compared to the Gaussian shape of Brownian motion (α = 0), we observe a depletion in regions of higher diffusivity. For subdiffusion, this causes the central dip of the PDF, while for superdiffusion probability is shifted towards the origin. Note that the shape of the PDFs is significantly different from those of CTRW processes. In particular, for the subdiffusive case, the CTRW-PDF has a pronounced cusp at the origin [5]. Curiously, the subdiffusive shape in figure 1 resembles the propagator for retarded wave motion [40].

Zoom In Zoom Out Reset image size

To characterize the HDP, we calculate the position autocorrelation (see appendix A) for the case t₂ > t₁,

$$\begin{eqnarray} &&\fl \langle x(t_1)x(t_2)\rangle = \frac{2^{p+1}}{\sqrt{\pi}}\frac{\Gamma(p+1) \Gamma(\frac{p}{2}+1)}{p^{2p}\Gamma(\frac{p}{2}+\frac{1}{2})}[D_0t_1]^{(p+1)/2}\nonumber\\ &&\times[D_0(t_2-t_1)]^{(p-1)/2}\,_2F_1\left(\frac{1-p}{2},\frac{p}{2}+1; \frac{3}{2};\frac{-t_{1}}{t_2-t_1}\right) \end{eqnarray} \tag{ 4 }$$

with the hypergeometric function ₂F₁(z). The Brownian limit 〈x(t)x(t + τ)〉 = 2D₀t follows for α = 0, while for t₁ = t₂ we recover equation (3). The asymptotic behaviour of expression (4), 〈x(t)x(t + τ)〉 ≃ τ^(p−1)/2t^(p+1)/2 for τ/t ≫ 1 implies that the correlations decay with τ for subdiffusive processes whereas they increase in the case of superdiffusion, similar to fractional Brownian motion (FBM) [23]. We also obtain the correlation of two consecutive increments along x(t). For τ/t ≪ 1, we find the simple expressions for p = 1/2 and 2:

$$\begin{equation} \langle[x(t)-x(t-\tau)][x(t+\tau)-x(t)]\rangle\sim\left\{\begin{array}{@{}ll} -\tau\sqrt{D_0/t}, & p=\frac{1}{2},\\ (D_0\tau)^2, & p=2.\end{array}\right. \end{equation} \tag{ 5 }$$

Indeed the occurrence of antipersistence (negative correlations) for subdiffusion and persistence for superdiffusion holds for all values of p (appendix A), again similar to FBM. We present exact results for the velocity autocorrelation function in appendix A.

To connect to single particle tracking experiments we now turn to the time averaged MSD of a trajectory x(t),

$$\begin{equation} \left\langle\,\overline{\delta^2(\Delta)}\,\right\rangle=\frac{1}{T-\Delta}\int\nolimits_{0}^{ T-\Delta}\left<\left[x(t+\Delta)-x(t)\right]^2\right>\mathrm{d}t, \end{equation} \tag{ 6 }$$

where Δ is the lag time and T the length of the time series x(t). In equation (6) we introduced the additional average over individual trajectories, $\langle \overline {\delta ^2(\Delta )}\rangle =N^{-1}\sum _{i=1}^N\overline {\delta ^2_i(\Delta )}$ [29]. For Δ ≪ T we find the linear Δ-dependence (compare appendix B)

$$\begin{equation} \left\langle\,\overline{\delta^2(\Delta)}\,\right\rangle =\frac{\Gamma(p+1/2)}{\pi^{1/2}} \left(\frac{2}{p}\right)^{2p}D_0^p\frac{\Delta}{T^{1-p}}. \end{equation} \tag{ 7 }$$

Remarkably, this result holds for both sub- and superdiffusive regimes, and we find the exact match

$$\begin{eqnarray*} &&\left\langle\,\overline {\delta ^2(\Delta )}\,\right\rangle =(\Delta /T)^{1-p}\langle x^2(\Delta )\rangle. \end{eqnarray*}$$

Equation (7) is our central result. Despite the simple nature of the Markovian HDP, we find that it displays weak non-ergodicity, i.e. the scaling of time and ensemble averages is different, $\langle \overline {\delta ^2(\Delta )}\rangle \neq \langle x^2(\Delta )\rangle$. Linear scaling of the type of equation (7) was found for subdiffusive CTRWs [29], but also for correlated CTRWs [41], ageing, non-renewal CTRW processes [42] and time-scaled Brownian motion [43]. It contrasts the ultraweak non-ergodicity of superdiffusive CTRWs [44]. In particular, the dependence on the length T of the time series is identical to CTRW and correlated CTRW subdiffusion: with increasing T the effective diffusivity D_eff ≃ T^p−1 decreases. In the superdiffusive regime, this D_eff increases with time, the particle becomes more mobile. For HDP processes the ageing dependence on T arises when the particle continues to venture into regions of changing diffusivities, for subdiffusive CTRWs this is due to increasingly longer sojourn times.

Exploring HDPs in more detail numerically, we note that the Langevin equation in the Stratonovich interpretation leads to an implicit mid-point iterative scheme for the particle displacement. At step i + 1, we thus have $x_{i+1}-x_i=\sqrt {2D([x_{i+1}+x_i]/2)}(y_{i+1}-y_{i})$, where the increments (y_i+1 − y_i) of the Wiener process represent a δ-correlated Gaussian noise with unit variance. Unit time intervals separate consecutive iteration steps. In our numerical study we consider the two generic cases p = 1/2 and 2. To avoid divergencies in the discrete scheme, for subdiffusion we regularize the diffusivity at x = 0, choosing D(x) = D₀A/(A + x²), such that for x² ≫ A we recover the original relation (1). For superdiffusion, to avoid trapping at x = 0 and thus prohibitively expensive simulations, we use the shifted form D(x) = D₀(|x| + 1), again with the correct scaling for large |x|.

From the generated trajectories, we compute the PDF (figure 1) as well as the ensemble and time averaged MSDs. The MSDs are shown in figure 2. For the ensemble averaged MSD (blue curve) we observe excellent agreement with equation (3). In the subdiffusion case, the predicted behaviour (3) is reached after less than a dozen steps, such that the rectification of D(x) and the choice of the off-centre initial position x₀ are practically negligible. For superdiffusion, an extended region of almost normal diffusion is observed for small particle displacements, as long as |x| < 1 in the rectified D(x). The long time behaviour shows good agreement between theory and simulations. The time averaged MSDs and their mean $\langle \overline {\delta ^2}\rangle$ follow nicely the predicted linear scaling. Only when the lag time Δ approaches T the linear behaviour levels off due to the finite trajectory length T.

Zoom In Zoom Out Reset image size

Next we study the apparent diffusion coefficients and scaling exponents from the individual $\overline {\delta ^2}$. The points of each trace $\overline {\delta ^2}$ were logarithmically binned. The initial part of the trajectory, Δ = 1–10² steps, was fitted by a linear law with diffusivity D₁. The long-time regime, Δ ∼ 10²–10⁴ steps, was fitted with two parameters, $\overline {\delta ^2}=D_{\beta }\Delta ^{\beta }$. The distributions of D₁, D_β and β were then fitted with the three-parameter gamma distribution, g(z) ∼ z^ν−1 e^−z/b e^−a/z, the parameters being fixed by the first three moments obtained from the histograms. This gamma distribution was recently proposed for the spread with Δ of $\overline {\delta ^2}$ of a Brownian walker [45]. In addition, we extracted the amplitude scatter PDF $\phi (\xi =\overline {\delta ^2}/\langle \overline {\delta ^2}\rangle )$ of individual $\overline {\delta ^2}$ around the mean $\langle \overline {\delta ^2}\rangle$.

Figures 3 and 4 show the results of this analysis, together with the fits to the gamma distribution, g(z). We see that the amplitude scatter shows a pronounced spread around the ergodic value ξ = 1 in both sub- and superdiffusive regimes. The shape of ϕ(ξ) changes only marginally with Δ in the linear region of $\overline {\delta ^2}$. Broad distributions of ϕ(ξ) are known from both sub- and superdiffusive CTRWs. The differences are, however, that for subdiffusive CTRWs ϕ(ξ) has a finite value at ξ = 0 due to completely immobilized particles [29]. For superdiffusive CTRWs ϕ(ξ) decays to zero at ξ = 0, but changes significantly with Δ [44]. The scatter distribution ϕ(ξ) is thus a good criterion to distinguish HDPs and CTRWs.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

More specifically, for subdiffusion the distributions of D₁ and the amplitude scatter ϕ(ξ) are wider, while for superdiffusion they show a sharper peak with a maximum shifted to the mean at ξ = 1. The width of ϕ(ξ) is roughly unity for subdiffusion and 1/2 for superdiffusion. It changes only slightly with T in the range T = 10³–10⁵, see figure 3. The distributions for D_β typically exhibit a long power-law tail with f(D_β) ∼ D^−1.2...1.4_β. The distributions for the apparent long-time exponent β are similar for sub- and superdiffusion, with 〈β〉_p=1/2 ≈ 0.86 and 〈β〉_p=2 ≈ 0.92. Finite T-effects lead to the undershoot of 〈β〉 compared to the theoretical value of unity.

Consider now the ergodicity breaking parameter [29]

$$\begin{equation} \mathrm{EB}=\lim_{T\to\infty}\frac{\left(\left\langle (\overline{\delta^2})^2 \right\rangle -\left\langle\,\overline{\delta^2}\,\right\rangle^2\right)}{\left\langle\,\overline{\delta^2}\,\right\rangle^2}\end{equation} \tag{ 8 }$$

that provides a good measure for the dispersion of time averaged MSDs for different types of diffusion processes [46]. The sufficient condition of ergodicity of a process is EB = 0. A necessary condition is that the ratio of the time and ensemble averaged MSD equals unity. Here we observe that this condition is not fulfilled for p ≠ 1,

$$\begin{eqnarray*} &&\mathcal {EB}=\left\langle\,\overline {\delta ^2(\Delta )}\,\right\rangle/\langle x^2(\Delta )\rangle =(\Delta /T)^{1-p}. \end{eqnarray*}$$

The ergodicity breaking parameters EB and $\mathcal {EB}$ extracted from simulations data at Δ/T ≪ 1 is EB ≈ 0.4 for sub- and EB ≈ 1.4 for superdiffusion, indicating the weakly non-ergodic nature of HDP processes in both regimes. The parameter $\mathcal {EB}$ follows the predicted scaling for large Δ and T values, i.e. when both ensemble and time averaged MSDs have converged to the theoretical scaling (figures 2 and 5). For Brownian diffusion our simulations yield the correct finite-time scaling $\mathrm {EB}(\Delta )=\frac {4}{3}\frac {\Delta }{T}$ and sharp amplitude scatter ϕ(ξ) ∼ δ(ξ − 1) at T → ∞, indicating the ergodic behaviour for long traces and the self-averaging property of normal diffusion. For p = 2 at Δ/T → 0 the estimated theoretical value is EB = 2/3 (equation (C.1)) close to the value from simulations.

Zoom In Zoom Out Reset image size

The simulations show that for subdiffusion, $\mathcal {EB}$ increases with Δ/T, while for superdiffusion it decreases with Δ/T (figure 5). In both cases $\mathcal {EB}$ approaches unity as Δ → T. At Δ/T ≪ 1, $\mathcal {EB}$ for p < 1 assumes smaller values, indicating weaker deviations from ergodicity, in contrast to superdiffusion, see figure 5. The parameter EB defined via the fourth moment is overall a more robust characteristic of the process. For instance, it varies only slightly with T and Δ for Δ/T ≪ 1 (not shown).

Anomalous diffusion is a widely observed phenomenon across many scales and disciplines. Following the dramatic increase of single particle tracking studies, the question arises of how we can physically interpret the recorded trajectories in terms of time averages of observables such as the MSD. Namely, in most cases the available physical theories provide results for the ensemble averages of physical observables. Once a system exhibits weak ergodicity breaking, however, the inequivalence of ensemble from time averages prohibits the application of such ensemble theories to the measured, time-averaged quantities. The quantitative interpretation of time averages therefore requires knowledge about whether the system is ergodic and the standard results can be used to fit the data, or whether the system is weakly non-ergodic. In that case knowledge of the time averaged observables is indispensable.

Here we discussed the seemingly simple Markovian HDP with power-law space-dependent diffusion coefficient, defined in terms of a Langevin equation, which is fully local in both space and time. Despite this locality HDPs are weakly non-ergodic: the time averaged MSD $\overline {\delta ^2 (\Delta )}$ scales linearly with the lag time Δ for both sub- and superdiffusion. Such a behaviour was previously known for subdiffusion from more complex processes: CTRWs with diverging characteristic waiting time [29], CTRWs with correlated waiting times [41], non-renewal, ageing CTRWs [42] and time-scaled Brownian motion [43]. For superdiffusion only an ultraweak ergodicity breaking was reported in which the coefficients of the time and ensemble MSDs differ [44], while space-correlated CTRWs feature a Richardson-type hyperdiffusive scaling of the time averaged MSD [41].

The amplitudes of $\overline {\delta ^2 (\Delta )}$ of different trajectories show a broad and asymmetric distribution, another sign of non-ergodicity. In contrast to CTRW subdiffusion, however, in HDPs there is no contribution from completely immobile trajectories. Concurrently, HDPs exhibit (anti)persistence of the increment correlation functions, such that for subdiffusion (superdiffusion) subsequent increments preferentially have opposite (equal) direction. This property is similar to the noise correlator of FBM, moreover the velocity autocorrelator of both HDP and FBM processes is strikingly similar.

HDPs represent new tools for the description of weakly non-ergodic dynamics. Due to their intuitive formulation in terms of space-dependent diffusivities this dynamic behaviour of HDPs directly follows from the physical properties of the environment. HDPs may be distinguished from other processes due to the difference in the amplitude scatter distribution of time averaged MSDs, the increment autocorrelation function, as well as the two ergodicity breaking parameters. It will be interesting to see whether two-dimensional generalizations of the HDP model may grasp some of the interesting dynamical patterns observed in single particle motion in biological cells.

How general are the results obtained herein for power-law forms of the position dependent diffusion coefficient? A preliminary numerical analysis shows that similar forms of weak ergodicity breaking is observed when we consider diffusivities with both stronger and weaker x-dependence, for instance, exponential or logarithmic forms. Details will be presented in a separate publication.

We thank E Barkai and H Buening for discussions and acknowledge funding from the Academy of Finland (FiDiPro scheme to RM) and Deutsche Forschungsgemeinschaft (grant no. CH 707/5-1 to AGC).

To solve the Langevin equation $\dot {x}(t)=\sqrt {2D(x)}\zeta (t)$ in the Stratonovich sense we introduce the new variable [36]

$$\begin{equation} y(x)=\int^x\frac{\mathrm{d}x'}{\sqrt{2D(x')}}=\frac{\sqrt{2/A}}{2-\alpha} |x|^{(2-\alpha)/2}\,\mathrm{sign}(x), \end{equation} \tag{ A.1 }$$

where y(t) is the Wiener process with the known Gaussian PDF for the initial condition y(0) = 0, namely,

$$\begin{equation} P(y,t)=\frac{1}{\sqrt{2\pi t}}\exp\left(-\frac{y^2}{2t}\right). \end{equation} \tag{ A.2 }$$

Returning to the x-variable yields the PDF (2).

With the back-transformation y → x, from equation (A.1) we find that x(t) = 2^p/2p^−p|y|^psign(y), and the autocorrelation function is evaluated via the two-point PDF for the Wiener process

$$\begin{eqnarray} &&\fl \left\langle x(t_1)x(t_2)\right\rangle=\frac{2^p}{p^{2p}}\int_{-\infty}^{ \infty}\mathrm{d}y_1\int_{-\infty}^{\infty}\mathrm{d}y_2P(y_1,t_1)\mathrm{sign}(y_1)\mathrm{sign}(y_2)|y_1|^p|y_2|^p \pi(y_2,t_2|y_1,t_1), \end{eqnarray} \tag{ A.3 }$$

where we used the transition probability

$$\begin{equation} \pi(y_2,t_2|y_1,t_1)=\frac{1}{\sqrt{2\pi(t_2-t_1)}}\exp\left(-\frac{(y_2-y_1) ^2}{2(t_2-t_1)}\right). \end{equation} \tag{ A.4 }$$

Evaluating the integral in equation (A.3) we arrive at equation (4).

To obtain the time averaged MSD, we calculate the auxiliary integral

$$\begin{equation} I(\tau)=\frac{1}{T-\tau}\int_0^{T-\tau}\langle x(t)x(t+\tau)\rangle\mathrm{ d}t. \end{equation} \tag{ B.1 }$$

To this end, we use the representation of the hypergeometric function via the Fox H-function (see [47], equation (1.131)) and then employ the integral over the H-function following  [48], equation (8.3.2.7), arriving at

$$\begin{eqnarray} &&\fl I(\tau)=\frac{2^pp^{-2p}\Gamma(p+1)(D_{0}\tau)^{(p-1)/2}(D_{0}(T-\tau))^{ (p+1)/2}}{\Gamma(p/2+1/2)\Gamma(1/2-p/2)}\nonumber\\ &&\times H_{3,3}^{3,1}\left[\frac{\tau}{T- \tau}\left|\begin{array}{ll}(1,1),\left(\frac{3}{2},1\right),\left(\frac{5}{2} +\frac{p}{2},1\right)\\ \left(\frac{3}{2}+\frac{p}{2},1\right),\left(\frac{1}{2} -\frac{p}{2},1\right),\left(1+\frac{p}{2},1\right)\end{array}\right.\right]. \end{eqnarray} \tag{ B.2 }$$

Finally, we use  [48], equation (8.3.2.3)] to find the result for the time averaged MSD, equation (7).

We now outline an approximative scheme to evaluate the ergodicity breaking parameter EB. Namely, from equation (A.1) we express the particle displacement through |x(t)| ∼ |y(t)|^p, where y(t) is the Wiener process. Then, we take 〈x²(t)〉 = (2D₀tp⁻²)^p, which at Δ/T ≪ 1 yields $\langle \overline {\delta ^2(\Delta )}\rangle \approx (2D_{0}p^{-2})^p\Delta /T^{1-p}$. Thus $\mathcal {EB}(\Delta )=\langle \overline {\delta ^2(\Delta )}\rangle /\langle x^2(\Delta )\rangle \approx (\Delta /T)^{1-p}$, identical to the result in the main text derived from the exact approach (A.1). The fourth moment of the time averaged MSD in this approximative scheme becomes

$$\begin{eqnarray} &&\fl \left\langle\,(\overline{\delta^2(\Delta)})^2\,\right\rangle =\frac{1}{(T-\Delta)^2} \int_{0}^{T-\Delta}\mathrm{d}t'\int_{0}^{T-\Delta}\mathrm{d}t''\langle[x(t'+\Delta)-x(t')]^2[x(t''+\Delta)-x(t'')]^2 \rangle, \end{eqnarray} \tag{ C.1 }$$

involving a combination of altogether nine terms, each containing a product of four noise terms ζ(t). Using 〈ζ(t₁)ζ(t₂)ζ(t₃)ζ(t₄)〉 = δ(t₁ − t₂)δ(t₃ − t₄) + δ(t₁ − t₃)δ(t₂ − t₄) + δ(t₁ − t₄)δ(t₂ − t₃), the integrals from the Wiener processes disappear and two integrations are left. These were evaluated numerically for a given value of p to extract the scaling of the ergodicity breaking parameter

$$\begin{equation} \mathrm{EB}(\Delta)=\frac{\left\langle(\overline{\delta^2})^2\right\rangle - \left\langle\,\overline{\delta^2}\,\right\rangle ^2}{\left\langle \overline{\delta^2}\right\rangle^2} \end{equation} \tag{ C.2 }$$

for Δ/T ≪ 1. Using the result $\left\langle\,\overline {\delta ^2(\Delta )}\,\right\rangle ^2\propto \Delta ^2$, we expand the fourth moment up to Δ² to compute EB(Δ/T → 0), corresponding to the limit T → ∞. For Brownian motion (p = 1) we find the correct scaling EB = (4/3)(Δ/T). The scheme for HDPs with p = 2 yields EB(0) = 2/3, close to the simulations results. Overall, this approximative scheme predicts that EB grows as the system deviates from the ergodic state p = 1 for both sub- and superdiffusion, see figure 5.

We now determine the normalized velocity autocorrelation function

$$\begin{equation} C_v^{(\delta)}(\tau)=\frac{\langle[x(\tau+\delta)-x(\tau)][x(\delta)-x(0)] \rangle}{\sqrt{\langle[x(\tau+\delta)-x(\tau)]^2}\sqrt{\langle[x(\delta)- x(0)]^2\rangle}}. \end{equation} \tag{ D.1 }$$

With the position autocorrelations (4), C^(δ)_v(τ) becomes a universal function of δ/τ. C^(δ)_v(τ) is plotted in figure C.1. Thus, for superdiffusion the correlations remain positive, for subdiffusion C^(δ)_v(τ) features a negative region mirroring the anti-persistence of the motion. In figure C.1 we also show the complete decay of the correlations within the increment δ for Brownian motion. Moreover, we include the velocity autocorrelation function of FBM for the sub- and superdiffusive cases. The general behaviours are quite similar.

Zoom In Zoom Out Reset image size