# Aging and nonergodicity beyond the Khinchin theorem

See allHide authors and affiliations

Edited by Lawrence A. Shepp, Rutgers, The State University of New Jersey, Piscataway, NJ, and approved June 23, 2010 (received for review March 22, 2010)

## Abstract

The Khinchin theorem provides the condition that a stationary process is ergodic, in terms of the behavior of the corresponding correlation function. Many physical systems are governed by nonstationary processes in which correlation functions exhibit aging. We classify the ergodic behavior of such systems and suggest a possible generalization of Khinchin’s theorem. Our work also quantifies deviations from ergodicity in terms of aging correlation functions. Using the framework of the fractional Fokker-Planck equation, we obtain a simple analytical expression for the two-time correlation function of the particle displacement in a general binding potential, revealing universality in the sense that the binding potential only enters into the prefactor through the first two moments of the corresponding Boltzmann distribution. We discuss applications to experimental data from systems exhibiting anomalous dynamics.

- anomalous diffusion
- ergodicity breaking
- single particle trajectories
- continuous time random walk
- irreversibility

Important tools for analyzing the generally complicated dynamics of physical observables *O* are their two-time correlation functions *C*_{O}(*t*_{2},*t*_{1}). These are pair products of dynamic observables that are averaged over an ensemble, as defined below. When such a correlation function describes a stationary process, *C*_{O}(*t*_{2},*t*_{1}) is a function of the time difference only; that is, . For such processes the correlation function contains the dynamical information on other physical observables via fundamental theorems (1). For instance, the Wiener-Khinchin theorem relates the power spectrum to the correlation function, or the fluctuation-dissipation theorem connects correlation functions to linear response functions. Another well-known example is Khinchin’s theorem (2), which provides a criterion for ergodicity of a process in terms of the corresponding stationary correlation function. However, the stationarity property is found to be violated in numerous systems (see below). These systems exhibit aging properties that are intimately connected to the nonstationary behavior. Correlation functions in aging systems behave very differently from their stationary counterparts (3–6). For instance, systems may be characterized by correlation functions of the type *C*_{O}(*t*_{2},*t*_{1}) = *h*_{O}(*t*_{1}/*t*_{2}) (if *t*_{2}≥*t*_{1}); i.e., the two times enter as a quotient rather than their difference. This nonstationary behavior is also connected to a breaking of ergodicity in the sense that long-time averages differ from ensemble averages of the same quantities (3, 7–10). The relation between correlation functions and ergodicity breaking can be quantified by the Edwards–Anderson parameter (3); see below.

We here present a possible generalization of the Khinchin theorem for aging systems, namely, provide the condition for ergodicity for systems exhibiting aging. Moreover, for systems with broken ergodicity we relate fluctuations of time averages to the corresponding correlation function. In particular, we derive an analytical expression for the two-time position correlation function in the presence of an external binding potential *U*(*x*) based on the general framework of the fractional Fokker–Planck equation (11–13). The latter describes systems in which the mean squared displacement in free space scales subdiffusively as *t*^{α}, the anomalous diffusion exponent α ranging in the interval 0 < *α* < 1 (12–14). Our results for the correlation function are quite general and could also be applied to continuous time random walk (CTRW) systems in confined geometries and mean field descriptions of models such as the quenched trap model (14). In particular, we show that for sufficiently long times the *correlation function behaves universally*, and the dependence on the potential *U*(*x*) enters only in the prefactor through the first two moments of the corresponding Boltzmann distribution. Moreover, we demonstrate agreement of our results with experimental data.

Physical systems displaying nonstationary behavior such as aging and ergodicity breaking traditionally included glassy systems such as spin glasses (3), colloidal glasses (15), gels (16), turbulent systems (17), or tracer dispersion in subsurface hydrology (18), among others. More recently advanced single-molecule experiments reveal other types of complex systems with similar behavior. These are systems exhibiting anomalous diffusion and slow, nonexponential relaxation dynamics (12–14, 19). For instance, they include blinking quantum dots (20–23), or biologically relevant systems. The latter include subdiffusion of tracer particles in living biological cells (24–29) or reconstituted crowding systems (30–33), protein conformational dynamics (34), or the motion of bacteria in a biofilm (35). Here we show how the knowledge of the aging correlation function allows us to quantify the nonergodic behavior of the underlying process.

## Correlation Functions and Ergodicity

Let *x* be the state of our system—a general point in the phase space whose dynamics is given by a one parameter flow on the state space, where *t* is time (the state space is assumed to be equipped with the Borel *σ*-algebra). We refer to an observable as a measurable function on state space, *O*(*x*), which in time *t* is therefore . We then define the “finite-time average” [1]and the ensemble average [2]where μ is a stationary ensemble measure. The ergodic theorem asserts that the infinite time limit exists for μ-almost every *x*, [3]where the function is constant on ergodic components, and that if there exists only one ergodic component, then [4]for μ-almost every *x*. The connection between the ergodic theorem for the observable *O* and its two-time correlation function was established by Khinchin (2), a work of considerable impact on statistical physics (36–38). Khinchin’s theorem states that an observable *O* is ergodic if the associated correlation function is “irreversible,” in the sense that if *O* fulfills [5]then the process is ergodic, namely, Eq. **4** holds. In the derivation of Khinchin’s theorem, it is assumed that the process is stationary and the system reached a steady state, . Note that for stationary systems (39), irreversibility is a broader concept than ergodicity, and therefore Khinchin’s theorem cannot always be reversed.

Having in mind aforementioned nonstationary systems, the following question arises: Does irreversibility imply ergodicity also for aging systems? And if not, what theorem replaces Khinchin’s? To start answering these questions, we quantify ergodicity in terms of second moments, namely, if the dynamics is ergodic, as originally pointed out by Khinchin, and holds for processes where converges to a constant. From Eq. **1** we find for the second moment [6]Assuming a symmetric correlation function *C*_{O}(*t*_{1},*t*_{2}) = *C*_{O}(*t*_{2},*t*_{1}) and aging behavior of the above form *C*_{O}(*t*_{2},*t*_{1}) = *h*_{O}(*t*_{1}/*t*_{2}) (*t*_{2}≥*t*_{1}) (i.e., we use the aging regime as the “steady state” of the system similarly to invoking stationarity in Khinchin’s theorem), Eq. **6** becomes [7]After the substitution (*t*_{1},*t*_{2}) → (*z*,*t*_{2}), where *z* = *t*_{1}/*t*_{2}, we find the time-independent result [8]Eq. **8** is a general relation between the fluctuation of the time-averaged and ensemble-averaged correlation functions. For ergodicity to hold, we require for μ-almost every *x* that due to Eqs. **3** and **4**, and so, because 〈〈*O*〉^{2}〉 = 〈*O*〉^{2}, the condition [9]must be fulfilled. Simultaneously we can rewrite the irreversibility condition in the form [10]Thus, to fulfill Eq. **9** the condition **10** is not sufficient, and Khinchin’s theorem does not hold for aging systems. Namely, irreversibility does not imply ergodicity in an aging system. We propose a generalization of Khinchin’s theorem stating that in the case of irreversible dynamics the condition for ergodicity in an aging system is given by [11]Thus the knowledge of the correlation function resolves the question of ergodicity also for an aging system. More importantly, with the help of Eq. **8** the correlation function allows one to quantify the fluctuations of the time averages and ergodicity breaking: [12]where we assumed irreversibility. For the case when the aging regime starts only for long enough time this relation holds for *t* comparable or larger than the starting time of the aging regime. We now turn to the calculation of the correlation function in an aging system and prove that it is irreversible but still nonergodic. The result will be shown to exhibit universal features that, together with Eq. **12**, imply a generic behavior of the fluctuations of time averages.

## Model

### Fractional Dynamics.

We consider fractional dynamics (see below) and focus on the state variable as the observable, namely, . Note that **x**(0) = *x*. To be specific we consider anomalous diffusion in an external potential *U*(*x*) = -∫^{x}*F*(*x*^{′})*dx*^{′}, because this problem attracted considerable attention in physical sciences (11, 12). The probability density function of **x**(*t*) is governed by the fractional Fokker–Planck equation (11, 12) [13]in which the Fokker–Planck operator includes drift and diffusion terms (40), *K*_{α} is the anomalous diffusion coefficient of dimension m^{2}/ sec^{α}, and *k*_{B}*T* the thermal energy. The fractional Riemann–Liouville operator [14]involves long-range memory effects (12, 13). Eq. **13** describes the time evolution of the single-time probability density function *f*(*x*,*t*) and has been studied extensively for different potentials (12, 13, 41). Eq. **13** can be derived as the long-time limit of a continuous time random walk model, in which the local probability to jump left and right is biased by the external potential *U*(*x*) (42, 43).

The fractional Fokker–Planck equation (**13**) can be rephrased in terms of the Langevin equations (44–47) [15]where *s* is an internal time (unitless) and *t* is the physical (laboratory) time. In Eq. **15** *η*(*s*) is standard Brownian motion; *K* > 0 is the diffusivity for the normal diffusion process in internal time *s*. Conversely *ω*(*s*) represents an asymmetric Lévy-stable process of order *α* such that the probability density function of *s* is (41, 45, 46) [16]*L*_{α} is a one-sided Lévy-stable probability density function with Laplace transform exp(-*λ*^{α}) (41). The representation of the fractional dynamics by means of the coupled Langevin equations (**15**) is usually termed subordination (44). With Eq. **16** it is easy to show (41) that statistical properties of **x**(*t*) are independent of the choice of *K* (one may set *K* = 1, for example). Equations of the form (**15**) are useful simulation tools (46, 47) and were used to investigate multiple-time probability density functions (45, 48, 49) and correlation functions (50, 51).

### Two-Point Probability Density Functions.

The relation between the solution of the fractional Fokker–Planck equation and its Brownian counterpart via subordination can be used to derive the *m*-point probability density functions for a subdiffusion process (45, 48, 49). In particular, the two-point probability density function is given by (45) [17]where *n*(*s*_{2},*t*_{2}; *s*_{1},*t*_{1}) represents the two-point probability density function of the inverse Lévy-stable process *s*(*t*), and *f*_{M}(*x*_{2},*s*_{2}; *x*_{1},*s*_{1}) is the two-point probability density function of the corresponding Markovian process **x**(*s*) defined in Eq. **15**(a). The Laplace transform of *n* is given by (45) [18]where Λ ≡ *λ*_{1} + *λ*_{2} and Θ is a step function. Here and in the following, the Laplace transform for the variable pairs *t*_{1}↔*λ*_{1} and *t*_{2}↔*λ*_{2} is denoted by the explicit dependence on the respective variable.

Knowing that **x**(*s*) is a Markov process, the two-point probability density function *f*_{M}(*x*_{2},*s*_{2}; *x*_{1},*s*_{1}) of the process **x**(*s*) is obtained from the well-known property *f*_{M}(*x*_{2},*s*_{2}; *x*_{1},*s*_{1}) = *P*_{M}(*x*_{2},*s*_{2}|*x*_{1},*s*_{1})*P*_{M}(*x*_{1},*s*_{1}|*x*_{0},0), where *x*_{0} is the initial position of the particle and *P*_{M}(*x*_{1},*s*_{1}|*x*_{0},0) is single-time probability function for the process **x**(*s*) found by solving Eq. **15**(a). We decompose this expression in the form [19]From Eqs. **17**–**19**, we obtain the subdiffusive probability density function in terms of the Markovian two-point probability density function *P*_{M} in Laplace space via [20]Here denotes the Laplace transform () of *P*_{M}(*x*,*s*|*x*_{0},0), and similarly for Λ^{α}. Eq. **20** is the final result connecting the two-point probability density function *f* in Laplace space to its Markovian counterpart.

## Results

### Position Autocorrelation.

Eq. **20** allows one to calculate general two-point correlation functions for subdiffusive systems governed by the fractional Fokker–Planck equation (**13**). In particular, the position-position correlation function is [21]in Laplace space. Laplace inversion (see *SI Text*) delivers the final result for the two-time correlation function [22]valid for *t*_{2}≥*t*_{1}≫(1/*K*_{α}*γ*_{1})^{1/α}, *γ*_{1} being the smallest nonzero eigenvalue of the Fokker–Planck operator (40), and [23]is the incomplete Beta function (52). The symbol 〈·〉_{B} denotes an average over the Boltzmann distribution (i.e., the stationary measure μ in our example) [24]Fig. 1 shows the sigmoidal behavior of the position autocorrelation function. It is important to note that **22** implies that the process is irreversible because in the limit Δ → ∞, *B*(*t*/*t* + Δ,*α*,1 - *α*) → 0 and .

