# Non-ergodicity, fluctuations, and criticality in heterogeneous diffusion processes

###### Abstract

We study the stochastic behavior of heterogeneous diffusion processes with the power-law dependence of the generalized diffusion coefficient encompassing sub- and superdiffusive anomalous diffusion. Based on statistical measures such as the amplitude scatter of the time averaged mean squared displacement of individual realizations, the ergodicity breaking and non-Gaussianity parameters, as well as the probability density function we analyze the weakly non-ergodic character of the heterogeneous diffusion process and, particularly, the degree of irreproducibility of individual realization. As we show, the fluctuations between individual realizations increase with growing modulus of the scaling exponent. The fluctuations appear to diverge when the critical value is approached, while for even larger the fluctuations decrease, again. At criticality, the power-law behavior of the mean squared displacement changes to an exponentially fast growth, and the fluctuations of the time averaged mean squared displacement do not seem to converge for increasing number of realizations. From a systematic comparison we observe some striking similarities of the heterogeneous diffusion process with the familiar subdiffusive continuous time random walk process with power-law waiting time distribution and diverging characteristic waiting time.

###### pacs:

02.50.-r,05.40.-a,05.10.Gg,87.10.Mn## I Introduction

Over the last fifteen years there has been a surge of studies in anomalous diffusion, characterized by the deviation of the mean squared displacement (MSD)

(1) |

of a stochastic process with the probability density function to find the particle at position at time , from the linear time dependence of ordinary Brownian motion vankampen (). Anomalous diffusion is usually characterized in terms of the power-law form

(2) |

with the anomalous diffusion coefficient of physical dimension and the anomalous diffusion exponent . Depending on the value of we distinguish subdiffusion () and superdiffusion () report (); bouchaud ().

Since the milestone discoveries of superdiffusion in turbulence already in 1926 richardson () and of subdiffusion in amorphous semiconductors scher () the recent vast increase of interest in anomalous diffusion is due to its discovery in numerous microscopic systems, in particular in biological contexts. The cytoplasm of biological cells is heavily crowded with various obstacles, including proteins, nucleic acids, ribosomes, the cytoskeleton, as well as internal membranes compartmentalizing the cell crowd (); crowd2 (). Diffusion of natural and artificial tracers in this complex environment is often subdiffusive. Similar situations are encountered in cell membranes zhou (); weiss (). The experimental evidence for subdiffusion in the crowded cytoplasm of living cells ranges from the motion of small labeled proteins lang-subdiff-nucleus-fcs (); fradin05 () over mRNA molecules and chromosomal loci goldingcox (); weber (), lipid and insulin granules granules (); insulin (), virus particles brauch01 (); brauch02 (), as well as chromosomal telomeres telosubdiff () and Cajal bodies platani () inside the nucleus. Subdiffusion was observed for membrane resident proteins experimentally weiss (); channels () and for membrane lipid molecules in computer simulations ad-lipids (). In controlled in vitro experiments with artificial crowders, anomalous diffusion was consistently observed szymanski (); jae_njp (). On larger scales anomalous diffusion was observed, for instance, for the motion of bacteria in a biofilm biofilm ().

The observed subdiffusion was ascribed to various physical mechanisms saxton (); sokolov12-sm (); weiss-AD-SM (); pt (); franosch13 (). Apart from the apparent transient anomalous diffusion caused by a crossover from free normal diffusion to the plateau value of the MSD saxton (); berez (), typically, three main families of anomalous diffusion processes are considered: (i) diffusion in a fractal environment where dead ends and bottlenecks slow down the motion on all scales fractal (); (ii) motion in a viscoelastic environment, in which the effective anomalous motion of the tracer particle in the correlated many-body environment shows long-ranged antipersistent motion goychuk (). The latter process is associated with fractional Brownian motion (FBM) and generalized Langevin equation motion with a power-law memory form of the friction kernel fbm (); fle (). (iii) And continuous time random walk (CTRW) models, in which the moving particle is successively trapped by binding events to the environment or caging effects for waiting times distributed like a power-law with scher (); montroll (). All three mechanisms lead to the power-law MSD (2) and they were indeed identified as processes generating the motion of different tracers in different cellular environments weber (); granules (); insulin (); channels (); ad-lipids (); szymanski (); jae_njp (); sokolov12-sm (); weiss-AD-SM (); pt (); franosch13 () or in colloidal systems in vitro xu ().

These three anomalous diffusion mechanisms characterized by anomalous diffusion with a constant in time generalized diffusivity were identified in experiments involving fairly large endogenous as well as artificial tracers. Recent experiments on eukaryotic cells lang11 () using considerably smaller tracer proteins sampling over much larger subvolumes of the cell indicated a systematic variation of the cytoplasmic diffusivity with the separation from the cell nucleus. These spatial diffusivity gradients are due partly to the non-uniform distribution of crowders in the cytoplasm. This distribution can non-trivially affect the diffusion of tracers of different sizes in the cell cytoplasm. In vitro, fast gradients of the diffusivity can be realized, for instance, via a local variation of the temperature in thermophoresis experiments therm1 (); therm2 (), as sketched in Fig. 1. The diffusion of Brownian particles in explicit solvents with temperature gradients was recently studied by multi-particle collision dynamics ripoll-12 (). Fluctuation-dissipation relations and spurious drift effects in systems with spatially varying friction coefficient were recently treated theoretically in Ref. farago14-friction (). On larger scales, the diffusion of water molecules monitored by diffusive Magnetic Resonance Imaging in the brain white matter was demonstrated to be heterogeneous and strongly anisotropic brain-diffusion (). The anisotropy is due to the presence of some spatially-oriented structures in the tissue and obstacles which give rise to a tensorial character of the apparent diffusion coefficient. The existence of a population splitting into two pools of water molecules with slow and fast diffusivities was shown for brain white matter brain-diffusion (). Finally, spatial heterogeneities are also abundant in the completely different context of anomalous diffusion in subsurface hydrology hydro ().

