1/f^{\beta} noise in a model for weak ergodicity breaking

# 1/fβ noise in a model for weak ergodicity breaking

Markus Niemann Ivan G. Szendro Holger Kantz Max-Planck-Institut für Physik komplexer Systeme, Noethnitzer Straße 38, 01187 Dresden, Germany
###### Abstract

In systems with weak ergodicity breaking, the equivalence of time averages and ensemble averages is known to be broken. We study here the computation of the power spectrum from realizations of a specific process exhibiting noise, the Rebenshtok–Barkai model. We show that even the binned power spectrum does not converge in the limit of infinite time, but that instead the resulting value is a random variable stemming from a distribution with finite variance. However, due to the strong correlations in neighboring frequency bins of the spectrum, the exponent can be safely estimated by time averages of this type. Analytical calculations are illustrated by numerical simulations.

###### keywords:
noise, weak ergodicity breaking, spectral estimation
###### Pacs:
05.10.Gg, 05.40.-a, 05.45.Tp
journal: Chemical Physics

## 1 Introduction

In recent years, the study of long range correlations in experimental data has moved into the focus of interest. Whereas a large amount of empirical results has been obtained using detrended fluctuation analysis PengHavlinStanleyEtAl95 (), the same information is, in principle, also contained in the power spectrum. Observations of nontrivial power law decays of the power spectrum have a long tradition and are usually denoted by -noise, or, more precisely, with . Examples include noise in electric resistors, in semiconductor devices (flicker noise), and fluctuations in geophysical data Weissman88 (); HoogeKleinpenningVandamme81 (); BalascoLapennaTelesca02 ().

By weak ergodicity breaking one describes the behavior of a system which is ergodic in the sense that a single trajectory is able to explore completely some invariant component of the phase-space (it is indecomposable), but where the invariant measure is not normalizable, so that the time scales for this exploration might diverge (in other words, there is no finite microscopic time scale). This concept goes back to Bouchaud Bouchaud92 (). A consequence of this is that time averages computed on a single (infinite) trajectory will not converge to the ensemble mean. Instead, the time average is a random variable itself, as a function of the initial condition of the trajectory, which is drawn from a distribution with non-zero variance. It has been shown that Continuous-Time Random Walks BelBarkai05 (), but also specific deterministic dynamical systems BelBarkai06a () exhibit this behavior. A natural consequence of it is the non-existence of a finite correlation time, reflected by a power law decay of the auto-correlation function, which translates itself into a -behavior of the power spectrum. This type of spectrum has been derived for a two state process of this type by Margolin and Barkai (MargolinBarkai06, , Eq. (33)).

As said, the crucial aspect of systems with weak ergodicity breaking is that a time average in the limit of an infinite trajectory does not converge to a sharp value but is a random variable with a distribution of finite width. This should also hold for the power spectrum. Therefore, it is not evident that an estimation of from a numerically computed power spectrum of a single trajectory is reliable. We therefore investigate the estimation of the power spectrum from single realizations of a specific process with weak ergodicity breaking.

A relevant remark is necessary: If one estimates the power spectrum by a discrete Fourier transform from a time series with sampling interval (i.e., uses the periodogram), then the variance of the estimator for the power contained in the th discrete frequency is . However, if the correlations in the time domain decay exponentially fast, then the errors in adjacent frequencies and are sufficiently weakly correlated such that the error of a binned spectrum where one takes averages over adjacent frequencies decays like (PressTeukolskyVetterlingEtAl86, , Sect.13.4). In the limit of an infinite time series, the estimation error of the binned power spectrum of an ergodic process therefore decays to zero, i.e., the spectrum assumes uniquely defined values. In the following, we will discuss results for binned power spectra of processes with weak ergodicity breaking. The binning implies that in the analytical calculations we are allowed to ignore the time discretization which is present in every numerical data analysis.

The next Sections are devoted to an analytical derivation of the fact that the binned power spectrum of a weakly non-ergodic process is a random variable whose distribution has a non-zero variance. In addition, we compute the correlations between different bins. We will show that because of these correlations, further binning cannot reduce the uncertainty about the true power, and moreover, that the whole uncertainty about the power spectrum reduces to a random normalization factor. This has the important consequence that the power law exponent can indeed be estimated numerically regardless of the difficulties to estimate the spectrum itself.

In Section 5 we illustrate the analytical (asymptotic) results by numerical simulations, also opposing time averages to ensemble averages. Note that because of the long range correlations, cutting a long trajectory into pieces and treat these as an ensemble is not valid, since this ensemble would not represent an independent sample of initial conditions.

## 2 The model and basic approach