The result (**22**) has some remarkable properties. Thus the external potential *U*(*x*) enters only via the prefactor, and only through the first two moments of the corresponding Boltzmann distribution. The long-time behavior of *C*_{x}(*t*_{2},*t*_{1}) is universal and depends only on the ratio *t*_{2}/*t*_{1}.

We note that beta function behavior for a correlation function was found previously for a simple two-state renewal model with power-law sojourn times on both states (3, 20, 22, 23, 53, 54). While our process is clearly not a two-state process, the universal behavior of expression **22** is due to the separation between the physical process in space and the associated temporal process. Such a separation is exactly the idea behind the subordination of time, Eq. **15**(b). The temporal process *t*(*s*) yields the time evolution governed by the waiting times between successive jumps. Due to the assumption of annealed disorder it is independent of the current particle position. It converges as a function of the number of jumps (∝ *s*) due to the generalized central limit theorem, corresponding to the long-time limit *s* → ∞ in Eq. **15**(b). The limiting behavior of *t*(*s*) is therefore a Lévy-stable law that is underlying the subdiffusion dynamics. Conversely, the spatial process explores the external potential and is not affected by the disorder if observed as a function of the internal time *s*. In fact, as a function of *s* the process corresponds to normal diffusion in an external field, and so the process converges to Boltzmann statistics characterized by the binding properties of the external potential. We note that the result (**22**) for the correlation function mirrors the convergence of both temporal and spatial processes and is independent of the microscopic properties of the model [e.g., the shape of *U*(*x*)]. To obtain the behavior when one of the processes has not converged one needs to use the full correlation function with a nontrivial time dependence involving all eigen-values of .

### Properties of the Two-Time Position Correlation.

The correlation function (**22**) displays a number of noteworthy features:

#### (i) Aging behavior.

The correlation function *C*_{x}(*t* + Δ,*t*) exhibits aging because its time dependence is of the form Δ/*t*. Aging behavior was indeed observed in many complex systems (3, 16, 55), for instance, in thermoremanent magnetization experiments (55, 56), in which the measured relaxation of the magnetization *M*(*t* + Δ,*t*) is proportional to the spin correlation function, according to generalized fluctuation-dissipation theorems (3, 54, 57). Accepting our stochastic theory as an approximation for the spin system behavior we used **22** to fit the thermoremanent magnetization data from ref. 55. The result of the fit is presented in Fig. 1. We observe good agreement between the data and our beta function results over the entire time window, with a slight discrepancy at short times. We note that the use of a nonzero value for the fitting parameter in **22** for the zero external field relaxation of the magnetization is consistent with observed asymmetric magnetic fluctuations in thermoremanent magnetization experiments (58) as opposed to the naively expected zero average behavior. Fitting with the beta function for the measured correlation function does not necessarily yield insight into the physics of the system, but classification of aging with particular fitting functions might be a profitable step (as the well-known functions, such as Cole–Cole functions, are useful in the classification of dielectric relaxation).

#### (ii) Time-averaged position.

We quantify ergodicity, or the departure from ergodicity, of the system by measuring fluctuations of the finite-time-averaged position (note that both the left- and right-hand side depend implicitly on the initial position *x*, which is suppressed in our notation **x**(*t*) and, hence, suppressed on the left-hand side as well). Combining Eq. **8** and **22**, we see that [25]where we used the relation . This result was previously obtained for a CTRW process (59). Clearly when *α* < 1; thus we observe weakly nonergodic behavior. In the limit *α* = 1 ergodicity is restored.

#### (iii) Edwards–Anderson parameter.

The Edwards–Anderson parameter was previously used to quantify the degree of weak ergodicity breaking in the context of spin glasses (3). It is defined as [26]In our current framework, for the case of a symmetric potential the Edwards–Anderson parameter becomes [27]reflecting irreversibility of our process that is still nonergodic. Conversely, interchanging the limits we find that [28]reflecting the aging character of the system. Note that Eq. **27** indicates that *q*_{EA} is determined by the Boltzmann distribution and is independent of α in the nonergodic phase.

#### (iv) Time-averaged mean-squared displacement.

From **22** we also obtain the behavior of the ensemble average of the time-averaged mean-squared displacement (60, 61). Namely, from a time series **x**(*t*^{′}) recorded in single particle tracking experiments one can define the time-averaged mean-squared displacement [29]where *t* is the overall measurement time, whereas the dependence on *x* enters implicitly on the right-hand side through **x**(*t*). To unburden the notation, in the following we drop the arguments (*x*, *t* - Δ) from . At finite measurement time *t* even in the Brownian limit the quantity is a random quantity depending on the particular trajectory. Performing an additional ensemble average, for a Brownian system the role of the lag time Δ in the long measurement time limit is completely interchangeable with the regular *t* dependence in the corresponding ensemble-averaged mean-squared displacement, for example, when *U*(*x*) = 0, . In the presence of a confining potential, one would naively expect the mean-squared displacement to saturate, as observed for a Brownian system. However, evaluating the ensemble average of ,[30]we find the a priori surprising result: For regular diffusion in a binding potential one obtains a saturation for long times, as for anomalous motion in the case of the ensemble average (12, 13). In contrast, for the time-averaged mean-squared displacement in our anomalous system from Eq. **30**, we find from the beta function expansion the behavior [31]valid in the limit (Δ/*t*) ≪ 1 and for Δ≫(1/*K*_{α}*γ*_{1})^{1/α}. Instead of the naively expected saturation, the time-averaged mean-squared displacement grows as a power with exponent (1 - *α*). Only when the lag time approaches the measurement time *t* this power-law growth stops, and the function dips to the ensemble averaged value. We note that the Δ^{1-α} scaling was recently reported for the case of a particle in a box (62).

#### (v) Numerical analysis of position autocorrelation.

Fig. 2 shows, over a large time span, the time-averaged mean-squared displacement of a subdiffusing particle (i) in an harmonic potential and (ii) in a box with reflecting boundaries. The initial particle position was chosen to be at the bottom of the potential and in the center of the box, respectively. At short lag times Δ we observe the linear scaling [32]with the lag time Δ. In this result only the dependence on the overall measurement time *t* bears witness to the fact that the underlying stochastic process is subdiffusive. Seemingly paradox, the lag time Δ enters linearly, in contrast to the associated ensemble-averaged mean-squared displacement 〈(**x**(*t*))^{2}〉 ∼ 2*K*_{α}*t*^{α}/Γ(1 + *α*). However, this is the result of the weak ergodicity breaking of the process, as shown in refs. 60 and 61. The free particle behavior at short Δ is an expected result, which can be obtained from scaling arguments or explicitly from the full correlation function: At sufficiently short times the particle does not yet feel the confinement due to the reflecting boundaries, or it does not yet experience the restoring force of the potential, respectively. The regime holds for scales of the lag time Δ that fulfill *K*_{α}Δ^{α} ≪ *L*^{2} in the example of the box, where *L* is the size of the box. For a general confining potential, the turnover time is nonuniversal and is dependent on all nonzero eigenvalues *γ*_{n} of the Fokker–Planck operator (40). Thus, at times Δ≫(1/*K*_{α}*γ*_{1})^{1/α} a transition occurs to the regime, **31**. We stress again that, in contrast to normal diffusion, no saturation is found at long lag times, and continues to grow for any Δ < *t*. Only as Δ approaches to the measurement time *t*, we obtain the convergence .

In Fig. 3 we show the simulations result for a number of individual trajectories in an harmonic binding potential, displaying significant scatter. This scatter between individual trajectories is, again, a result of the weak ergodicity breaking of the underlying stochastic process and resembles qualitatively the behavior observed in single particle tracking experiments (24–29). The amplitude of the scatter can be quantified in terms of the dimensionless random variable , the relative scatter of the time average with respect to its ensemble mean. Using arguments similar to ref. 61, it can be shown that the distribution of *ξ* is given in terms of a one-sided Lévy-stable law in the form (61) [33]where the stable density *L*_{α}(*t*) is the same as for Eq. **16**. Note that the random variable *ξ* is in the denominator of *L*_{α}, and therefore the associated width is finite. For instance, for the case *α* = 1/2 used in Figs. 2 and 3 we find the Gaussian law lim _{t→∞}*ϕ*_{1/2}(*ξ*) = (2/*π*) exp(-*ξ*^{2}/2). Finally, in the Brownian limit *α* = 1, the distribution converges to the sharp behavior *ϕ*_{1}(*ξ*) = *δ*(*ξ* - 1), restoring ergodicity in the sense that no more scatter occurs. The distribution of *ξ*, given by Eq. **33**, is the same for both unbounded and bounded anomalous diffusion. This is simply due to the mentioned fact that temporal and spatial process are uncoupled.

## Discussion

Correlation functions are a standard tool to experimentally probe the temporal evolution of a system. They provide information on how the present value of a physical observable influences its value in the future and are therefore important indicators to the specific process that governs the system’s dynamics. The significance of the correlation function behavior for the fundamental concepts in physics is revealed through Khinchin’s theorem, which provides a condition for a stationary process to be ergodic in terms of the corresponding correlation function. Herein we proposed a possible generalization of Khinchin’s theorem for a class of nonstationary aging process. We provide not only a generalized condition for ergodicity for such processes in terms of the corresponding aging correlation function but also quantify the deviations from ergodicity. For the broad class of nonstationary processes described by Eq. **13**, we derived analytically the time dependence of two-point correlation functions for subdiffusing particles under situations of confinement. In particular we revealed a universal behavior for the two-time position correlation function involving the incomplete beta function. All features of the confining potential enter the correlation function solely through the prefactor in terms of the first and second moments of the associated Boltzmann distribution. Of course, the expression for the correlation function in **22** is not restricted to a position correlation function, but can be used to describe a very general class of quantities, e.g., a potential energy correlation function. The generality of our results is a direct consequence of the convergence of the Markovian process **x**(*s*) in the jump space, and the ubiquitous role of Lévy statistics due to the generalized central limit theorem.

## Acknowledgments

We thank G. G. Kenning for the permission to use the experimental data published in ref. 55. We are grateful to I. Eliazar for helpful discussions. This work was supported by the Israel Science Foundation and the Deutsche Forschungsgemeinschaft.

## Footnotes

^{1}To whom correspondence should be addressed. E-mail: metz{at}ph.tum.de.Author contributions: S.B., R.M., and E.B. designed research; S.B., R.M., and E.B. performed research; S.B., R.M., and E.B. contributed new reagents/analytic tools; S.B., R.M., and E.B. analyzed data; and S.B., R.M., and E.B. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1003693107/-/DCSupplemental.

## References

- ↵
- Forster D

- ↵
- Khinchin AI

- ↵
- Yound AP

- Bouchaud JP,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Mattson J,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Stefani FD,
- Hoogenboom JP,
- Barkai E

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Gal N,
- Weihs D

- ↵
- ↵
- Seisenberger G,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- Yang H,
- et al.

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Risken H

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Barkai E,
- Sokolov IM

- ↵
- ↵
- ↵
- ↵
- Abramowitz M,
- Stegun I

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Physics