In the present study we examine the effects of the strength of the diffusivity gradient as defined by the scaling exponent and of the initial particle position in heterogeneous diffusion processes (HDPs) with space-dependent diffusion coefficient . We pay particular attention to a phenomenon, that recently received considerable attention for its immediate relevance to the surging field of single particle tracking experiments in microscopic systems, namely, the so-called weak ergodicity breaking. This is the distinct disparity between physical observables depending on whether they are evaluated in the conventional ensemble sense or, from measured time series of the particle position, in terms of time averages russians-eb (); bouchaud-eb (); pt (); sokolov12-sm (); franosch13 (). We find that the HDP process gives rise to weakly non-ergodic behavior with a pronounced amplitude scatter of the time averaged MSD of individual trajectories. We obtain details of the distribution of this scatter as well as the frequently used ergodicity breaking and non-Gaussianity parameters. Our analysis uncovers remarkable similarities of the HDP process with those of CTRW motion. In particular, we study the behavior of the HDP at the critical value of the scaling exponent of the diffusivity.

The paper is organized as follows. We introduce the HDP model with -dependent diffusivity in Sec. II. The main results for the evolution of the MSD, the time averaged MSD, the probability density function, the ergodicity breaking and non-Gaussianity parameters in the whole range of the model parameters are then presented in Sec. III. We discuss our results and point out the directions for future research in Sec. IV.

## Ii Heterogeneous diffusion processes

HDPs are defined in terms of the multiplicative yet Markovian Langevin equation hdp ()

(3) |

where is the position-dependent diffusion coefficient and represents white Gaussian noise. In what follows we concentrate on the power-law form

(4) |

for the diffusivity, where the amplitude has dimension . Logarithmic and exponential forms for were considered in Ref. hdp (). The offset in Eq. (4) avoids either divergencies of () or stalling of the particle () around in the simulations. In the following calculations we use the scaling form . We interpret the Langevin equation (3) in the Stratonovich sense vankampen (); hdp ().

The MSD following from the stochastic equation (3) with diffusivity (4) takes on the power-law form hdp ()

(5) |

where we introduced the scaling exponent

(6) |

which denotes superdiffusion for and subdiffusion for . For the diffusivity grows away from the origin, leading to a progressive acceleration of the particle as it ventures into more distant regions from the origin, and vice versa for . In the special case , the theory developed in Ref. hdp () breaks down as the scaling exponent (6) diverges. In Refs. fulinski (); lubensky07 () an exponential growth of the MSD was found, see the discussion below. For even larger values of the anomalous diffusion exponent becomes negative, i.e., we observe a strong localization, see the discussion below. The probability density function (PDF) of the HDP is given by the stretched or compressed Gaussian hdp ()

(7) |

For positive it exhibits a cusp, while for negative the PDF features a dip to zero at the origin.

In single particle tracking experiments one measures the time series of the particle position for a time span . It is usually evaluated in terms of the time averaged MSD pt (); sokolov12-sm (); franosch13 ()

(8) |

where is the lag time. Brownian motion is ergodic and for sufficiently long measurement times the equality holds pccp (); he (). Anomalous diffusion (2) described by FBM or fractional Langevin equation motion is asymptotically ergodic deng (); goychuk () but may feature transiently non-ergodic behavior jae_pre (); jae_njp () in confinement as well as transient aging jochen (), the explicit dependence on the time span elapsing between initial preparation of the system and start of the recording of the position, .

CTRW processes with diverging time scales of the waiting time distribution show weakly non-ergodic behavior for all . Namely, despite the scaling (2) of the MSD, scales linearly with . More precisely, if we average over sufficiently many trajectories, the quantity

(9) |

for CTRWs scales like lubelski (); he (); pccp (); neusius (). This linear dependence on is preserved for aging CTRWs johannes (). Apart from this linear -dependence we also observe the dependence on the process time , a signature of aging pt ().

Interestingly, the HDP with diffusivity (4) displays weak ergodicity breaking of the form hdp (); fulinski (); spain ()

(10) |

This behavior is analogous to that of scale-free CTRW motion, despite of the fact that the increment correlation function of HDPs are (anti)persistent in analogy to the ergodic FBM hdp (); deng (). We note that also other processes such as aging and correlated CTRWs vincent (); lomholt () as well as scaled Brownian motion with time-dependent diffusivity novikov (); sbm-eb () exhibit the duality between (2) and a linear -dependence of . Eq. (10) shows the -scaling as function of the process time . For subdiffusion with , that is, the effective diffusivity of the process decays over time, as the particle ventures into low-diffusivity areas. Conversely, for the diffusivity increases over time, the particle discovers areas with increasingly higher .