Rebenshtok and Barkai introduced the following model for a thermodynamic system showing weak ergodicity breaking RebenshtokBarkai07 (); RebenshtokBarkai08 (): In its simplest form it is characterized by two distributions, one waiting time distribution with density and one distribution of an observable according to the probability density . In the model let be i.i.d. random variables distributed according to and let be i.i.d. waiting times distributed according to . The process takes the value in the time and for the next time interval of length . In general

 X(t)=χi for Ti−1≤t

with

 Ti=i−1∑j=0τj. (2)

We assume that the first four moments of are finite, i.e.,

 μi=∫dxxiκ(x)<∞for i=1,2,3,4. (3)

Moreover, we assume that the distribution is centered, i.e., . The waiting time distribution should be in the domain of normal attraction of an one-sided Lévy stable distribution with exponent (). Therefore, the Laplace transform of ,

 ^ϕ(λ)=∫dte−λtϕ(t), (4)

can be expanded as

 ^ϕ(λ)=1−(νλ)α+o(λα)as λ→0+ (5)

with the scaling parameter .

The Fourier transform of the time series up to time is given by

 FT(ω)=∫T0dteiωtX(t). (6)

The spectrum may be estimated from

 ST(ω)=1n(T)FT(ω)FT(−ω). (7)

The function denotes a normalization. Normally, it takes the form , but this model requires to take for a non trivial spectrum. For finite time series one would have to add correction terms for the fact that the basis for Eq. (7) is a biased estimate of the correlation function (e.g., see (KammeyerKroschel89, , chapter 8.1)). These terms vanish for the asymptotic case which we are considering here. If we have an ensemble of realizations, one can look at the unbinned spectrum

 Sub(ω)=limT→∞ST(ω). (8)

Denoting the ensemble average by , we are interested in the expectation value of the spectrum

 ⟨Sub(ω)⟩=limT→∞1Tα⟨FT(ω)FT(−ω)⟩ (9)

and the covariances

 ⟨Sub(ω1)Sub(ω2)⟩=limT→∞1T2α⟨FT(ω1)FT(−ω1)FT(ω2)FT(−ω2)⟩. (10)

The limits Eqs. (9) and (10) can be rewritten to

 ⟨Sub(ω)⟩=limr→∞1rα⟨FrT1(ω)FrT2(−ω)⟩∣∣∣T1=T2=1,⟨Sub(ω1)Sub(ω2)⟩=limr→∞1r2α⟨FrT1(ω1)FrT2(−ω1)FrT3(ω2)FrT4(−ω2)⟩∣∣∣T1=T2=1T3=T4=1. (11)

In the following, it is helpful to define the two- and four point correlations

 C2(t1,t2)=⟨X(t1)X(t2)⟩,C4(t1,t2,t3,t4)=⟨X(t1)X(t2)X(t3)X(t4)⟩ (12)

and their (double resp. quadruple) Laplace transforms

 ^C2(λ1,λ2)=L[C2(t1,t2)]=∫dt1∫dt2e−λ1t1−λ2t2C2(t1,t2), (13)
 ^C4(λ1,λ2,λ3,λ4)=L[C4(t1,t2,t3,t4)]=∫d4te−λtC4(t). (14)

The double resp. quadruple Laplace transforms of the expressions in Eq. (11) are

 L[⟨FT1(ω)FT2(−ω)⟩]=1λ1λ2^C2(λ1−iω,λ2+iω),L[⟨FT1(ω1)FT2(−ω1)FT3(ω2)FT4(−ω2)⟩]=1λ1λ2λ3λ4^C4(λ1−iω1,λ2+iω1,λ3−iω2,λ4+iω2). (15)

The limits in Eq. (11) can also be expressed in Laplace space using a formulation following the multidimensional Tauberian theorem by Drozhzhinov and Zav’jalov DrozhzhinovZavjalov80 (); Drozhzhinov83 () as

 L[limr→∞1rα⟨FrT1(ω)FrT2(−ω)⟩]=limζ→0+ζαλ1λ2^C2(ζλ1−iω,ζλ2+iω), (16) (17)

In the next section, we will show how one can calculate the Laplace transforms of the multi point correlations .

However, when only one time series is given, taking from Eq. (7) as estimate of the spectrum will lead to unsatisfactory results as the value of will fluctuate except for some pathological cases. One often relies to the concept of binning. Instead of considering a single frequency , one averages all available frequencies in an given interval . For a time series of length the Fourier transform returns values at frequencies which are evenly spaced with distance . For large the sum can be approximated by an integral and we can consider as observable for the binned spectrum

 ST(ω,Δω)=1Δω∫ω+12Δωω−12Δωdω′ST(ω′). (18)

