Anomalous diffusion and ergodicity breaking in heterogeneous diffusion processes

We demonstrate the non-ergodicity of a simple Markovian stochastic processes with space-dependent diffusion coefficient $D(x)$. For power-law forms $D(x) \simeq|x|^{\alpha}$, this process yield anomalous diffusion of the form $\\simeq t^{2/(2-\alpha)}$. Interestingly, in both the sub- and superdiffusive regimes we observe weak ergodicity breaking: the scaling of the time averaged mean squared displacement $\{\delta^2}$ remains \emph{linear} and thus differs from the corresponding ensemble average $\$. We analyze the non-ergodic behavior of this process in terms of the ergodicity breaking parameters and the distribution of amplitude scatter of $\{\delta^2}$. This model represents an alternative approach to non-ergodic, anomalous diffusion that might be particularly relevant for diffusion in heterogeneous media.

Since the early systematic studies of Perrin and Nordlund [1], single particle tracking has become a routine method to trace individual particles in microscopic systems [2], but also of 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 2 (t) ≃ t p [5]. Subdiffusion (0 < p < 1) was found for the motion of tracers in living biological cells [6][7][8][9][10] or of bacteria in biofilms [11]. Superdiffusion (p > 1) is observed for active motion in biological cells [10], of mussels [12], plant lice [13], or higher animals and humans [3,4].
Of particular current interest are the ergodic properties of a stochastic process: is the information obtained from time averages of a single or few trajectories representative for the ensemble [14,15]? Of the various stochastic models for anomalous diffusion, diffusion on random fractals has been shown to be ergodic [16]. Fractional Brownian and fractional Langevin equation motion driven by long-range correlated Gaussian noise are transiently non-ergodic [17][18][19]. In contrast, subdiffusive Scher-Montroll continuous time random walks (CTRW) [20] with diverging characteristic waiting time exhibit weak ergodicity breaking [21], and time averages remain random even in the long time limit [22]. Weak ergodicity breaking was indeed observed for the blinking dynamics of quantum dots [23], protein motion in human cell walls [7], and lipid granule diffusion in yeast cells [8].
Can we come up with a simple stochastic model that still features such intricate non-ergodic behavior? Here we show that the MSD of heterogeneous diffusion processes (HDPs) with space dependent diffusion constant D(x) ≃ |x| α scales like x 2 (t) ≃ t p with p = 2/(2 − α), while the time averaged MSD δ 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 distribu-tion of δ 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 [24] as well as mesoscopic approaches to transport in heterogenous porous media [25] and on random fractals [26]. Recently, maps of the local cytoplasmic diffusivities in bacterial and eukaryotic cells showd a highly heterogeneous landscape [27], recalling the strongly time-varying diffusion coefficients of tracers in cells, see, e.g., Ref. [28].
Model. 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, dx/dt = 2D(x)ζ(t). 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 [29], to ensure the correct limiting transition for the noise with infinitely short correlation times [30]. Specifically, for D(x) we employ the power-law form [31] The dimension of D 0 is [D 0 ] = cm 2−α sec −1 . In the Stratonovich interpretation, we introduce the substitu- is the Wiener process [29]. For the initial condition P (x, 0) = δ(x) a compressed (α < 0) or stretched (α > 0) Gaussian PDF emerges (see Supplementary Material (SM) [32]), with p = 2/(2 − α). The ensemble averaged MSD x 2 (t) = x 2 P (x, t)dx becomes Thus, for α < 0 we find subdiffusion, and superdiffusion for α > 0. Brownian motion with x 2 (t) = 2D 0 t emerges for α = 0, and α = 1 produces ballistic motion. The diffusion becomes increasingly fast when α tends to 2. Fig. 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 Fig. 1 resembles the propagator for retarded wave motion [33].
To characterize the HDP, we calculate the position autocorrelation (see SM [32]) for the case t 2 > t 1 , with the hypergeometric function 2 F 1 (z). The Brownian limit x(t)x(t + τ ) = 2D 0 t follows for α = 0, while for t 1 = t 2 we recover Eq. (3). The asymptotic behavior of expression (4), implies that the correlations decay with τ for subdiffusive processes whereas they increase in the case of superdiffusion, similar to fractional Brownian motion (FBM) [17]. 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 p = 2, To connect to single particle tracking experiments we now turn to the time averaged MSD of a trajectory x(t), where ∆ is the lag time and T the length of the time series x(t). In Eq. (6) we introduced the additional average over individual trajectories, [22]. For ∆ ≪ T we find the linear ∆-dependence [32] Remarkably, this result holds for both sub-and superdiffusive regimes, and we find the exact match δ 2 (∆) = (∆/T ) 1−p x 2 (∆) . Eq. (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, δ 2 (∆) = x 2 (∆) . Linear scaling of the type of Eq. (7) was found for subdiffusive CTRWs [22], but also for correlated CTRWs [34]. It contrasts the ultraweak non-ergodicity of superdiffusive CTRWs [35]. 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 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 p = 2. To avoid divergencies in the discrete scheme, for subdiffusion we regularize the diffusivity at x = 0, choosing D(x) = D 0 A/(A + x 2 ), such that for x 2 ≫ 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 0 (|x| + 1), again with the correct scaling for large |x|.
Numerical results. From the generated trajectories, we compute the PDF (Fig. 1) as well as the ensemble and time averaged MSDs. The MSDs are shown in Fig. 2. For the ensemble averaged MSD (blue curve) we observe excellent agreement with Eq. (3). In the subdiffusion case, the predicted behavior (3) is reached after less than a dozen steps, such that the rectification of D(x) and the choice of the off-center initial position x 0 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 behavior shows good agreement between theory and simulations. The time averaged MSDs and their mean δ 2 follow nicely the predicted linear scaling. Only when the lag time ∆ approaches T the linear behavior levels off due to the finite trajectory length T . Next we study the apparent diffusion coefficients and scaling exponents from the individual δ 2 . The points of each trace δ 2 were logarithmically binned. The initial part of the trajectory, ∆ = 1 . . . 10 2 steps, was fitted by a linear law with diffusivity D 1 . The long-time regime, ∆ ∼ 10 2 . . . 10 4 steps, was fitted with two parameters, δ 2 = D β ∆ β . The distributions of D 1 , D β , and β were then fitted with the 3-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 δ 2 of a Brownian walker [36]. In addition, we extracted the amplitude scatter PDF φ(ξ = δ 2 / δ 2 ) of individual δ 2 around the mean δ 2 . Figs. 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 δ 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 [22]. For superdiffusive CTRWs φ(ξ) decays to zero at ξ = 0, but changes significantly with ∆   Figure 4: The same as in Fig. 3 but for superdiffusive HDPs. [35]. The scatter distribution φ(ξ) is thus a good criterion to distinguish HDPs and CTRWs. More specifically, for subdiffusion the distributions of D 1 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 3 . . .
that provides a good measure for the dispersion of time averaged MSDs for different types of diffusion processes [37]. 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, EB = δ 2 (∆) The simulations show that for subfdiffusion, EB increases with ∆/T , while for superdiffusion it decreases with ∆/T (Fig. 5). In both cases EB approaches unity as ∆ → T . At ∆/T ≪ 1, EB for p < 1 assumes smaller values, indicating weaker deviations from ergodicity, in contrast to superdiffusion, see Fig. 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).
Conclusions. The seemingly simple Markovian HDP with power-law space-dependent diffusion coefficient is 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 δ 2 (∆) scales linearly with the lag time ∆ for both sub-and superdiffusion. Such a behavior was previously known for subdiffusion from more complex processes: CTRWs with diverging characteristic waiting time [22] and CTRWs with correlated waiting times [34]. For superdiffusion only an ultraweak ergodicity breaking was reported in which the coefficients of the time and ensemble MSDs differ [35], while space-correlated CTRWs feature a quadratic scaling of the time averaged MSD [34].
The amplitudes of δ 2 (∆) 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 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 behavior 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. At the same time our observations for HDPs pose the interesting question for the minimal necessary condition for the occurrence of weak ergodicity breaking.
The authors thank E. Barkai and H. Buening for discussions and acknowledge funding from the Academy of Finland (FiDiPro scheme to RM) and Deutsche Forschungsgemeinschaft (Grant CH 707/5-1 to AGC).