In the current paper we perform a detailed analysis of the MSD and the time averaged MSD in the entire range of , including the critical value . We are particularly interested in the extent of the weakly non-ergodic behavior, especially the fluctuations of the time averaged MSD around the mean value characterized by the ergodicity breaking parameter. Moreover we analyze the non-Gaussianity of the process. In our analysis we study effects of the scaling exponent of as well as the initial position of the particle. The latter is known to affect the time-scales at which anomalous diffusion becomes significant hdp (); kehr-87 (). We combine analytical and numerical approaches. The simulations scheme for HDPs was introduced in Ref. hdp ().

## Iii Results

In this Section we start with the analysis of the MSD and the time averaged MSD with its amplitude fluctuations, the latter being quantified by the corresponding scatter distribution. We then analyze the ergodicity breaking and non-Gaussianity parameters. Finally, we study the PDF of the function that quantifies the degree of particle dispersion.

### iii.1 MSD and time averages MSD

Fig. 2 shows the results from computer simulations for the MSD and the time averaged MSD from individual realizations along with the average taken over all single time traces. The values for the scaling exponent of the diffusivity (4) studied in Fig. 2 cover both the sub- and superdiffusive domains and include, in particular, the critical value where the scaling exponent of the MSD (2) diverges. For each we show single trajectories. In all cases, apart from the critical point, the scaling of the MSD and both individual () and mean () time averaged MSDs agree well with the expected analytical behavior shown by the dashed lines: the scaling exponent of the MSD varies consistently with , while that of the time averaged MSD remains unity throughout. The initial discrepancy between theory and MSDs is due to the choice for the initial position , whose influence relaxes on a time scale depending on both and . The deviations of individual time traces at long lag times from the predicted behavior is due to unavoidable, bad statistics when the lag time gets close to the overall length of the time series.

Irreproducibility of time averages of physical observables such as the MSD is an intrinsic property of weakly non-ergodic processes pt (); he (); pccp (); lubelski (). For CTRW processes with scale-free waiting time distribution individual traces contain different, few dominating waiting time events that cause the amplitude scatter between different trajectories or, in other words, fluctuations of the apparent effective diffusion constant. This phenomenon occurs no matter how long the measurement time is taken. For HDPs the scatter is due to the difference in the extent of excursions to regions of significantly different diffusivity.

The amplitude scatter between individual realizations in our HDP shown in Fig. 2 varies significantly with the value of the scaling exponent : away from the Brownian value in both subdiffusive () and superdiffusive () cases the fluctuations of become more pronounced when the modulus of increases. At the critical point the fluctuations of appear to diverge, while beyond this critical point, the fluctuations decrease again. Moreover, a population splitting in a faster (steeper slope of ) and slower (shallower slope) fraction of trajectories appears hdp (), see especially the panel for . In terms of the dimensionless variable

(11) |

the amplitude scatter distribution reflects the randomness of individual time averages of the MSD. For a sub- and a superdiffusive it was analyzed in Ref. hdp (). As shown here for a whole spectrum of values, there is a clear trend towards extreme fluctuations at the critical point , but even for considerably smaller values such as the fluctuations around the mean are significant. The width of the fluctuations of in each panel varies only marginally with the lag time , apart from the behavior at . As studied in Ref. hdp (), the relative amplitude scatter distribution can be fitted with the generalized Gamma-distribution. In particular, it tends to zero at , in contrast to subdiffusive CTRW processes, for which is always positive, indicating completely stalled trajectories he (); jae_jpa (); pt (); pccp ().

We note that when the exponent approaches the critical value , the number of steps necessary to approach the theoretically predicted asymptote (5) increases significantly. Thus, as seen from Fig. 2, for only few simulation steps, are needed for negative with larger modulus. It increases to steps for , for , and already for .

Once the anomalous diffusion exponent becomes negative, that is, for values of larger than the critical value , the MSD (5) becomes a decreasing function of time. This agrees with our simulations if is chosen sufficiently large to enable the relaxation to the theoretical asymptote. For instance, for the extreme value and values of the initial position of and above the MSD indeed follows the theoretical prediction (not shown). In this region the diffusivity grows fast away from the origin the decreasing MSD corresponds to the localization of particles in regions of slow diffusivity, see the detailed discussion in Sec. III.4 below.

In the limit , due to the huge spread one or few extremely large amplitudes of the time averaged MSD may substantially affect the mean . The MSD in this limit follows an exponential growth , as indicated by the dashed curve in the panel for in Fig. 2. Such a fast increase of the MSD is consistent with the divergence of the scaling exponent as function of in Eq. (5) as well as with the exponential MSD growth predicted for a parabolically space-varying diffusivity in Refs. lubensky07 (); fulinski (). This property can be straightforwardly inferred from the diffusion equation,

(12) |

Multiplying both sides with and integrating over one arrives at lubensky07 ()

(13) |

This also rationalizes the observation that the scatter of is maximal for as the particles perform extremely far-reaching excursions relative to other -values.

### iii.2 Ergodicity breaking parameters

As introduced in Refs. russians-eb (); he () the fluctuations of the time averaged MSD can be quantified by their variance, the ergodicity breaking parameter

(14) |