The spectrum can be estimated from

 Sbin(ω,Δω)=limT→∞ST(ω,Δω). (19)

Analytically, it is helpful to let the bin size go to zero as a last step

 Sbin(ω)=limΔω→0S(ω,Δω). (20)

It is important to take the last limit at the very end. In most of the important cases, the value of will converge almost surely to the value of the spectrum.

In this paper, we want to consider the situation where the stochastic process is the model by Rebenshtok and Barkai. The estimation of the spectrum is connected with averaging over the time series, therefore the question arises how the weak ergodicity breaking affects the observable . Does the value of converge to a single value or does it converge to a non trivial probability distribution? It will turn out that the latter is the case. Unlike the observable discussed by Rebenshtok and Barkai RebenshtokBarkai07 (); RebenshtokBarkai08 (), it seems not to be possible to determine directly the probability distribution of . But one can obtain important information from the moments and .

The first step is to look at with :

 ⟨Sbin(ω,Δω)⟩=limT→∞T−αΔω∫ω+12Δωω−12Δωdω′⟨FT(ω′)FT(−ω′)⟩. (21)

We proceed similarly to the unbinned case and formulate the limit in Laplace space (note that the Laplace transform and the -integral commute by Fubini’s theorem and the estimate where is defined in Eq. (3))

 L[limr→∞r−α1Δω∫ω+12Δωω−12Δωdω′⟨FrT1(ω′)FrT2(−ω′)⟩]=limζ→0+ζαΔω∫ω+12Δωω−12Δωdω′1λ1λ2^C2(ζλ1−iω′,ζλ2+iω′). (22)

Similarly, with and can be obtained from

 ⟨Sbin(ω1,Δω1)Sbin(ω2,Δω2)⟩=limr→∞r−2α1Δω1Δω2∫ω1+12Δω1ω1−12Δω1dω′1∫ω2+12Δω2ω2−12Δω2dω′2⟨FrT1(ω′1)FrT2(−ω′1)FrT3(ω′2)FrT4(−ω′2)⟩∣∣∣T1=T2=1T3=T4=1. (23)

One gets

 L[limr→∞r−2αΔω1Δω2∫ω1+12Δω1ω1−12Δω1dω′1∫ω2+12Δω2ω2−12Δω2d% ω′2⟨FrT1(ω′1)FrT2(−ω′1)FrT3(ω′2)FrT4(−ω′2)⟩]=limζ→0+ζ2αΔω1Δω2∫ω1+12Δω1ω1−12Δω1dω′1∫ω2+12Δω2ω2−12Δω2dω′21λ1λ2λ3λ4^C4(ζλ1−iω′1,ζλ2+iω′1,ζλ3−iω′2,ζλ4+iω′2). (24)

Therefore the problem can be reduced to the determination of the Laplace transforms of the point correlations functions. We describe in the next section how to obtain them.

## 3 A diagrammatic approach to the multi point correlation functions

