Spurious ergodicity breaking in normal and fractional Ornstein–Uhlenbeck process

The Ornstein–Uhlenbeck process is a stationary and ergodic Gaussian process, that is fully determined by its covariance function and mean. We show here that the generic definitions of the ensemble- and time-averaged mean squared displacements fail to capture these properties consistently, leading to a spurious ergodicity breaking. We propose to remedy this failure by redefining the mean squared displacements such that they reflect unambiguously the statistical properties of any stochastic process. In particular we study the effect of the initial condition in the Ornstein–Uhlenbeck process and its fractional extension. For the fractional Ornstein–Uhlenbeck process representing typical experimental situations in crowded environments such as living biological cells, we show that the stationarity of the process delicately depends on the initial condition.


Introduction
The Ornstein-Uhlenbeck process is one of the most fundamental physical processes, originally devised to describe the velocity distribution and relaxation of a Brownian particle under the influence of a velocity-dependent friction. The Ornstein-Uhlenbeck process belongs to the class of Gaussian and Markovian processes, and it is described in terms of the stochastic Langevin equation [1][2][3] 4 dx + λxdt = σdB t · (1) Here dB t is the increment of the well-known Brownian motion (Wiener process) B t , and λ and σ are positive constants. 1/λ defines a natural dynamic time scale, and σ is the intensity of the fluctuations. Under certain conditions discussed below the Ornstein-Uhlenbeck process is the only non-trivial process in the class of Gauss-Markov processes that has a stationary solution [5]. Physically, overdamped Brownian particles in an optical tweezers trap [6] or tethered to an anchor by a flexible polymer [7] are adequately described in terms of an Ornstein-Uhlenbeck process. The Ornstein-Uhlenbeck process is also used as a phenomenological model for the confinement observed in the tracer diffusion in critical random environments [8]. A wide field of applications of the Ornstein-Uhlenbeck process lies in finance. The Ornstein-Uhlenbeck process was adopted in 1970s by Vašíček to model the evolution of the interest rate of financial markets [4]. Extending this Vašíček model, Hull and White took into account explicitly time dependent drift μ and λ [9]. There are other variants of the Vašíček model, for instance, the jump-extended Vašíček model in which an exponential jump noise following a Poisson distribution is added to equation (1) [10]. There also exist extensions of the Ornstein-Uhlenbeck process to non-Gaussian processes with applications in finance [11], including option pricing [12], commodity derivative pricing [13] and electricity pricing [14]. Such models have also been utilised to model neural activity [15] or to study the statistics of neuron spikes [16]. The Ornstein-Uhlenbeck process corresponds to the continuous-time analogue of a discrete-time autoregressive AR(1)-process [17][18][19].
In a direct extension of the Ornstein-Uhlenbeck process (1) one replaces the white Gaussian noise dB t by power-law correlated fractional Gaussian noise [20]. In absence of the damping term this so-called fractional Brownian motion captures the motion of diffusive particles in viscoelastic environments, such as artificially crowded media [21][22][23], lipid bilayer membranes [24][25][26], or the cytoplasm of living biological cells [27][28][29]. The correlations in the noise effect anomalous diffusion of the form x 2 (t) t α [30,31]. Combined with a Hookean restoring force exerted by optical tweezers a tracer particle in a biological cell [28,32] then follows the fractional Ornstein-Uhlenbeck process [33]. Formally, the fractional Ornstein-Uhlenbeck process is still Gaussian and stationary, yet it is strongly non-Markovian. As we will see this causes fundamental differences. We note that in finance power-law correlations are frequently observed in the dynamics of stock structure and price dynamics [34], commodity prices [35], and return of the closing values of the financial indices [36][37][38]. There exist several studies modelling such long-ranged correlated processes with ARFIMA, GARCH, and FIARCH processes, and to quantify the cross-correlation of mutually dependent processes [38,39].
With modern microscopic techniques it is possible to track single sub-micron tracer particles and even single molecules through complex media such as live biological cells [6,40]. The time-series extracted from such single-particle trajectories are typically evaluated in terms of time-averaged physical observables [41,42]. To address the motion of a Brownian or fractional Brownian particle under the action of an external potential by analysing a single trajectory of its movement, it is essential to understand whether the physical process governing the motion of the particle is ergodic, or not [31,43,44]. To infer the ergodic property of a given Gaussian process it is sufficient that the associated two-time covariance function solely depends on the difference of the two times [45]. This property rests on the fact that for Gaussian processes all properties can be deduced from the mean and covariance function [46,47]. An indirect approach to deduce the ergodic property of the process is to compare the behaviour of the mean squared displacement (MSD) and the time averaged MSD [31,[48][49][50].
We here scrutinise the exact ergodic and stationary behaviour of the regular and fractional Ornstein-Uhlenbeck processes and show that they fundamentally differ in some of their behaviour, despite of the fact that both are ergodic. In particular, we elucidate the precise role of the initial condition and invalidate the general belief that the assertion of an equilibrium initial condition necessarily recovers the stationary property of the process. We first analyse the detailed statistical properties from the covariance of the Ornstein-Uhlenbeck process in section 2, including the ensemble-and time-averaged MSDs and the effect of the initial condition. Section 3 provides an analogous analysis for the fractional Ornstein-Uhlenbeck process. In section 4 we discuss our results and conclude. Some mathematical details are deferred to the appendix.

Ornstein-Uhlenbeck process
We define the Ornstein-Uhlenbeck process in terms of the stochastic differential equation (1), in which dB t is the increment of Brownian motion B t with the covariance function [51] Cov In this formulation, Gaussian white noise corresponds to the time derivative of the increment, dB t /dt. After solving the stochastic differential equation (1), x(t) is formally obtained as where x 0 = x(t = 0) defines the initial condition. Since B t is a continuous process, via integration by parts the above equation is recast into with B 0 = 0. The MSD for a random process x(t) is defined as where The MSD can also be written in terms of the covariance function, For the Ornstein-Uhlenbeck process the MSD then assumes the following expression, where Var(X) stands for the variance of a random variable X. Note that in the limit λ → 0 of free Brownian motion this notation leads to the MSD lim λ→0 Ω 2 (t) = σ 2 t. The time-averaged MSD (TAMSD) is defined as [31,41,48] where T is the total measurement time and Δ is called the lag time. For the Ornstein-Uhlenbeck process, the TAMSD yields in the form

Properties of the Ornstein-Uhlenbeck process
Since the Ornstein-Uhlenbeck process is a Gaussian process, it suffices to know the covariance function and mean to infer its properties. The mean of x(t) according to equation (4) is Recall that a Gaussian process is stationary and ergodic if the covariance function at two times exclusively depends on the time difference, that is, Cov(x(t 1 ), x(t 2 )) = G(|t 2 − t 1 |) in terms of the continuous function G. The covariance (10) satisfies the requirements of stationarity if (i) t 1 or t 2 are significantly larger than 1/λ, or (ii) if Var(x 0 ) = σ 2 2λ . The first condition is asymptotic with respect to 1/λ: the process loses the memory of its initial condition after the correlation time 1/λ. The second condition is valid for all times, it corresponds to starting the process with the equilibrium distribution.
The equilibrium stationary distribution can be deduced from the Fokker-Planck equation of the Ornstein-Uhlenbeck process [52] ∂P(x, t) where P(x, t) is the probability density function of the process. The solution for P(x, t) is [3,52] P(x, t) = λ πσ 2 In the stationary limit t 1/λ the stationary probability density function is given by P(x) = [λ/(πσ 2 )] −1/2 exp(−λx 2 /σ 2 ), for which the variance becomes Assume that the distribution of x 0 satisfies the stationary distribution, Var(x 0 ) = σ 2 /(2λ), from equations (7) and (9) one arrives at The fact that MSD and TAMSD are equivalent for an equilibrium stationary initial distribution is the direct consequence of the stationary property of the process, that can be directly inferred from the covariance (10). Thus, MSD and TAMSD indeed coincide. Yet there exists an intrinsic problem regarding the way MSD and TAMSD are defined, and the equivalency between the two is only valid under the strict conditions that the equilibrium initial condition is met and that the process is both Gaussian and Markovian, see the next section for the non-Markovian fractional Ornstein-Uhlenbeck process.
Compare the pairs of equations (7) and (9) as well as (14) and (15). In the first pair, (7) and (9), we note the existence of two time scales in the TAMSD, Δ and T, the latter of which does not exist in the definition of the MSD. This effects a disparity in inferring the stationary state in a consistent way from MSD and TAMSD. From the MSD, the stationary state is reached when the expression ceases to depend on t, that is when t 1/λ. In contrast, for the TAMSD the stationarity condition depends on the interplay between lag time Δ and measurement time T. The observation time T identifies the total time the process has been monitored to evolve, and one identifies the stationary state of the process when T 1/λ. Yet Δ signifies the magnitude of the time window in the sliding average, comparing two instances of the process. Necessarily, Δ < T, however, also the lag time Δ needs to be compared with the natural dynamical time scale imposed by 1/λ. There exist two distinct regimes: (i) if Δ is much smaller than 1/λ the fluctuations present in the system during this time interval have not relaxed. Therefore any statistical inference cannot be justified although the overall process has reached stationarity for T 1/λ. (ii) Stationarity is reached when T 1/λ and Δ 1/λ, as long as Δ T is simultaneously fulfilled. Obviously, for the trivial case T < 1/λ the process cannot be stationary. When the initial condition is chosen to be the equilibrium distribution, Var(x 0 ) = σ 2 /2λ, we see from equations (14) and (15) that the situation is different: here stationarity is reached once t 1/λ for the MSD and Δ 1/λ for the TAMSD. Note that the signature of T disappears (an indication of stationarity). The caveat here is that, for the MSD, asserting the equilibrium initial condition, which implies the stationary property of the process, does not imply the independence of the MSD of t, in contrast to the case of the TAMSD, in which the dependency on T disappears.
This discrepancy also manifests itself in the asymptotes of MSD and TAMSD when the equilibrium initial condition is not asserted. The asymptotes of MSD and TAMSD in the stationary state read Indeed, for the MSD the stationary value asymptote depends on the variance Var(x 0 ) of the chosen initial distribution. This contradicts the common intuition that, once the process reaches its stationary state, any trace of the initial condition must have vanished. In contrast, Var(x 0 ) is absent from the limiting value of the TAMSD. Knowing that the Ornstein-Uhlenbeck process is stationary and ergodic, these observables, suggesting non-ergodic behaviour, are thus unsuitable. In particular, the above difference could potentially lead to wrong conclusions for the ratio of noise strength σ 2 and trap strength λ depending on which measure is chosen for the evaluation of an experiment. This discussion elucidates the fundamental difference between the generic definitions of the MSD and the TAMSD, essentially quantifying different properties of a random process. Thus, while the MSD quantifies the dispersal of an ensemble of walkers at a given time instant t with respect to the initial condition, the TAMSD quantifies how increments of the process evolve as function of the lag time. We now embark for modified definitions of these most widely used physical observables for stochastic processes for the case of the Ornstein-Uhlenbeck process.

Generalised definitions of the ensemble-averaged MSD
We propose to recalibrate the definition of the MSD in the generalised form where the subscript Δ indicates the generalisation. This modified MSD describes the dispersal of the process from time t to t + Δ; in other words, the dispersal of increments in which the mean effect of the initial condition and the drift are removed. We can rewrite expression (18) in terms of the covariance function in the form In the limit λ → 0 of free Brownian motion, this definition produces lim λ→0 Ω 2 Δ (t) = σ 2 Δ, which is the same expression as obtained for the classical definition, albeit with t replaced by Δ. In this generalised formulation the integrand of the TAMSD (8) is exactly the generalised expression of the MSD given by equation (18), that is, (20) which readily yields equation (9). We observe that for equilibrium initial conditions, Var(x 0 ) = σ 2 /2λ, the generalised expressions for MSD and TAMSD yield exactly the same result (σ 2 /λ) [

Fractional Ornstein-Uhlenbeck processes
The fractional Ornstein-Uhlenbeck process is the extension of the normal Ornstein-Uhlenbeck process (1), in which the increments of Brownian motion are substituted by the increments of fractional Brownian motion, B H t . Here H is the Hurst exponent, which is allowed to vary in the interval H ∈ (0, 1] [20]. The fractional Ornstein-Uhlenbeck process is therefore given by the stochastic differential equation [33] Here d B H t is the increment of fractional Brownian motion B H t . 7 The tilde is introduced here to denote the extension of fractional Brownian motion to the negative time domain, such that Fractional Brownian motion with Hurst parameter H ∈ (0, 1], is a continuous centred Gaussian process defined by the covariance function [20] Cov For H = 1/2, B  (21) is Since fractional Brownian motion is a continuous process this integral exists [59]. Integrating by parts the equation above can be rewritten in terms of B H t in the form Then the covariance function becomes We note that while free fractional Brownian motion, corresponding to the limit λ → 0, is ergodic [49,60,61], transient non-ergodicity occurs when the process is confined. Namely, for an harmonic external confinement (the Ornstein-Uhlenbeck process, that is) it was shown analytically and experimentally that the relaxation of the MSD is exponential while a slower power-law relaxation is observed for the TAMSD [23,62]. When x 0 is fixed or the distribution of x 0 is independent of the fractional Brownian motion all terms involving x 0 B H t vanish. Therefore the covariance function simplifies to After calculating the integrals (see appendix A for details) the covariance of the fractional Ornstein-Uhlenbeck process reads where M (a, b, z) is Kummer's function of the first kind (the confluent hypergeometric function of the first kind [63]). The integral representation of this function is given by For H = 1/2 the covariance function (28) consistently reduces to expression (10) of the regular Ornstein-Uhlenbeck process (note that M(2, 3, x) = 2(1 − e x + xe x )/x 2 ). On closer inspection of the covariance function, unlike for the case of the regular Ornstein-Uhlenbeck process above, in which the equilibrium distribution of the initial condition yields a stationary covariance function (see equation (10)), we notice that there is no possible form for Var(x 0 ) such that the covariance function (28) would exclusively depend only on the time difference between the two time points of the process. In other words, there is no initial condition, such that Cov(x(t 1 ), x(t 2 )) = G(|t 2 − t 1 |) for any given t 1 and t 2 . Asserting an equilibrium initial condition does not fulfil the requirement of an ergodic and stationary process for any t 0. Indeed, let us assume that the x 0 have an equilibrium distribution corresponding to the normal distribution N (0, ξ 2 ) with variance [33,53,59], Here the integral is calculated in appendix D. This result can also be obtained by recalling that Integration by part leads to the expression which yields the result (30). Observe that by substituting the variance of x 0 in equation (28) the covariance function would still depend on the absolute times t 1 and t 2 . To provide a hint why this is the case, recall our earlier assumption on x 0 . Our assumption that x 0 and B H t are not correlated yielded a covariance function which is not stationary for finite t 1 and t 2 . It asymptotically approaches the stationary covariance function when t 1 and t 2 tend to infinity. Furthermore, observe that x 0 and B H t are correlated in the case of fractional Ornstein-Uhlenbeck process, since the driving noise has a long-range memory. This is also reflected in the generalised MSD and TAMSD. Since the closed analytical expressions for the generalised MSD and TAMSD are too cumbersome to be presented here, we refer to Appendix B and observe that indeed the generalised expressions for MSD and TAMSD differ from one another. As we show now, in the stationary state ergodicity is indeed fulfilled.
To proceed, we note that the fractional Ornstein-Uhlenbeck process has the stationary solution [53] x indicated by the subscript s. Note that to achieve this stationary solution the domain of t has been changed to t ∈ (−∞, ∞). For this case from which it is inferred that every stationary solution x s (t) of the Langevin equation (21) has the same distribution as x(t) in the long-time limit. Consequently, we deduce that the covariance function for the stationary solution is given by (see appendix C for details) Obviously, the covariance function for the stationary solution depends only on the time difference between the two time points. With the use of equations (19) and (35) the generalised MSD and TAMSD are given by From this equivalency we conclude that the fractional Ornstein-Uhlenbeck process is ergodic in the sense of the generalised MSD. Figure 1 details the functional behaviour of the different MSDs. In the left panels for the non-stationary case, as expected the disparity between the generic MSD (5) and the TAMSD (8) is distinct. In contrast, using the generalised MSD (18) for the stationary solution the expected ergodic behaviour is restored.
For completeness, figure 2 shows how the two different versions of the MSD and the TAMSD approach the plateau value for different values of H. As can be seen for normal diffusion with H = 1/2 the relaxation is always exponential. In contrast, we recover a power-law relaxation for the TAMSD and for the generalised definition of the MSD. While this power-law form for the TAMSD was discussed earlier [62] and verified experimentally [23], the full agreement between the TAMSD and the generalised MSD is a distinct behaviour following from our definition (18) here.

Conclusions
It is commonly assumed that asserting equilibrium initial condition is sufficient and necessary for a confined stochastic process to remain stationary at all times t 0. We here demonstrated that for the case of the fractional Ornstein-Uhlenbeck process this is in fact not true. Generally, for any process which is not a Markov process one should bear in mind that due to long range correlations the assumption that the process is stationary requires one to take into account the entire history of the system. Therefore, asserting any assumption on the initial condition of the process would perturb the stationary state of the process, even in the case when this initial condition is the equilibrium distribution.
Moreover, we revealed another subtle point on how to define the stationary state of the process based on generalised definitions of the MSD and the TAMSD. While it is often believed that the sufficient condition to infer that the process has reached its stationary state is given when in the TAMSD the observation time tends to infinity. In this statement, though, it is neglected that Δ needs to be considered, as well. Indeed, while the lag time should remain significantly below the observation time, Δ T, the lag time needs to be much larger than the natural dynamic time scale of the process, Δ 1/λ. The Ornstein-Uhlenbeck process and its fractional extension are essential in modelling physical systems in the presence of an external potential. They are Gaussian processes with the difference that in the former case, the correlations are short-lived (Markov process) while in the latter case the correlations are long-ranged. It was further demonstrated that the Ornstein-Uhlenbeck process is stationary for all t 0 if the equilibrium initial condition is asserted. In contrast, this does not hold true for the fractional Ornstein-Uhlenbeck process due to the fact that the process is not Markovian.
These results will also be important for the correct analysis of measured trajectories of generic processes driven by fractional Gaussian noise in terms of the TAMSD, for instance, under confinement [64]. Moreover, the finite-time ergodic properties of the normal Ornstein-Uhlenbeck process as studied in [65,66] should be considered in view of the generalised definitions of the MSD and TAMSD provided here.

Acknowledgments
We acknowledge funding from the Deutsche Forschungsgemeinschaft (DFG), Grant Number ME 1535/7-1. RM acknowledges the Foundation for Polish Science (Fundacja na rzecz Nauki Polskiej, FNP) for support within an Alexander von Humboldt Honorary Polish Research Scholarship. We acknowledge the support of the German Research Foundation (DFG) and the Open Access Publication Fund of Potsdam University.

Appendix A. Covariance function of the non-stationary fractional Ornstein-Uhlenbeck process
To calculate the covariance function given by equation (27) two types of integrals need to be calculated. One are the single integrals with respect to either t 1 or t 2 . The other is the double integral with respect to the t 1 and t 2 . Throughout the integrations, it is always assumed that t 2 > t 1 for simplicity. Whenever the difference between the two times is relevant, the result is written in terms of the modulus.
The single integral with respect to t 1 is given by Following the same procedure one arrives at a similar expression for the integral with respect to t 2 The second type of the integrals appearing in the covariance of fractional Ornstein-Uhlenbeck process is given by Hence we arrive at the following expression for the double integral,

Appendix B. MSD and TAMSD of the non-stationary fractional Ornstein-Uhlenbeck process
After arriving at the covariance function for the non-stationary solution (28) the MSD of the fractional Ornstein-Uhlenbeck processes is deduced from Before calculating the TAMSD it is worthwhile checking that in the long time limit the expression above coincides with the earlier equation (36). In the limit t → ∞ the expression in the last square brackets remains unchanged while for the first and second square brackets it is only the second term which contributes to a non-zero value, namely, (t + Δ) 2H+1 M(2H + 1, 2H + 2, −λ(t + Δ)) and t 2H+1 M(2H + 1, 2H + 2, −λt). Considering the latter in the aforementioned long-time limit, where γ(s, x) is the lower incomplete Gamma function. Therefore, in the limit t → ∞ one indeed consistently recovers the expression of the generalised MSD for the stationary solution of the fractional Ornstein-Uhlenbeck process, equation (36). The complexity in the integration of the generalised MSD for the TAMSD is due to terms of the kind t 2H+1 M(2H + 1, 2H + 2, −λt) and t 2H+1 e −2λt M(2H + 1, 2H + 2, λt). The integration of such terms can be achieved as follows, After changing the order of the integration, Similarly, the second type of integration can be performed along the following steps, Analogously, And lastly, The final result for the TAMSD is then given by While the expressions for the TAMSD and the generalised MSD differ, both share the same asymptote in the long-time limit t, T → ∞. Moreover, the disparity between both is expected when the system has not yet reached the stationarity.

Appendix C. Covariance of the stationary fractional Ornstein-Uhlenbeck process
For the derivation of the covariance (35) we calculate the following integrals Following the same procedure one obtains the expression below for the integral with the differential dt 1 , The second type of integrals appearing in the covariance function is given by e λ(t 1 +t 2 ) (t 2 − t 1 ) 2H dt 1 dt 2 · The first two integrals are easily evaluated. In the last two integrals, by changing of variable t 1 − t 2 = q one arrives at The evaluation of the first double integral is straightforward and it yields e 2λt 1 2λ 2H+2 Γ(2H + 1). Evaluating the second double-integral requires changing the order of integration. For the last integral one arrives at the following, − e 2λt 2 2λ(2H + 1) (t 2 − t 1 ) 2H+1 M(2H + 1, 2H + 2, λ(t 1 − t 2 )) + e 2λt 1 2λ(2H + 1) (t 2 − t 1 ) 2H+1 M(2H + 1, 2H + 2, λ(t 2 − t 1 )).
Summing up all the calculations, we obtain