When , it means that the process is perfectly reproducible and all time averages over sufficiently long trajectories give the same value. This case corresponds to the sharp scatter distribution pt (); he (); pccp (); jae_jpa (). For the canonical Brownian Motion the ergodicity breaking parameter reaches the zero value in the form at finite ratio . Weakly non-ergodic processes have a positive value of . An alternative, weaker condition for ergodicity is when the MSD (2) and the time averaged MSD (9) coincide. To measure the relative deviations from ergodicity, the parameter was introduced levywalks1 ()

(15) |

For ergodic dynamics its value is unity.

To calculate the parameter analytically is not always an easy task. To see this, note that the MSD follows from the change of the stochastic variable to the standard Wiener process in the form hdp (). For the time averaged MSD the calculation is already more complicated, as it involves the two-point position correlation function. The latter is expressed via Fox -functions for the HDP hdp (). The analytical derivation of the ergodicity breaking parameter , however, requires fourth-order moments, whose calculation is a formidable task. So far only approximate methods are known for HDPs hdp (). For subdiffusive CTRW processes it is possible to obtain the parameter more easily from the conjectured and numerically proven equivalence of with the ratio of the number of steps in an individual realization and the average he (); johannes (). As for the limiting distribution for is known hughes (), this allows straightforward calculation of and its moments he (); johannes (). For the multiplicative process studied here such a scheme does not work. We also note that another parameter involving the fourth moment of the particle displacement is the non-Gaussianity parameter franosch13 () discussed below. Due to the lack of an analytical theory, a major reason for the current simulations study is to explore the behavior of these two parameters in the whole range of the model parameters.

Figs. 3 and 4 display the dependence of the ergodicity breaking parameters on the lag time , the initial position , as well as the scaling exponent . We find that for standard Brownian motion with the ergodicity breaking parameter EB follows the known asymptote from above. Concurrently as expected for ergodic motion, except for very small values because of the initial relaxation of the influence of the initial position . As we depart from the value of Brownian Motion, the magnitude of grows and its functional dependence on the lag time becomes less pronounced, see Fig. 3a. Approaching the critical point , due to the huge fluctuations of and the ensuing possibility of extreme events the ergodicity breaking parameter also explicitly depends on the number of traces used for the averaging and reaches values of and higher. On both sides of the critical point, corresponding to large positive or negative values of , the values of the parameter approach one another, as seen in Fig. 4a, while exhibits a jump, Fig. 4b.

The ratio of the different MSDs given by also depends on the power exponent and the initial value . As anticipated already from Fig. 2, for the chosen initial condition the magnitude of decreases as gets progressively negative, while grows as increases from 0 to 2, reaching very large values at . This is illustrated in Fig. 3b which is also consistent with the scaling

(16) |

with the lag time at different values. Because of the initial condition , similar to our statements for this asymptote is approached later when .

In Fig. 4 we also analyze the effects of the initial position and show the dependence of the ergodicity breaking parameters for long traces or short lag times, i.e., when which is practically equivalent to , versus the anomalous diffusion exponent and the scaling exponent . We observe that in the region the values of obtained from simulations is in good agreement with the analytical estimate from Ref. hdp (), represented by the dashed curve in Fig. 4a and b. These approximations are based on the asymptotic scaling of and do not consider the effect of the initial positions , assumed to be relaxed in the relevant limit. To see their actual impact on the dynamics we study the ergodicity breaking parameters numerically.

We observe that the value of grows from the small pre-asymptotic Brownian value at to progressively larger values at larger modulus of , as approaches the critical value from below and above. The functional dependence of obtained from our simulations follows quite well the predictions from the approximate calculation of hdp (), particularly in the range , see Fig. 4b. The deviations for even more pronounced variation of as tends to 2 are likely due to insufficient statistics. As we show in Fig. 4a by the circles and triangles, the value reveals measurable deviations for as compared to traces used for the averaging, while for less extreme values the two sets yield nearly identical results. This is a further indication towards a divergence of the fluctuations at the critical point. Namely, the chance to find an even larger amplitude of in a given realization grows with the number of simulated trajectories.

As already anticipated in Ref. hdp () the HDP process considered here is in some aspects reminiscent of subdiffusive CTRW processes with scale-free, power-law waiting time distributions with regard to the scaling of the MSD and the time averaged MSD, compare also Sec. IV. Apparently, the ergodicity breaking properties of these two processes also reveal similar features: the fluctuations of increase when the anomalous diffusion exponent becomes successively smaller than the Brownian value 1. Indeed, as shown in Fig. 4b for the relevant region there exists even a close agreement of the quantitative behavior of CTRW and HDP. We show the ergodicity breaking parameter obtained for the CTRW process he ()

(17) |

where is the exponent in the PDF of the waiting times report (). For the CTRW process the exponent also occurs in the MSD (2). We note that in contrast to HPDs in superdiffusive CTRW processes the nonergodic behavior is merely ultraweak, effecting a different prefactor in the time averaged MSD compared to the MSD levywalks (); levywalks1 ().

For we find that for initial positions close to the origin, in agreement with the results presented in Ref. hdp (). For or the analytical model of Ref. hdp () no longer applies but the simulations yield another region of growth for the value of . Note that also depends on the trace length (not shown), as discussed for two-dimensional HDPs hdp (). As we show in Fig. 4c, the magnitude of strongly varies with . In the region the ratio of the MSDs grows as indicated by the dashed asymptote in Fig. 4c, while for large negative we observe . We finally note that reveals a much weaker dependence on the number of simulated traces, see the circles and triangles in Fig. 4c. The reason is that involves only the second moment of the time averaged MSD, while the is defined in terms of the fourth order moment, which is more sensitive to large variations.