In NiemannKantz08 () a diagrammatic method was introduced how to determine the joint probability distributions and multi point correlation functions of a continuous-time random walk. In this section we are adapting this method to determine . For any set of natural numbers , one defines

 Cn[q](t)=⟨n∏k=1χqk1[Tqk,Tqk+1[(tk)⟩ (25)

with the indicator function ( being any subset )

 1A(t)={1for t∈A,0for t∉A. (26)

The can be considered as partial correlations which only contribute, when the time is in the th waiting time (for ). Therefore,

 Cn(t)=⟨X(t1)⋯X(tn)⟩=∞∑q1,…,qn=0Cn[q](t). (27)

Define the contribution of one step as

 η(j)n[q](t)=∏i:qi=jχj1[0,τj[(ti)∏i:qi>jδ(ti−τj)∏i:qi

One obtains by induction

 Cn[q](t)=⟨★jη(j)n[q](t)⟩ (29)

where denotes the convolution in the time variable and runs from to any natural number larger than .

Transforming into Laplace space gives

 ^Cn[q](λ)=L[Cn[q](t)]=∏j⟨^η(j)n[q](λ)⟩. (30)

Using the following notations from NiemannKantz08 () (suppressing the dependence on )

 Vj={i:qj=i},Ej={i:qij}. (31)

The set (“vertex”) contains the indices which are at position , the set (“earlier”) contains the indices which are before that position and the set (“later”) contains the indices which are after that position. With the notation

 ΛJ=∑j∈Jλj (32)

one gets

 ⟨^η(j)n[q](λ)⟩=⟨∫dnte−λt∏i∈Vjχj1[0,τj[(ti)∏i∈Ejδ(ti−τj)∏i∈Ljδ(ti)⟩=μ|Vj|∏v∈Vjλv∑J∈P(Vj)(−1)|J|^ϕ(ΛJ∪Lj). (33)

Here, the power set of is denoted by and is the number of elements of . Therefore, the last sum goes over all subsets of indices of .

As in NiemannKantz08 () we introduce diagrams which denote the relative ordering of the s. A generic diagram is given in Fig. 1. The main idea behind this consists in grouping all steps with the same , and . If then there can be several successive steps with the same and . These successive steps are combined to a horizontal line. If then there can be only one step with these sets of indices. This is represented by a vertex with the indices leaving. Therefore the diagram Fig. 1 represents all with the following properties: First there is any number of steps with , and . A single step with , and follows. At next is any number of steps with , and . The rest in interpreted analogously. For each point correlation there is only a finite number of diagrams and summing over all represented by a diagram turns out to be easy. It turns out that one can associate to each part of the diagram a factor (for a more comprehensive derivation of these results, we refer to NiemannKantz08 ()). One obtains the following rules (where is defined in Eq. (3))

 γvertex i=μ|Vi|∏v∈Viλv∑J∈P(Vi)(−1)|J|^ϕ(ΛJ∪Li),γline i=∞∑j=0^ϕj(ΛLi)=11−^ϕ(ΛLi). (34)

Since we assumed that , we do not need to consider diagrams which contain vertices with only one leaving line. Therefore a single diagram remains for which drawn in Fig. 2a. This gives

 ^C2(λ1,λ2)=γline,L={1,2}γvertex,V={1,2},L={}=μ2λ1λ21−^ϕ(λ1)−^ϕ(λ2)+^ϕ(λ1+λ2)1−^ϕ(λ1+λ2). (35)

For the four point correlation two types of diagrams are relevant, which are drawn in Fig. 2b. The contributions are

 DI(λ)=γline,L={1,2,3,4}γvertex,V={3,4},L={1,2}=×γline,L={1,2}γ% vertex,V={1,2},L={}=μ22λ1λ2λ3λ4^ϕ(Λ{1,2})−^ϕ(Λ{1,2,3})−^ϕ(Λ{1,2,4})+^ϕ(Λ{1,2,3,4})1−^ϕ(Λ{1,2,3,4})=×1−^ϕ(λ1)−^ϕ(λ2)+^ϕ(λ1+λ2)1−^ϕ(λ1+λ2) (36)
 DII(λ)=γline,L={1,2,3,4}γvertex,V={1,2,3,4},L={}=μ4λ1λ2λ3λ4∑J∈P({1,2,3,4})(−1)|J|^ϕ(ΛJ)1−^ϕ(Λ{1,2,3,4}). (37)

The function is obtained by summing over all possible combinations of indices at the vertices, i.e.,

 ^C4(λ)=DI(λ1,λ2,λ3,λ4)+DI(λ1,λ3,λ2,λ4)+DI(λ1,λ4,λ2,λ3)+DI(λ2,λ3,λ1,λ4)+DI(λ2,λ4,λ1,λ3)+DI(λ3,λ4,λ1,λ2)+DII(λ1,λ2,λ3,λ4). (38)

## 4 Evaluating the spectral observables

### 4.1 The unbinned observables

We can use the results of the last section to determine the spectral observables. Using the expansion Eq. (5) and the two point correlations Eq. (35) to calculate the limit Eq. (16)

 limζ→0+ζαλ1λ2^C2(ζλ1−iω,ζλ2+iω)=1λ1λ2(λ1+λ2)αμ2να2−^ϕ(iω)−^ϕ(−iω)ω2 (39)

The Laplace transform can be inverted by using

 L[min(T1,T2)αΓ(1+α)]=1λ1λ2(λ1+λ2)α. (40)

Therefore, one has

 ⟨Sub(ω)⟩=μ2ναΓ(1+α)2−^ϕ(iω)−^ϕ(−iω)ω2 (41)

with the behavior near

 ⟨Sub(ω)⟩=2μ2cos(π2α)ναΓ(1+α)1|ω|2−α+o(1|ω|2−α). (42)

This last equation shows that the spectrum has a form near the origin.

Following (Feller71, , Lemma XV.1.3), there exists an such that

 ^ϕ(iω)=1⇔ω∈ΩZ. (43)

The case appears if and only if the waiting times can take only the values . Therefore, more generally,

 ^ϕ(z)=^ϕ(z+iω)for Re(z)≥0 and ω∈ΩZ. (44)

In this paper we concentrate on the spectrum for frequencies . The frequencies behave differently, e.g., the behavior of is completely described by the result of Rebenshtok and Barkai RebenshtokBarkai07 (); RebenshtokBarkai08 (). Now, we can proceed to calculate the limit Eq. (17) by looking at the terms Eqs. (36) and (37) with the parameters given in Eq. (38). We have

 limζ→0+ζ2α(DI(ζλ1−iω1,ζλ2+iω1,ζλ3−iω2,ζλ4+iω2)+DI(ζλ3−iω2,ζλ4+iω2,ζλ1−iω1,ζλ2+iω1))=Γ(1+α)21Λα{1,2,3,4}(1Λα{1,2}+1Λα{3,4})⟨Sub(ω1)⟩⟨Sub(ω2)⟩. (45)

Using the fact , yields similarly

 limζ→0+ζ2α(DI(ζλ1−iω1,ζλ4+iω2,ζλ2+iω1,ζλ3−iω2)+DI(ζλ2+iω1,ζλ3−iω2,ζλ1−iω1,ζλ4+iω2))=Γ(1+α)21Λα{1,2,3,4}(1Λα{2,3}+1Λα{1,4})⟨Sub(ω1)⟩⟨Sub(ω2)⟩=×1ΩZ(ω1−ω2) (46)

and

 limζ→0+ζ2α(DI(ζλ1−iω1,ζλ3−iω2,ζλ2+iω1,ζλ4+iω2)+DI(ζλ2+iω1,ζλ4+iω2,ζλ1−iω1,ζλ3−iω2))=Γ(1+α)21Λα{1,2,3,4}(1Λα{1,3}+1Λα{2,4})⟨Sub(ω1)⟩⟨Sub(ω2)⟩=×1ΩZ(ω1+ω2). (47)

Finally

 limζ→0+ζ2αDII(ζλ1−iω1,ζλ4+iω2,ζλ2+iω1,ζλ3−iω2)=0. (48)

In A, we derive the following Laplace transform

 L[min(T1,T2,T3,T4)αmin(T1,T2)α×F(α,−α;1+α;min(T1,T2,T3,T4)min(T1,T2))]=Γ(1+α)2λ1λ2λ3λ41Λα{1,2,3,4}1Λα{1,2} (49)

with the hypergeometric function . As

 F(α,−α;1+α;1)=Γ(1+α)2Γ(1+2α) (50)

we obtain for the limit Eq. (11)

 ⟨Sub(ω1)Sub(ω2)⟩=2Γ(1+α)2Γ(1+2α)⟨Sub(ω1)⟩⟨Sub(ω2)⟩×(1+1ΩZ(ω1−ω2)+1ΩZ(ω1+ω2)). (51)

Eq. (51) will play an important role in determining the properties of the binned case. In particular, the last equation determines the variance of (for )

 Var[Sub(ω)]=Vub(α)⟨Sub(ω)⟩2 (52)

with the function

 Vub(α)=4Γ(1+α)2Γ(1+2α)−1. (53)

As already mentioned in the introduction, even in the ergodic limit , the value is larger than zero and therefore the observable does fluctuate.

The correlation coefficient between and (for and ; here denotes the periodicity of , see Eq. (43))

 ρ[Sub(ω1),Sub(ω2)]={1%forω1−ω2∈ΩZR(α)for ω1−ω2∉ΩZ (54)

with

 R(α)=2Γ(1+α)2−Γ(1+2α)4Γ(1+α)2−Γ(1+2α). (55)

In many cases we have , i.e.,

 ρ[Sub(ω1),Sub(ω2)]=R(α)% for |ω1|≠|ω2|. (56)

In the ergodic limit, we have . Therefore, the observations and of the spectrum at two different frequencies and are uncorrelated.

### 4.2 The binned observables

Here, we are going to use the results from the unbinned case to evaluate the limits Eqs. (22) and (24). The limits on the right hand sides of these equations can be commuted with the frequency integrals describing the binning. As this statement is not obvious, we are giving an argument.

First, we are considering with . Then we have the simple estimate

 |^C2(ζλ1−iω,ζλ2+iω)|≤1ω241−^ϕ(ζ(λ1+λ2)). (57)

The commutativity follows by dominated convergence. Therefore, we get

 ⟨Sbin(ω,Δω)⟩=1Δω∫ω+12Δωω−12Δωdω′⟨Sub(ω′)⟩ (58)

and accordingly

 ⟨Sbin(ω)⟩=⟨Sub(ω)⟩. (59)

It is not surprising that the binned estimate yields also the correct expectation value for an ensemble of realizations.

The argument for (with and ) needs more work. First notice that we have the following inequality

 |1−^ϕ(λ+i