Scaled Brownian motion as a mean field model for continuous time random walks
SScaled Brownian motion as a mean field model for continuous time random walks
Felix Thiel ∗ and Igor M. Sokolov † Institut f¨ur Physik, Humboldt-Universit¨at zu Berlin, Newtonstrasse 15, D–12489 Berlin, Germany (Dated: September 16, 2018)We consider scaled Brownian motion (sBm), a random process described by a diffusion equationwith explicitly time-dependent diffusion coefficient D ( t ) = D t α − (Batchelor’s equation) which,for α <
1, is often used for fitting experimental data for subdiffusion of unclear genesis. We showthat this process is a close relative of subdiffusive continuous-time random walks and describes themotion of the center of mass of a cloud of independent walkers. It shares with subdiffusive CTRWits non-stationary and non-ergodic properties. The non-ergodicity of sBm does not however gohand in hand with strong difference between its different realizations: its heterogeneity (“ergodicitybreaking”) parameter tends to zero for long trajectories.
PACS numbers: 05.40.Fb,05.10.Gg
Anomalous diffusion is a generic name for a class oftransport processes which are close to diffusion in theirorigin (i.e. can be represented via generalized randomwalk schemes or Langevin equations) but do not lead tothe mean squared displacement growing as the first powerof time (cid:104) r ( t ) (cid:105) = 2 dDt (1)(with D being the diffusion coefficient and d the dimen-sion of space), as predicted by the Fick’s laws. Within therandom walk schemes such deviations from the normaldiffusion picture can arise either due to broad distribu-tions of the waiting times between the steps (continuoustime random walk models, CTRW), or due to slow decayof correlations between steps, or both, see [1] for a re-view, leading to the change of the power law in the timedependence of the mean squared displacement, (cid:104) r ( t ) (cid:105) ∝ t α . The processes with α < α > D in Eq.(1) vanishesin the long time limit; in the second case it diverges.The single trajectory dynamics in normal diffusion isdescribed by the Langevin equation˙ x = √ Dξ ( t )with white, delta-correlated Gaussian noise ξ ( t ), (cid:104) ξ ( t ) (cid:105) =0, (cid:104) ξ ( t ) ξ ( t (cid:48) ) (cid:105) = δ ( t − t (cid:48) ); the time-dependence of the prob-ability density function (PDF) of the process or the oneof its transition probabilities is given by the Fick’s secondlaw (diffusion equation) ∂∂t p ( x, t ) = D ∂ ∂x p ( x, t )(both equations given here in one dimension). The de-scription of anomalous diffusion of different origins oftenfollows by modification of one of the equations above. In experiments, many processes of anomalous diffusionof unknown origin, i.e. when the observable of interestwhich cannot be fitted to the solutions of Eq.(1), are fit-ted to the results obtained for the so-called scaled Brow-nian motion (sBm) [2], a diffusion process with explicitlytime-dependent diffusion coefficient D ( t ) = D t α − . Nu-merical simulations of Ref. [3] show that, at least for thecase of FRAP (fluorescence recovery after photobleach-ing), the fits may be astonishingly good, independenton the true nature of the simulated process (percolation,CTRW, etc.). This means that the form of FRAP recov-ery curves, if they hint onto anomalous diffusion, hardlydepends on the origin of the corresponding anomaly. Us-ing sBm model for calculating other properties may how-ever be dangerous, as long as the nature of the sBm modelitself and the one of the process under investigation arenot well-understood [1].The PDF and the transition probabilities in sBm aregiven by the Batchelor’s equation [4], ∂∂t p ( x, t ) = D t α − ∂ ∂x p ( x, t ) , (2)leading to the mean squared displacement (cid:104) x ( t ) (cid:105) ∝ t α (the original one was for α = 3). Initially, the Eq.(2)was proposed for description of superdiffusive turbulentdispersion, as an alternative to the Richardson’s diffu-sion equation with the distance-dependent diffusion co-efficient. Some of its shortcomings for description of theturbulent dispersion were clear to Batchelor himself, seeRef. [5] for more detailed discussion. The equations ofBatchelor’s type are typically postulated and do not fol-low from any explicit physical model, one of seldom ex-clusions being [6].The Langevin description of the corresponding processis given by ˙ x = (cid:112) D αt α − ξ ( t ) (3)with white, delta-correlated Gaussian noise ξ ( t ), (cid:104) ξ ( t ) (cid:105) =0, (cid:104) ξ ( t ) ξ ( t (cid:48) ) (cid:105) = δ ( t − t (cid:48) ). Both, the Langevin and the dif-fusion equation for sBm can be reduced to the ones for the a r X i v : . [ phy s i c s . d a t a - a n ] N ov normal one by the time rescaling t → τ = t α . Thus, sBmis a random process which is subordinated to a Brown-ian motion (Wiener process) under the deterministic timechange given above and as such is a relative of the CTRWwith the only (but important) difference that in CTRWthe time change is stochastic. The qualitative discus-sion of the relation between the sBm and the CTRW wasgiven in [1] without proofs and calculations. Here weclose this gap and provide the deeper analysis of sBm.Thus, we show that sBm can be considered as a homog-enized (mean field) approximation to CTRW and sharesits property of aging and ergodicity breaking, but thisergodicity breaking does not go hand in hand with inho-mogeneity (non-convergence in distribution of the meansquared displacement as obtained by the moving time av-erage). This stresses that the absence of ergodicity (dueto non-stationarity) does not imply the non-zero valueof the “ergodicity breaking parameter” which character-izes the heterogeneity of realizations. In sBm, althougha close relative of CTRW, this heterogeneity is removedby the pre-averaging procedure. In what follows we con-centrate on the subdiffusive case α <
1, as considered in[3].Let us first discuss a general situation and considera random process x i ( t ) in continuous time t . Let usassume that the process possesses all necessary single-point and cross-moments. We will associate this ran-dom process with the coordinate of the i -th walker attime t . The process will be taken to possess zero mean.Let us now consider the behavior of the center of massof m independent walkers (i.e. the mean coordinate of m independent random processes x i ( t ), i = 1 , , ..., m ) X ( t ) = (1 /m ) (cid:80) mi =1 x i ( t ). Note that (cid:104) x i ( t ) (cid:105) = 0 dueto the symmetry of the process, (cid:104) x i ( t ) x j ( t ) (cid:105) = 0 (for i (cid:54) = j ) due to independence of different realizations,and (cid:104) x i ( t ) x i ( t ) (cid:105) = σ ( t ), so that (cid:104) X ( t ) (cid:105) = σ ( t ) /m .In order to prove the corresponding limit theorems weneed to define the mean position in such a way that itssecond moment does not depend on m i.e. to rescale x ( m ) ( t ) → X ( t ) √ m , or, in other words, to redefine the“center of mass” position as x ( m ) ( t ) = 1 √ m m (cid:88) i =1 x i ( t ) . Let us note that the random process x ( m ) ( t ) giving therescaled center of mass position retains the correlationfunction of the initial x i ( t ) process. Note that for sym-metric processes with correlation function (cid:104) x i ( t ) x j ( t (cid:48) ) (cid:105) = C ( t, t (cid:48) ) δ ij , x ( m ) ( t ) has the same correlation func-tion: (cid:104) x ( m ) ( t ) x ( m ) ( t (cid:48) ) (cid:105) = m (cid:80) mi,j =1 (cid:104) x i ( t ) x j ( t (cid:48) ) (cid:105) = m (cid:80) mi,j =1 C ( t, t (cid:48) ) δ ij = m (cid:80) mi =1 C ( t, t (cid:48) ) = C ( t, t (cid:48) ). Themean squared change in x ( m ) , given by C ( t, t ), doesnot depend on m . In the limiting case, x ( t ) =lim m →∞ x ( m ) ( t ) will correspond to the “mean field co-ordinate” (MF-position) of the walker. The position of the center of mass of the cloud is obtained by inverserescaling.We now show that the position x ( t ) and all possiblevectors ( x ( t ) , x ( t ) , ..., x ( t n )) comprising the walkers po-sitions at different times converge to Gaussian, so thatthe process x ( m ) ( t ) converges to a Gaussian process for m → ∞ .Let us consider a sequence of times t , t ,..., t n anda corresponding vector X i = ( x i ( t ) , x i ( t ) , ..., x i ( t n )).This vector (stemming from probing the position of the i -th walker at time instants t , t , ..., t n ) does possess amean (cid:104) X i (cid:105) = 0 and a covariation matrix between its com-ponents with finite elements C jk = C ( t j , t k ). Accordingto the central limit theorem for multivariate (vector) dis-tributions the corresponding mean X ( m ) = √ m (cid:80) mi =1 X i converges in distribution to a multivariate Gaussian withzero mean and the same covariation matrix, see e.g. [7, 8].Since this happens for any sequence of observation times,the limit process, the one describing the “center-of-massmotion” of a cloud of walkers, is a Gaussian process withthe correlation function inherited from the single realiza-tion x i ( t ). This process can be considered as a kind of“mean-field approximation” for our initial process x i ( t ).We have seen that pooling (superimposing many sta-tistical copies of the initial process) leads us to a Gaus-sian process with the correlation function inherited froma single copy. Now we turn to x i ( t ) being a CTRWwith the power-law waiting time density ψ ( t ) (cid:39) τ α t − − α and the mean squared displacement σ ( t ) (cid:39) D t α with D being the combination of mean squared displace-ment per step a and typical waiting time τ [9]. Letus subdivide the time axis in short intervals of dura-tion dt , and resample each of the CTRW processes asa simple random walk with the step duration dt andwith three possible step lengths s i of zero and ± i = t/dt ; the double steps during the dt intervals canbe neglected provided dt is small enough). The steps ofthis process are not independent (if the leading processis not a Poissonian one, i.e. if the waiting time distri-butions are not exponential), but always uncorrelated,because of the symmetry: (cid:104) s i s j (cid:105) = 0 for i (cid:54) = j i.e. forthe resampled process (cid:104) s i s j (cid:105) = δ ij . The displacementsof a walker during two non-intersecting time intervals∆ t = t − t , ∆ X = (cid:80) [ t /dt ] i =[ t /dt ] s i , and ∆ t = t − t ,∆ X = (cid:80) [ t /dt ] j =[ t /dt ] s j (the square brackets here denotethe whole part of the corresponding number) are non-correlated since, for all possible i and j , i (cid:54) = j . Thus,for two non-intersecting time intervals (cid:104) ∆ X ∆ X (cid:105) = 0.This observation allows us to get the position-positioncorrelation function C ( t, t (cid:48) ) = (cid:104) x ( t ) x ( t (cid:48) ) (cid:105) in CTRW. Tak-ing t (cid:48) > t on can put C ( t, t (cid:48) ) = (cid:104) x ( t )[ x ( t ) + ∆ x ( t (cid:48) − t ) (cid:105) = (cid:104) x ( t ) (cid:105) + (cid:104) ∆ x ( t (cid:48) − t ) x ( t ) (cid:105) with ∆ x ( t (cid:48) − t ) being the walker’sdisplacement during the time interval of duration t (cid:48) − t starting at t . The mean (cid:104) ∆ x ( t (cid:48) − t ) x ( t ) (cid:105) vanishes as dis-cussed above. Therefore the correlation functions C ( t, t (cid:48) )for CTRW are given by C ( t, t (cid:48) ) = (cid:104) x (min( t, t (cid:48) )) (cid:105) = 2 D [min( t, t (cid:48) )] α . (4)We now consider the superposition of m independentCTRWs. When the number m of independent pooled (su-perimposed) random processes tends to infinity, so doesalso the number of events within each interval of a fixedlength, and the displacement ∆ x of a pooled process dur-ing this interval tends to a Gaussian. The displacementsat different ∆ t -intervals are non-correlated Gaussian ran-dom variables and are therefore independent [7]: the de-pendence present in the values of x i ( t ) gets “dissolved”when the number of processes tends to infinity. The meansquared displacement during the ∆ t -interval starting at t is given by (cid:104) ∆ x (cid:105) = 2 D ( t + ∆ t ) α − D t α (cid:39) D αt α − ∆ t. (5)We then can write x ( t + ∆ t ) = x ( t ) + ∆ x ( t ) and inter-pret it as a finite-difference approximation to a Langevinequation dxdt = η ( t ) , i.e. as the result of integrating this equation over the fi-nite time interval ∆ t starting at t . Since the distributionof ∆ x ( t ) is Gaussian, the values of ∆ x at different timeintervals are uncorrelated, and (cid:104) ∆ x (cid:105) = 2 D αt α − ∆ t ,the properties of the noise follow: this noise is Gaussianwith (cid:104) η ( t ) (cid:105) = 0 and (cid:104) η ( t ) η ( t (cid:48) ) (cid:105) = 2 D αt α − δ ( t − t (cid:48) ), i.e.exactly as in Eq.(3). This diffusion process with explic-itly time-dependent diffusion coefficient is the sBm, andthe Batchelor’s equation, Eq.(2) appears as a Fokker-Planck equation for the Langevin equation above. Themean field process shares with the CTRW the same time-dependence of the mean squared displacement; (cid:104) x ( t ) (cid:105) ∝ t α in the ensemble average. Therefore, the double, en-semble and moving time average, goes as (cid:104) x ( t ) (cid:105) = 1 T − t (cid:90) T − t (cid:104) [ x ( t + τ ) − x ( τ )] (cid:105) dτ = 2 D T − t (cid:90) T − t [( t + τ ) α − τ α ] dτ ≈ D T α − t for t (cid:28) T , just like in CTRW [10], showing ergodicitybreaking (a trivial one, due to the non-stationarity ofthe process). At difference with CTRW, the “ergodicitybreaking parameter” for this process vanishes, showingthat its different realizations are extremely similar.Let x ( t ) = T − t (cid:82) T − t [ x ( τ + t ) − x ( τ )] dτ be the mov-ing time averaged mean squared displacement. In CTRWthis one is a random variable: its distribution does notnarrow for T → in the sense that the heterogeneity (“er-godicity breaking”) parameter EB [11] EB = (cid:68) x ( t ) (cid:69) − (cid:104) x ( t ) (cid:105) (cid:104) x ( t ) (cid:105) (6) tends to a finite limit EB = 2Γ (1 + α ) / Γ(1 + 2 α ) − T → ∞ . Expressing EB via correlation functionswe see that vanishing of this parameter for a stationary process leads to ergodicity (i.e. equality of the ensemblemean and the moving time average) as a consequence ofBirkhoff-Khinchin theorem, stressing that ergodicity im-plies stationarity and some amount of “homogeneity”.Scaled Brownian motion (Batchelor’s process) is non-stationary (since it is explicitly time-inhomogeneous),and non-ergodic, but shows high amount of homogene-ity among its trajectories, as we proceed to show by ex-plicit calculation of EB and showing that it vanishes for T → ∞ .The only thing we need to calculate is (cid:104) x ( t ) (cid:105) : (cid:68) x ( t ) (cid:69) = 1( T − t ) (cid:42)(cid:90) T − t dτ [ x ( τ + t ) − x ( τ )] × (cid:90) T − t dτ [ x ( τ + t ) − x ( τ )] (cid:43) = 2( T − t ) (cid:42)(cid:90) T − t dτ d (cid:90) T − t − τ d dτ m [ x ( τ m + t ) − x ( τ m )] [ x ( τ m + τ d + t ) − x ( τ m + τ d )] (cid:29) , (7)where we changed the variables of integration to τ m =min( τ , τ ), τ d = | τ − τ | and used the symmetry ofthe expression. The order of integration and ensem-ble averaging can be interchanged and the mathemati-cal expectation is readily evaluated by using the Gaus-sian property of x ( t ). Let us denote by ˜ C ( τ , τ ) = (cid:104) [ x ( τ + t ) − x ( τ )][ x ( τ + t ) − x ( τ )] (cid:105) the correlation func-tion of the increments. Expressed in variables of Eq.(7)˜ C reads: ˜ C ( τ m , τ m + τ d ) = 2 D [( τ m + t ) α − (min( τ m + t, τ m + τ d )) α ] = 2 D [( τ m + t ) α − ( τ m + τ d ) α ] I [0 ,t ) ( τ d ),where I A ( x ) is an indicator function that equals unity if x ∈ A and vanishes otherwise. Note, that the incrementprocess is also Gaussian, so that its higher correlators de-compose into products of ˜ C ( τ , τ ). Thus the integrandin the last expression of Eq.(7) reads: (cid:68) [ x ( τ m + t ) − x ( τ m )] [ x ( τ m + τ d + t ) − x ( τ m + τ d )] (cid:69) = 2 ˜ C ( τ m , τ m + τ d ) + ˜ C ( τ m , τ m ) ˜ C ( τ m + τ d , τ m + τ d )= 8 D [( τ m + t ) α − ( τ m + τ d ) α ] I [0 ,t ) ( τ d )+ ˜ C ( τ m , τ m ) ˜ C ( τ m + τ d , τ m + τ d )The last summand can be identified with the (cid:104) x ( t ) (cid:105) term in the numerator of Eq. (6) and therefore cancelsout. The remaining summand vanishes, if τ d > t . Thisis the main difference to usual CTRW, where the exis-tence of such longer correlations is responsible for non-vanishing EB-parameter.Now let us assume T to be sufficiently large, i.e. T − t (cid:29) t , and proceed. (cid:68) x ( t ) (cid:69) − (cid:104) x ( t ) (cid:105) = 16 D ( T − t ) × (cid:90) t dτ d (cid:90) T − t − τ d dτ m [( τ m + t ) α − ( τ m + τ d ) α ] . (8)For large τ m the integrand is asymptotically equal to α τ α − m ( t − τ d ) . For α > / T (cid:29) t is: (cid:68) x ( t ) (cid:69) − (cid:104) x ( t ) (cid:105) = 16 α D t α − T − t ) − α . Inserting the corresponding expression into Eq.(6) showsthat EB approaches zero as T − .For α ≤ / τ m leads to the diver-gence of the inner integral in Eq.(8) at the lower limitof integration. Let us split the integral into two partsat some intermediate time t i > (cid:82) t dτ d (cid:82) T − t − τ d dτ m [( τ m + t ) α − ( τ m + τ d ) α ] = (cid:82) t i dτ m [( τ m + t ) α − ( τ m + τ d ) α ] + α (cid:82) T − t − τ d t i dτ m ( t − τ d ) τ α − m . The first summand is bounded from above by t i [ t α − τ αd ] , and the second one converges on the upperlimit, which thus can be taken to go to infinity. Thereforethe inner integral (and thus the whole double integral inEq.(8)) tends for T large to a function independent on T . We thus conclude that (cid:68) x ( t ) (cid:69) − (cid:104) x ( t ) (cid:105) (cid:104) x ( t ) (cid:105) (cid:39) Z α (cid:18) tT (cid:19) α , α ≤ α α − tT , α >
12 (9)with Z α = (cid:82) dy (cid:82) ∞ dx [( x + 1) α − ( x + y ) α ] .We see, that for any positive value of α EB willvanish, and shows a crossover between two types of T -dependence at α = 1 /
2. Figure 1 shows simulated valuesof EB for different values of α and T in a double loga-rithmic plot together with the predictions of Eq(9).In conclusion: Scaled Brownian motion (Batchelor’sprocess), described by a diffusion equation with explic-itly time-dependent diffusion coefficient and often used infitting data for experimental situations showing anoma-lous diffusion of unclear origin, is, in subdiffusive case α <
1, a close relative of CTRW and describes the motionof the center of mass of a cloud independent continuoustime random walkers. The model shows the same kindof ergodicity breaking due to non-stationarity as the cor-responding CTRW process. This one, however, does not go hand in hand with strong heterogeneity of its differ-ent realizations: the corresponding “ergodicity breaking”parameter vanishes in the limit of long trajectories. -5 -4 -3 -2 -1 H e t e r og e n e i t y P a r a m e t e r EB Observation Time T Scaling Behaviour
FIG. 1: Scaling Behaviour of EB . The Batchelor’s Processwas simulated and EB was calculated for different values of α and T . α -values range from 0 . . × ), 0 . (cid:3) ), 0 . (cid:12) ),and 1 . (cid:52) ). Each point is calculated from 5000 trajectories.The given power laws are (from upper to lower line ) T − . , T − . , and T − . The authors acknowledge financial support by DFGwithin IRTG 1740 research and training group project. ∗ Electronic address: [email protected] † Electronic address: [email protected][1] I.M. Sokolov, Soft Matter, , 9043 (2012)[2] S. C. Lim and S. V. Muniandy, Phys. Rev. E , 021114(2002)[3] M. J. Saxton, Biophys. J. , 2226 (2001)[4] G.K. Batchelor, Math. Proc. Cambridge Phil. Soc., ,345 (1952)[5] A.S. Monin, A.M. Yaglom, Statistical Fluid Mechanics:Mechanics of Turbulence , Volume 1, M.I.T. Press, Cam-bridge, Mass. (1971)[6] E.B. Postnikov and I.M. Sokolov, Physica A , 5095 -5101 (2012)[7] W. Feller,
An Introduction to Probability Theory and ItsApplications , Wiley, ... Vol. II Chapter III.[8] P. Mukhopadhyay,
Multivariate Statistical Analysis ,World Scientific, Singapore, 2009 (see Sec. 3.5).[9] J. Klafter and I.M. Sokolov,
First Steps in RandomWalks: From Tools to Applications , Oxford universityPress (2011)[10] A. Lubelski, I.M. Sokolov and J. Klafter, Phys. Rev. Lett. , 250602 (2008)[11] Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev.Lett.101