### iii.3 Non-Gaussianity parameter

The non-Gaussianity parameter is a sensitive experimental indicator that often enables one to distinguish the type of particular diffusion processes observed in single-particle tracking experiments weiss14-pccp (). It is related to the stationarity of increments of the diffusion process and involves the fourth moment along the time-averaged trajectory. For the diffusion process in an embedding space of dimension it is defined via the experimentally relevant time averaged quantities as franosch13 ()

(18) |

where, in analogy to Eq. (8), the fourth moment is defined via

(19) |

For Brownian and fractional Brownian motion, both Gaussian processes, one consistently finds . For diffusion processes revealing a transient anomalous diffusion behavior, this parameter deviates substantially from zero, see, e.g., the discussion in Ref. weiss14-pccp ().

We systematically examine the behavior of in Fig. 5a. Corroborating the results for the ergodicity breaking parameter, for Brownian motion we obtain approximately zero values. As the power in deviates from zero, the non-Gaussianity parameter reveals a rich behavior as a function of the lag time and the initial position .

More specifically, in the region the non-Gaussianity parameter progressively grows and reaches considerably large values for large , see the curve for in Fig. 5a. The value systematically decreases with the lag time along the trace for negative . In the region the value of also grows with but, in contrast to the case , the function stays rather constant with . As we approach the critical the non-Gaussianity parameter reaches very high values. In the region the value of decreases again. These features of the functional behavior of correspond with the properties of , compare the curves in Fig. 5a and 3a. Note that similar to the statistics required to obtain a smooth curve for strongly depends on : the result from traces in Fig. 5a acquires pronouncedly higher fluctuations when we approach the critical value .

We find that variations in the initial positions only have a marginal effect on . In Fig. 5b we illustrate the range of values for exponents in the range . Even for the wide range from to the value of does not reveal any appreciable variation, in stark contrast to the strong -dependence of . This indicates that for the HDP the non-Gaussianity parameter is more robust than with respect to the choice of the initial condition. In Fig. 5b the dashed line indicates the proportionality to , which nicely matches the measured shape of .

### iii.4 Probability density function and particle focusing

We finally analyze the time evolution of the PDF of the particle position in Fig. 6 for three different values of the scaling exponent , distinguishing superdiffusion (top row for , i.e., ballistic motion with ), normal diffusion (middle row for and ), and subdiffusion (bottom row for , i.e., ). In each row the leftmost panel represents a very early evolution of the PDF close to the initial condition, while in the rightmost panel the PDF in some of the cases has almost reached the diffusion limit, in which the analytical asymptote (7) is valid.

We observe that for HDPs with positive , i.e., when the diffusivity grows with the modulus of the position, the PDF for particles with an initial position far away from the origin, an asymmetric shape of the PDF is effected. Namely, they progressively move towards regions of small diffusivity and accumulate there. This focusing due to the quenched nature of the diffusivity erases any memory of the initial condition and the common asymptote (7) of stretched Gaussian shape () is approached. In the opposite case with negative the focusing of particles in lower diffusivity regions applies again, albeit at longer times than shown here. This time, however, instead of the cusp at the origin for positive , the PDF drops down to zero at the origin and acquires the bimodal shape predicted by (7), a compressed Gaussian. For the Brownian case with vanishing , no focusing takes place. In the final panel the distribution is already so broad that the individual, shifted PDFs appear on top of each other.

Similarly to the evolution of the MSD shown in Fig. 2, the relaxation time required for the system of diffusing particles to approach this long-time limit, depends on the power of the diffusivity in Eq. (4). For more negative values the traces are not yet relaxed to the theoretical long-time shape even after steps. The particles starting at larger remain trapped for the whole length of the simulation. Conversely, when the diffusivity increases fast away from the origin (moderate positive values), the equilibration to the long-time PDF is relatively fast, see the panel for in Fig. 6.

## Iv Discussion and Conclusions

For HPD processes whose diffusivity varies in power-law form with the distance from the origin and which give rise to anomalous diffusion we analyzed in detail the weakly non-ergodic behavior for varying power exponents and initial positions of the particles. In particular, we examined the functional dependencies and magnitudes of the variation of the ensemble and time averaged characteristics of such HDPs with these parameters. The fluctuations of the amplitude of the time averaged MSD of individual time traces were shown to systematically grow with increasing departure of the anomalous diffusion exponent from the Brownian value , corresponding to the scaling exponent in the power-law form for . The fluctuations increase dramatically when the scaling exponent approaches its critical value . At that critical point we corroborated the turnover from the power-law scaling in time of the MSD to an exponential growth. Beyond the critical point the fluctuations of the time averaged MSD become smaller again. Concurrently we observe a distinct population splitting into a faster and slower fraction of time averaged MSDs.

We paid particular emphasis on several parameters used to classify the departure from ergodic behavior, namely the two ergodicity breaking parameters and as well as the non-Gaussianity parameter . While is simply defined as the ratio of the mean time averaged MSD versus the MSD, both and are based on the fourth order moments of the particle position and are known analytically only from approximate theories. A detailed numerical analysis of their properties was therefore used to obtain more concrete information on their behavior for different values of the scaling exponent and the initial particle position in the heterogeneous environment REM (). Within the analyzed parameter range we find that the behavior of is indeed in agreement with the heuristic theoretical analysis from Ref. hdp (). The parameter portraits for the ergodicity breaking and non-Gaussianity parameters obtained from our numerical analysis will be useful for actual data evaluation. We also explored the detailed behavior of the HDP dynamics at the critical value and its vicinity, in particular with respect to the dramatic values reached by reflecting the dramatic fluctuations of the time averaged MSD. A systematic numerical analysis of the particle PDF for varying initial conditions sheds additional light on HDPs with different exponent , in particular, the visualization of the particle focusing in low diffusivity regions.

In HDPs the non-ergodic behavior arises due to the heterogeneity of the environment. Physically, this represents a space-dependent mobility or temperature. It is not a property of the particle, and each time the particle revisits a given point in space it has the same diffusivity . In a random walk sense, this scenario could also be translated into a local dependence of the waiting time for a jump event. In that interpretation the HDP corresponds to a motion in a quenched energy landscape bouchaud (), albeit a deterministic (in contrast to random) one. As such, the process is intrinsically different from renewal CTRW processes, which correspond to the motion in an annealed environment bouchaud ().

Nevertheless we observed that subdiffusive HDPs share a number of features with subdiffusive CTRWs with scale-free, power-law waiting time distribution. These features include the scaling laws for the ensemble and time averaged MSDs, and the form of the ergodicity breaking parameter (in the relevant range ). From the sole analysis of the MSDs and as well as the ergodicity breaking parameters, a significant distinction between CTRW and HDP is therefore difficult. Yet there exist some crucial differences between the HDP and CTRW processes, that can be measured. Thus, in HDPs with the distribution of the relative amplitude of individual realizations decays to zero at , in contrast to the CTRW’s finite fraction of immobile particles reflected in the finite value . Indeed the scatter distribution was previously advocated as a good diagnosis tool for different stochastic processes pt (); jae_jpa (), complementing other stochastic analysis methods methods (); weiss14-pccp (); thiel ().

HDPs are physically meaningful alternatives to the standard stochastic models for anomalous diffusion processes, namely, the CTRW process with long-tailed waiting time distribution, FBM and fractional Langevin equation motion based on Gaussian yet long-ranged in time correlations, diffusion in fractal environments, or their combinations granules (); insulin (); channels (); thiel (). HPDs are weakly non-ergodic, sharing some yet not all properties of CTRW processes. It is therefore important to explore the properties of HDPs in more detail in the future as well as to develop more sophisticated tools to interpret measured time series and with confidence identify the underlying stochastic mechanism. In particular, we will study the properties of confined and aged HDPs as well as spatially varying diffusivities with a random component, and the coupling of HDPs with active motion.

###### Acknowledgements.

The authors acknowledge funding from the Academy of Finland (FiDiPro scheme to RM) and the German Research Council (DFG Grant CH 707/5-1 to AGC).## References

- (1) N. G. van Kampen, Stochastic processes in physics and chemistry (North Holland, Amsterdam, NL, 1981).
- (2) J.-P. Bouchaud and A. Georges, Phys. Rep. 195, 127 (1990).
- (3) R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
- (4) L. F. Richardson, Proc. Roy. Soc. London Ser. A 110, 709 (1926).
- (5) H. Scher and E. W. Montroll, Phys. Rev. B 12, 2455 (1975).
- (6) S. B. Zimmerman and A. P. Minton, Annu. Rev. Biophys. Biomol. Struct. 22, 27 (1993).
- (7) H. X. Zhou, G. Rivas, and A. P. Minton, Annual Rev. Biophys. 37, 375 (2008).
- (8) H.-X. Zhou, J. Phys. Chem. B 113, 7995 (2009).
- (9) M. Weiss, M. Elsner, F. Kartberg, and T. Nilsson, Biophys. J. 87, 3518 (2004).
- (10) M. Wachsmuth, W. Waldeck, and J. Langowski, J. Mol. Biol. 298, 677 (2000).
- (11) D. S. Banks and C. Fradin, Biophys. J. 89, 2960 (2005).
- (12) I. Golding and E. C. Cox, Phys. Rev. Lett. 96, 098102 (2006).
- (13) S. C. Weber, A. J. Spakowitz, and J. A. Theriot, Phys. Rev. Lett. 104, 238102 (2010).
- (14) J.-H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede, and R. Metzler, Phys. Rev. Lett. 106, 048103 (2011).
- (15) S. M. A. Tabei, S. Burov, H. Y. Kim, A. Kuznetsov, T. Huynh, J. Jureller, L. H. Philipson, A. R. Dinner, and N. F. Scherer, Proc. Natl. Acad. Sci. USA 110, 4911 (2013).
- (16) G. Seisenberger, M. U. Ried, T. Endreß, H. Büning, M. Hallek, and C. Bräuchle, Science 294, 1929 (2001).
- (17) C. Bräuchle, G. Seisenberger, T. Endreß, M. U. Ried, H. Büning, and M. Hallek, Chem. Phys. Chem. 3, 299 (2002).
- (18) I. Bronstein, Y. Israel, E. Kepten, S. Mai, Y. Shav-Tal, E. Barkai, and Y. Garini, Phys. Rev. Lett. 103, 018102 (2009); K. Burnecki, E. Kepten, J. Janczura, I. Bronshtein, Y. Garini, and A. Weron, Biophys. J. 103, 1839 (2012).
- (19) M. Platani, I. Goldberg, A. I. Lamond, and J. R. Swedlow, Nat. Cell Biol. 4, 502 (2002).
- (20) A. V. Weigel, B. Simon, M. M. Tamkun, and D. Krapf, Proc. Natl. Acad. Sci. USA 108, 6438 (2011).
- (21) G. R. Kneller, K. Baczynski, and M. Pasenkiewicz-Gierula, J. Chem. Phys. 135, 141105 (2011); J.-H. Jeon, H. Martinez-Seara Monne, M. Javanainen, and R. Metzler, Phys. Rev. Lett. 109, 188103 (2012); M. Javanainen, H. Hammaren, L. Monticelly, J.-H. Jeon, M. S. Miettinen, H. Martinez-Seara Monne, R. Metzler, and I. Vattulainen, Faraday Disc. 161, 397 (2013).
- (22) J. Szymanski and M. Weiss, Phys. Rev. Lett. 103, 038102 (2009); W. Pan, L. Filobelo, O. Galkin, V. V. Uzunova, and P. G. Vekilov, ibid. 102, 058101 (2009).
- (23) J.-H. Jeon, N. Leijnse, L. B. Oddershede, and R. Metzler, New J. Phys. 15, 045011 (2013).
- (24) S. S. Rogers, C. van der Walle, and T. A. Waigh, Langmuir 24, 13549 (2008).
- (25) M. J. Saxton and K. Jacobson, 26, 373 (1997).
- (26) D. Ernst, M. Hellmann, J. Kohler, and M. Weiss, Soft Matter 8, 4886 (2012).
- (27) I. M. Sokolov, Soft Matter 8, 9043 (2012).
- (28) F. Höfling and T. Franosch, Rep. Progr. Phys. 76, 046602 (2013).
- (29) E. Barkai, Y. Garini, and R. Metzler, Phys. Today 65(8) (8), 29 (2012).
- (30) A. M. Berezhkovsky, L. Dagdug, and S. M. Bezrukov, Biophys. J. 106, L09 (2014).
- (31) S. Havlin and D. Ben-Avraham, Adv. Phys. 36, 695 (1987); A. Klemm, R. Metzler, and R. Kimmich, Phys. Rev. E 65, 021112 (2002).
- (32) I. Goychuk, Phys. Rev. E 80, 046125 (2009); Adv. Chem. Phys. 150, 187 (2012).
- (33) A. N. Kolmogorov, Dokl. Acad. Sci. USSR 26, 115 (1940); B. B. Mandelbrot and J. W. van Ness, SIAM Rev. 1, 422 (1968).
- (34) S. C. Kou, Ann. Appl. Stat. 2, 501 (2008); E. Lutz, Phys. Rev. E 64, 051106 (2001); R. Kupferman, J. Stat. Phys. 114, 291 (2004); W. Min, G. Luo, B. J. Cherayil, S. C. Kou, and X. S. Xie, Phys. Rev. Lett. 94, 198302 (2005).
- (35) E. W. Montroll and G. H. Weiss, J. Math. Phys. 6, 167 (1965).
- (36) Q. Xu, L. Feng, R. Sha, N. C. Seeman, and P. M. Chaikin, Phys. Rev. Lett. 106, 228102 (2011); I. Y. Wong, M. L. Gardel, D. R. Reichman, E. R. Weeks, M. T. Valentine, A. R. Bausch, and D. A. Weitz, ibid. 92, 178101 (2004).
- (37) T. Kühn, T. O. Ihalainen, J. Hyväluoma, N. Dross, S. F. Willman, J. Langowski, M. Vihinen-Ranta, and J. Timonen, PLoS ONE 6, e22962 (2011).
- (38) C. B. Mast, S. Schink, U. Gerland, and D. Braun, Proc. Natl. Acad. Sci. USA 110, 8030 (2013).
- (39) Y. T. Maeda, T. Tlusty, and A. Libchaber, Proc. Natl. Acad. Sci. USA 109, 17972 (2012).
- (40) M. Yang and M. Ripoll, J. Chem. Phys. 136, 204508 (2012).
- (41) O. Farago and N. Grønbech-Jensen, E-print arXiv:1402.4598.
- (42) D. le Bihan, Nature Rev. Neurosci. 4, 469 (2003).
- (43) M. Dentz and D. Bolster, Phys. Rev. Lett. 105, 244301 (2010); B. Berkowitz, A. Cortis, M. Dentz, and H. Scher, Rev. Geophys. 44, RG2003 (2006).
- (44) S. M. Rytov, Yu. A. Kravtsov, and V. I. Tatarskii, Principles of Statistical Radiophysics 1: Elements of Random Process Theory (Springer, Heidelberg 1987).
- (45) J.-P. Bouchaud, J. Phys. I (Paris) 2, 1705 (1992); G. Bel and E. Barkai, Phys. Rev. Lett. 94, 240602 (2005); A. Rebenshtok and E. Barkai, ibid. 99, 210601 (2007); M. A. Lomholt, I. M. Zaid, and R. Metzler, Phys. Rev. Lett. 98, 200603 (2007).
- (46) A. G. Cherstvy, A. V. Chechkin, and R. Metzler, New J. Phys. 15, 083039 (2013); Soft Matter 10, 1591 (2014); A. G. Cherstvy and R. Metzler, Phys. Chem. Chem. Phys. 15, 20220 (2013).
- (47) A. W. C. Lau and T. C. Lubensky, Phys. Rev. E 76, 011123 (2007).
- (48) A. Fuliński, J. Chem. Phys. 138, 021101 (2013); Phys. Rev. E 83, 061140 (2011); Acta Phys. Polonica 44, 1137 (2013).
- (49) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
- (50) S. Burov, J.-H. Jeon, R. Metzler, and E. Barkai, Phys. Chem. Chem. Phys. 13, 1800 (2011); S. Burov, R. Metzler, and E. Barkai, Proc. Natl. Acad. Sci. USA 107, 13228 (2010).
- (51) W. Deng and E. Barkai, Phys. Rev. E. 79, 011112 (2009).
- (52) J.-H. Jeon and R. Metzler, Phys. Rev. E 85, 021147 (2012).
- (53) J. Kursawe, J. H. P. Schulz, and R. Metzler, Phys. Rev. E 88, 062124 (2013).
- (54) A. Lubelski, I. M. Sokolov, and J. Klafter, Phys. Rev. Lett. 100, 250602 (2008).
- (55) T. Neusius, I. M. Sokolov, and J. C. Smith, Phys. Rev. E 80, 011109 (2009).
- (56) J. H. P. Schulz, E. Barkai, and R. Metzler, Phys. Rev. Lett. 110, 020602 (2013); Phys. Rev. X 4, 011028 (2014).
- (57) Compare also the discussion in P. Massignan, C. Manzo, J. A. Torreno-Pina, M. F. García-Parajo, M. Lewenstein, G. J. Lapeyre Jr., E-print arXiv:1401.6110.
- (58) V. Tejedor and R. Metzler, J. Phys. A 43, 082002 (2010); M. Magdziarz, R. Metzler, W. Szczotka, and P. Zebrowski, Phys. Rev. E 85, 051103 (2012). See also the alternative approach in A. V. Chechkin, M. Hofmann, and I. M. Sokolov, Phys. Rev. E 80, 031112 (2009); J. H. P. Schulz, A. V. Chechkin, and R. Metzler, J. Phys. A. 46, 475001 (2013).
- (59) M. A. Lomholt, L. Lizana, R. Metzler, and T. Ambjörnsson, Phys. Rev. Lett. 110, 208301 (2013).
- (60) F. Thiel and I. M. Sokolov, Phys. Rev. E 89, 012115 (2014).
- (61) D. S. Novikov, E. Fieremans, J. H. Jensen, and J. A. Helpern, Nature Phys. 7, 508 (2011).
- (62) J. W. Haus and K. W. Kehr, Phys. Reports 150, 263 (1987).
- (63) J.-H. Jeon and R. Metzler, J. Phys. A 43, 252001 (2010).
- (64) A. Godec and R. Metzler, Phys. Rev. Lett. 110, 020603 (2013).
- (65) B. D. Hughes, Random Walks and Random Environments, Vol. I: Random Walks (Oxford University Press, Oxford, UK, 1995).
- (66) G. Zumofen and J. Klafter, Physica D 69, 436 (1993); A. Godec and R. Metzler, Phys. Rev. E 88, 012116 (2013); D. Froemberg and E. Barkai, ibid. 87, 030104(R) (2013); ibid. 88, 024101 (2013); T. Akimoto, Phys. Rev. Lett. 108, 164101 (2012).
- (67) D. Ernst, J. Koehler, and M. Weiss, Phys. Chem. Chem. Phys., 16, 7686 (2014).
- (68) Compare the discussion of the initial condition for anomalous diffusion in logarithmic potentials in A. Dechant, E. Lutz, D. A. Kessler, and E. Barkai, Phys. Rev. E 85, 051124 (2012).
- (69) F. Thiel, F. Flegel, and I. M. Sokolov, Phys. Rev. Lett. 111, 010601 (2013); Phys. Rev. E 89, 012136 (2014).
- (70) M. Magdziarz, A. Weron, K. Burnecki, and J. Klafter, Phys. Rev. Lett. 103, 180602 (2009); M. Magdziarz and J. Klafter, Phys. Rev. E 82, 011129 (2010); Y. Meroz, I. Eliazar, and J. Klafter, J. Phys. A 42, 434012 (2009); Y. Meroz, I. M. Sokolov, and J. Klafter, Phys. Rev. Lett. 110, 090601 (2013); S. Condamin, V. Tejedor, R. Voituriez, O. Bénichou, and J. Klafter, Proc. Natl. Acad. Sci. USA 105, 5675 (2008); V. Tejedor, O. Bénichou, R. Voituriez, R. Jungmann, F. Simmel, C. Selhuber-Unkel, L. Oddershede, and R. Metzler, Biophys. J. 98, 1364 (2010); M. Bauer, R. Valiullin, G. Radons, and J. Kaerger, J. Chem. Phys. 135, 144118 (2011); J.-H. Jeon, E. Barkai, and R. Metzler, J. Chem. Phys. 139, 121916 (2013).