Denis S Grebenkov

Time-averaged MSD for switching diffusion

Keywords: numbers: 02.50.-r, 05.40.-a, 02.70.Rr, 05.10.Gg switching diffusion, Kärger model, ergodicity breaking, time averaged MSD

We consider a classic two-state switching diffusion model from a single-particle tracking perspective. The mean and the variance of the time-averaged mean square displacement (TAMSD) are computed exactly. When the measurement time (i.e., the trajectory duration) is comparable to or smaller than the mean residence times in each state, the ergodicity breaking parameter is shown to take arbitrarily large values, suggesting an apparent weak ergodicity breaking for this ergodic model. In this regime, individual random trajectories are not representative while the related TAMSD curves exhibit a broad spread, in agreement with experimental observations in living cells and complex fluids. Switching diffusions can thus present, in some cases, an ergodic alternative to commonly used and inherently non-ergodic continuous-time random walks that capture similar features.

I. INTRODUCTION

A continuous-time random walk (CTRW) was originally proposed to model hopping processes in semiconductors with randomly distributed heavy tailed waiting times due to random energetic trapping [START_REF] Montroll | Random Walks on Lattices. II[END_REF][START_REF] Bouchaud | Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications[END_REF]. This model has been extensively studied and rapidly became an archetypical model of anomalous diffusion in different contexts, including the intracellular transport in microbiology [START_REF] Metzler | The random walk's guide to anomalous diffusion: a fractional dynamics approach[END_REF][START_REF] Bressloff | Stochastic models of intracellular transport[END_REF][START_REF] Metzler | Anomalous diffusion models and their properties: nonstationarity, non-ergodicity, and ageing at the centenary of single particle tracking[END_REF]. This is also an emblematic model of weak ergodicity breaking (WEB) and aging phenomena [START_REF] Bouchaud | Weak ergodicity breaking and aging in disordered systems[END_REF][START_REF] Burov | Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking[END_REF][START_REF] Schulz | Aging Effects and Population Splitting in Single-Particle Trajectory Averages[END_REF]. In fact, when the waiting time distribution has no first moment, there is no characteristic time scale for waiting events, so that the duration of a single stalling period may be comparable to the overall measurement time that prohibits self-averaging needed for ergodic behavior. Moreover, longer measurement times favor longer waiting times so that the statistical properties of the system depend on the measurement time. In the biological context, CTRWs were applied to mimic molecular caging in an overcrowded intracellular environment when a diffusing macromolecule, surrounded by other macromolecules and filaments, may wait a relatively long time in such an effective molecular "cage" before being able to jump to a next cage [START_REF] Barkai | Strange kinetics of single molecules in living cells[END_REF][START_REF] Wong | Anomalous Diffusion Probes Microstructure Dynamics of Entangled F-Actin Networks[END_REF][START_REF] Jeon | In Vivo Anomalous Diffusion and Weak Ergodicity Breaking of Lipid Granules[END_REF][START_REF] Weigel | Ergodic and nonergodic processes coexist in the plasma membrane as observed by single-molecule tracking[END_REF]. However, the validity of the considered heavy-tailed distribution of waiting times in the biological context remains debatable [START_REF] Szymanski | Elucidating the Origin of Anomalous Diffusion in Crowded Fluids[END_REF]. In particular, such a distribution may have an exponential cut-off so that anomalous features are transient, before getting into normal diffusion at long times [START_REF] Burov | Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking[END_REF][START_REF] Lanoiselée | Revealing nonergodic dynamics in living cells from a single particle trajectory[END_REF]. Although a CTRW with cut-off sounds more realistic, the underlying mathematical formalism is rather difficult.

In this light, a simpler model of random switching between several diffusion states, characterized by diffusion * Electronic address: denis.grebenkov@polytechnique.edu coefficient D 1 , . . . , D J , is appealing. A particle started in a state i undergoes normal diffusion with diffusivity D i for a random time, until it switches to another state j, with the switching rate k ij , and so on. Diffusion states can represent different conformations of a large macromolecule (e.g., globular versus filamentous structures) and thus distinct effective radii. Alternatively, diffusion states can account for eventual temporal binding of the diffusing particle to other molecules that may either slow down or even stop its motion. Such switching diffusions are often employed to describe the dynamics in biological systems [START_REF] Bressloff | Stochastic switching in biology: from genotype to phenotype[END_REF][START_REF] Sungkaworn | Single-molecule imaging reveals receptor-G protein interactions at cell surface hot spots[END_REF][START_REF] Weron | Ergodicity breaking on the neuronal surface emerges from random switching between diffusive states[END_REF] and used as simple models of intermittent processes [START_REF] Lanoiselée | Unravelling intermittent features in single particle trajectories by a local convex hull method[END_REF].

In this paper, we mainly focus on the two-state model with diffusion coefficients D 1 and D 2 and the switching rates k 12 and k 21 . Note that 1/k 12 and 1/k 21 are the mean residence times of the particle in the states 1 and 2, respectively. The molecular caging effect can be modeled by setting D 2 = 0, i.e., the particle does not move in the second state. This model has been extensively studied, in particular, in the nuclear magnetic resonance literature, in which it is known as the Kärger model [START_REF] Kärger | NMR self-diffusion studies in heterogeneous systems[END_REF][START_REF] Fieremans | Monte Carlo study of a two-compartment exchange model of diffusion[END_REF][START_REF] Moutal | The Kärger vs bi-exponential model: theoretical insights and experimental validations[END_REF]. Here, we look at this model from a single-particle tracking perspective. This simple model will allow us to investigate the reproducibility of measurements over individual random trajectories and the effect of their duration T . In particular, if both residence times 1/k 12 and 1/k 21 are small as compared to T , the particle switches very often between two states and manages to probe these states reliably during the measurement time. In this limit, the intermittent process is seen as an ergodic normal diffusion with some mean diffusivity D (see below). In contrast, if both residence times are much larger than T , the particle remains within a single state over the whole measurement with a high probability. As a consequence, such a single trajectory would bring information only about one state, whereas another trajectory may bring information only about the other state. In other words, individual trajec-tories are not representative of the whole dynamics and thus, from a practical point of view, not ergodic. Since the switching model is ergodic, we call this regime apparent WEB due to insufficient measurement time. Finally, the situation when one or both residence times are comparable to the measurement time is a borderline case. This model will thus help us to study the transition from one regime to the other that is quite typical for intracellular transport measurements.

In order to quantify the reliability of measurements and an eventual apparent WEB, we investigate the timeaveraged mean square displacement (TAMSD) of a particle undergoing such an intermittent motion:

δ 2 (t, T ) = 1 T -t T -t 0 dt 0 (X(t + t 0 ) -X(t 0 )) 2 , ( 1
)
where t is the lag time. We compute exactly both the mean and the variance of the random variable δ 2 (t, T ) to evaluate the ergodicity breaking (EB) parameter χ, also known as the squared coefficient of variation of the TAMSD [START_REF] He | Random Time-Scale Invariant Diffusion and Transport Coefficients[END_REF][START_REF] Burov | Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking[END_REF]:

χ = [δ 2 (t, T )] 2 δ 2 (t, T ) 2 -1 = var{δ 2 (t, T )} (mean{δ 2 (t, T )}) 2 , (2)
where . . . denotes the expectation, i.e., the ensemble average over the space of all possible trajectories X(t) of duration T . The paper is organized as follows. In Sec. II, we provide the propagator for multiple successive points in a matrix form which is then used to derive the mean and the variance of the TAMSD. Different asymptotic regimes of the resulting EB parameter are then analyzed. In Sec. III, we discuss an apparent WEB when the measurement time is comparable to or smaller than the mean residence times. In particular, we illustrate a large spread of TAMSD curves obtained via numerical simulations. We also mention extensions for multi-state switching diffusion and restricted diffusion in bounded domains. Section IV summarizes the results and concludes.

II. THEORETICAL RESULTS

A. Propagator

We consider a particle that diffuses on a real axis R and switches randomly between two diffusivities D 1 and D 2 with rates k 12 and k 21 (see Sec. III B and III C for extensions). Let P ii0 (x, t|x 0 , t 0 ) be the propagator of this particle, i.e., the probability density of finding the particle at point x in state i at time t given that it started from point x 0 in state i 0 at time t 0 . These four propagators satisfy four coupled partial differential equations

∂ ∂t P 1,i0 = D 1 ∂ 2 ∂x 2 P 1,i0 -k 12 P 1,i0 + k 21 P 2,i0 , (3a) ∂ ∂t P 2,i0 = D 2 ∂ 2 ∂x 2 P 2,i0 -k 21 P 2,i0 + k 12 P 1,i0 (3b)
(with i 0 = 1, 2), subject to the initial conditions: P i,i0 (x, t = t 0 |x 0 , t 0 ) = δ i,i0 δ(x -x 0 ) (a rigorous mathematical formulation of switching models and of some their properties can be found in [START_REF] Yin | Hybrid Switching Diffusions: Properties and Applications[END_REF][START_REF] Yin | Properties of solutions of stochastic differential equations with continuous-state-dependent switching[END_REF][START_REF] Baran | Feynman-Kac formula for switching diffusions: connections of systems of partial differential equations and stochastic differential equations[END_REF]). The solution of these equations is well known (e.g., see [START_REF] Kärger | NMR self-diffusion studies in heterogeneous systems[END_REF]). Here we recall the main formulas that will be needed for the analysis of the TAMSD. In the Fourier space,

P ii0 (x, t|x 0 , t 0 ) = ∞ -∞ dq 2π e -iq(x-x0) Pii0 (q, t -t 0 ), (4)
Eqs.

(3) are reduced to coupled first-order differential equations whose solution is easily found in a matrix form P(q, t) = P11 (q, t) P12 (q, t) P21 (q, t) P22 (q, t) = exp(-M q t), [START_REF] Metzler | Anomalous diffusion models and their properties: nonstationarity, non-ergodicity, and ageing at the centenary of single particle tracking[END_REF] where

M q = D 1 q 2 + k 12 -k 21 -k 12 D 2 q 2 + k 21 . (6)
The eigenvalues and eigenvectors of the matrix M q ,

M q V q = V q γ + q 0 0 γ - q .
are known explicitly:

γ ± q = 1 2 (D 1 + D 2 )q 2 + (k 12 + k 21 ) (7)
± ((D 2 -D 1 )q 2 + (k 21 -k 12 )) 2 + 4k 12 k 21 and V q =   k 21 k 12 + D 1 q 2 -γ + q k 21 k 12 + D 1 q 2 -γ - q 1 1   . ( 8
)
One gets thus P(q, t) = e -γ + q t µ - q -e -γ - q t µ + q k 21 (e -γ - q t -e -γ + q t ) k 12 (e -γ - q t -e -γ + q t ) e -γ - q t µ - q -e -γ + q t µ + q γ + q -γ - q , (9) with µ ± q = D 1 q 2 + k 12 -γ ± q . In the special case q = 0, one has γ + 0 = k 12 + k 21 and γ - 0 = 0. To compute the moments of the TAMSD, one needs the propagator for multiple successive points, which is simply the product of the above propagators due to the Markov property of the process:

P (x n , i n , t n ; x n-1 , i n-1 , t n-1 ; . . . ; x 1 , t 1 , i 1 ; x 0 , i 0 , 0) = P in,in-1 (x n , t n |x n-1 , t n-1 ) . . . P i1,i0 (x 1 , t 1 |x 0 , 0) = R n dq 1 2π . . . dq n 2π e -iqn(xn-xn-1)-...-iq1(x1-x0)
× Pin,in-1 (q n , t n -t n-1 ) . . . Pi1,i0 (q 1 , t 1 ). [START_REF] Wong | Anomalous Diffusion Probes Microstructure Dynamics of Entangled F-Actin Networks[END_REF] Denoting by p 1 (resp., p 2 = 1 -p 1 ) the probability of starting in the state 1 (resp. 2) at time 0, the marginal propagator averaged over the state variables i k reads

P (x n , t n ; x n-1 , t n-1 ; . . . ; x 1 , t 1 ; x 0 , 0) = R n dq 1 2π . . . dq n 2π e -iqn(xn-xn-1)-...-iq1(x1-x0) ×P n (q n , t n -t n-1 ; . . . ; q 2 , t 2 -t 1 ; q 1 , t 1 ), (11)
where

P n (q n , t n -t n-1 ; . . . ; q 2 , t 2 -t 1 ; q 1 , t 1 ) (12) = 1 1 † P(q n , t n -t n-1 ) . . . P(q 1 , t 1 ) p 1 p 2 .

B. Mean TAMSD

We first calculate the characteristic function of the displacement between times t 1 and t 2 such that 0 < t 1 < t 2 :

G(q) ≡ e iq(X(t2)-X(t1)) = ∞ -∞ dx 1 ∞ -∞ dx 2 e iq(x2-x1) P (x 2 , t 2 ; x 1 , t 1 ; x 0 , 0) = P 2 (q, t 2 -t 1 ; 0, t 1 ). ( 13
)
From this characteristic function, one can compute the moments of the displacement, in particular, the mean square displacement:

(X(t + t 0 ) -X(t 0 )) 2 = -lim q→0 ∂ 2 G(q) ∂q 2 (14) = 2 Dt + 2 (p 1 k 12 -p 2 k 21 )(D 1 -D 2 ) k 2 (1 -e -kt )e -kt0 ,
where

k = k 12 + k 21 (15)
and

D = D 1 k 21 + D 2 k 12 k (16)
is the mean diffusivity. In particular, the ensemble averaged MSD is

X 2 (t) = 2 Dt + 2 (p 1 k 12 -p 2 k 21 )(D 1 -D 2 ) k 2 (1 -e -kt ) . (17)
From Eq. ( 14), one can also deduce the mean TAMSD

δ 2 (t, T ) = 2 Dt (18) + 2(p 1 k 12 -p 2 k 21 )(D 1 -D 2 ) k 3 (T -t) (1 -e -kt )(1 -e -k(T -t) ) .
If the initial probabilities p 1 and p 2 are set to be from the equilibrium,

p 1 = p eq 1 = k 21 k , p 2 = p eq 2 = k 12 k , (19)
then

δ 2 (t, T ) = 2 Dt. ( 20
)
One cannot therefore reveal the intermittent character of this process from the mean TAMSD alone. Moreover, the mean value does not depend on the measurement time T , as for normal Brownian motion.

C. Variance of TAMSD

The computation of the variance of the TAMSD is much more involved as it requires the four-point correlation function. In fact, one has

[δ 2 ] 2 = 2 (T -t) 2 T -t 0 dt 0 T -t t0 dt ′ 0 (21) × (X(t 0 + t) -X(t 0 )) 2 (X(t ′ 0 + t) -X(t ′ 0 )) 2 ,
from which the variance follows as usual:

var{δ 2 } = [δ 2 ] 2 -δ 2 2 . ( 22
)
For computing the second moment of the TAMSD, we introduce the characteristic function

G(q, q ′ ) ≡ e iq(X(t0+t)-X(t0))+iq ′ (X(t ′ 0 +t)-X(t ′ 0 )) . (23)
Since we have t 0 < t ′ 0 in Eq. ( 21), there are two cases:

(i) for t 0 < t 0 + t < t ′ 0 < t ′ 0 + t, we get G(q, q ′ ) = P 4 (q ′ , t; 0, t ′ 0 -t 0 -t; q, t; 0, t 0 ); ( 24
) (ii) for t 0 < t ′ 0 < t 0 + t < t ′ 0 + t, we get G(q, q ′ ) = P 4 (q ′ , t ′ 0 -t 0 ; q + q ′ , t 0 + t -t ′ 0 ; q, t ′ 0 -t 0 ; 0, t 0 ). ( 25
)
In the evaluation of the integral in Eq. ( 21), we consider separately two cases: t < T /2 and t > T /2.

Case t < T /2

In this case, one can split the integral in Eq. ( 21) into three parts

[δ 2 ] 2 = 2 (T -t) 2 T -2t 0 dt 0 t 0 dt ′ 0 F 2 (t, t 0 , t ′ 0 ) (26) + T -2t 0 dt 0 T -2t-t0 0 dt ′ 0 F 1 (t, t 0 , t ′ 0 ) + t 0 dt 0 t0 0 dt ′ 0 F 2 (t, T -t -t 0 , t ′ 0 )
(note that the integration variables t 0 and t ′ 0 were shifted), with

F 1 (t, t 0 , t ′ 0 ) = lim q→0 q ′ →0 ∂ 4 P 4 (q ′ , t; 0, t ′ 0 ; q, t; 0, t 0 ) ∂q 2 ∂q ′2 , F 2 (t, t 0 , t ′ 0 ) = lim q→0 q ′ →0 ∂ 4 P 4 (q ′ , t ′ 0 ; q + q ′ , t -t ′ 0 ; q, t ′ 0 ; 0, t 0 ) ∂q 2 ∂q ′2 .
After long and cumbersome computations of these derivatives and integrals in Eq. ( 26), we derive the following expression for the variance of the TAMSD for t < T /2:

var{δ 2 } = 4 D2 t 3 (4T -5t) 3(T -t) 2 + 8k 12 k 21 (D 1 -D 2 ) 2 k 6 (T -t) 2 e -kT (1 -e kt ) 2 -2(3 + 2k(T -t))e -kt (27)
+ kT (3(kt) 2 -4kt + 4) -2(2(kt) 3 -3(kt) 2 + 5kt -3) ,
where the initial probabilities p i were set to their equilibrium values p eq i in Eq. ( 19) to get a more compact expression.

Case t > T /2

In this case, only the second option in Eq. ( 25) is possible so that one does not need to split the integral in Eq. ( 21) into three parts, and one gets

[δ 2 ] 2 = 2 (T -t) 2 T -t 0 dt 0 T -t-t0 0 dt ′ 0 F 2 (t, t 0 , t ′ 0 ), ( 28
)
where the integration variable was shifted in the second integral. The evaluation of this integral yields the variance of the TAMSD for t > T /2:

var{δ 2 } = 4 D2 (T 2 -6T t + 11t 2 ) 3 + 8k 12 k 21 (D 1 -D 2 ) 2 k 6 (T -t) 2 e -kT (1 -2e kt ) -2(3 + 2k(T -t))e -kt (29)
+ 5e k(T -2t) -k 3 T 3 -2k 2 T 2 (3kt -1) + kT (9k 2 t 2 -4kt + 2) -2(2k 3 t 3 -k 2 t 2 + kt + 1) ,
where the initial probabilities p i were again set to their equilibrium values p eq i for simplicity.

D. Ergodicity breaking parameter

The expressions [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF] for the variance of the TAMSD present the main computational result of this paper. The first term in these expressions is the variance of the TAMSD for Brownian motion with the mean diffusivity D [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF][START_REF] Grebenkov | Probability distribution of the timeaveraged mean-square displacement of a Gaussian process[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF][START_REF] Grebenkov | Optimal and suboptimal quadratic forms for noncentered Gaussian processes[END_REF]. One can check that the second term in Eqs. [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF] is nonnegative, i.e., the switching between two states can only increase the variance of the TAMSD.

As Eqs. [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF] for the variance are provided exclusively for the case p i = p eq i , for which the mean TAMSD in Eq. ( 20) is particularly simple, the analysis of the variance is equivalent, up to a simple multiplicative factor (2 Dt) 2 , to that of the ergodicity breaking parameter χ defined by Eq. ( 2). In the rest of the paper, we focus on this parameter.

As the trajectory duration T goes to infinity (for a fixed lag time t), the EB parameter vanishes asymptotically as

χ ≃ T -1 4t 3 + 2k 12 k 21 (D 1 -D 2 ) 2 k 3 D2 (30) × 3(kt) 2 -4kt + 4 -4e -kt (kt) 2 + O(T -2 ),
so that the switching process is ergodic, as expected. In contrast, in the double limit k 21 = p eq 1 k → 0 and k 12 = p eq 2 k → 0 (with fixed p eq i and k → 0), Eq. ( 27) yields

χ ≃ p eq 1 p eq 2 (D 1 -D 2 ) 2 (p eq 1 D 1 + p eq 2 D 2 ) 2 + O(T -1 ). (31)
In this limit, the particle stays infinitely long in either of two states, i.e., the process is not ergodic, and the variance of the TAMSD does not vanish in the limit T → ∞. This singular situation can also describe two populations of particles with distinct diffusivities D 1 and D 2 , and the leading term of the EB parameter in Eq. ( 31) is a consequence of their mixture with relative fractions p eq 1 and p eq 2 . This limit could alternatively be obtained by setting

δ 2 = α ζ 1 + (1 -α) ζ 2
, where ζ i is the TAMSD for the i-th population (that differs by the factor D i ), and a random selection between two populations is realized by a Bernoulli random variable α taking the value 1 with probability p eq 1 and 0 with probability p eq 2 = 1 -p eq 1 . Moreover, the EB parameter χ remains close to the limiting expression (31) when T ≪ 1/k. In other words, if the measurement time T is short as compared to the residence time 1/k, the system exhibits an apparent WEB. Clearly, the order of two limits, k 12 , k 21 → 0 and T → ∞, does matter here: sending T → ∞ for fixed k 12 and k 21 yields the zero variance, as expected.

In the limit t → 0, Eq. ( 27) gives

χ ≃ 2k 12 k 21 (D 1 -D 2 ) 2 (kT -1 + e -kT ) k 4 D2 T 2 + O(t). ( 32
)
We first consider the particular case, in which two switching rates are identical: k 12 = k 21 = k/2. For our illustrative purposes, we use dimensionless units for all parameters. We set the measurement time T = 1000 (i.e., the trajectory with a thousand steps). We recall that the EB parameter for Brownian motion is a monotonously growing function of the lag time t, so that the smallest available lag time t = 1 provides the most accurate estimation of the TAMSD (see, e.g., [START_REF] Grebenkov | Probability distribution of the timeaveraged mean-square displacement of a Gaussian process[END_REF]). One can check that the same property holds for two-state switching diffusion. As a consequence, we select the lag time t = 1 to be in the optimal situation. Figure 1 shows the behavior of the EB parameter as a function of the residence time 1/k. In the case of equal diffusion coefficients, D 2 /D 1 = 1, the states are identical, and the switching model is reduced to normal diffusion. The EB parameter does not depend on the switching rate and is equal (in the leading order) to 4t/(3T ) according to Eq. [START_REF] Grebenkov | Optimal and suboptimal quadratic forms for noncentered Gaussian processes[END_REF]. In turn, the stronger the difference between two states (i.e., smaller D 2 /D 1 ), the larger the EB parameter at large residence times 1/k. In the extreme case of one immobile state (with D 2 = 0), the EB parameter reaches the value 1 according to Eq. [START_REF] Golding | Physical Nature of Bacterial Cytoplasm[END_REF]. In contrast, a fast switching even between very distinct states leads again to normal diffusion with the mean diffusivity D and thus the EB parameter remains close to 4t/(3T ). Note that in the considered case of equal switching rates, the EB parameter does not exceed 1. We conclude that, due to slow switching between states, the distribution of TAMSD even for the smallest lag time (t = 1) can be relatively broad, i.e., the standard deviation can be comparable to the mean value: χ ∼ 1. In this situation, an increase of the trajectory duration T can improve the estimation only when T exceeds the mean residence time 1/k. This is in contrast with the case of normal diffusion, for which the EB parameter is of the order of t/T and can thus be made very small, allowing for accurate estimations of the diffusion coefficient for long enough trajectories.

A richer insight on the EB parameter is provided in Fig. 2 that shows the dependence on two residence times 1/k 12 and 1/k 21 . In this situation, there are four relevant time scales: t, T , 1/k 12 and 1/k 21 and thus various regimes. As previously, we fix t = 1 and T = 1000. First, we present in Fig. 2(a) the case of an immobile state with D 2 = 0 (and D 1 = 1). When 1/k 21 is the smallest time scale, the particle does not almost stay in the immobile state, and the EB parameter is close to its value 4t/(3T ) for normal diffusion. When 1/k 21 is getting comparable to the lag time t, the switching is not fast enough any more, and the presence of the immobile phase increases the EB parameter. Finally, when 1/k 21 increases further, the EB parameter can become much larger than 1. In fact, in the limit k 21 → 0 (with fixed k 12 ), one gets

χ ≃ 2A k 3 12 (T -t) 2 t 2 k -1 21 + O(1), (33)
where The situation is considerably different for a particle undergoing slow diffusion in the second state (i.e., D 2 is small but not strictly zero). Figure 2(b) shows an example for D 2 = 0.01. While the behavior of the EB parameter for small 1/k 21 is expectedly similar to the former case with D 2 = 0, there is significant difference for large 1/k 21 . First, the values of the EB parameter, which can still be large, are much smaller than those shown in Fig. 2(a). Second, the EB parameter exhibits a maximum as a function of the mean residence time 1/k 12 for a fixed 1/k 21 . When 1/k 12 is small (with 1/k 21 large), the particle stays most of the time in the second state and thus undergoes almost normal diffusion with diffusivity D 2 , so that the EB parameter again reaches the small value 4t/(3T ). This is a significant difference with respect to the case D 2 = 0.

A = e -k12T (

III. DISCUSSION

When the measurement time T is comparable to or smaller than the residence times 1/k 12 1/k 21 , the particle does not have enough time to probe the phase space. An individual random trajectory of the particle is thus not representative of the ensemble, suggesting an apparent weak ergodicity breaking. In particular, the EB parameter can be of order of 1 and even much larger, particularly if the diffusion coefficient in one of the states is very small (or zero). As a consequence, one can expect a large spread of TAMSD curves evaluated from individual trajectories. This feature, which is often observed in experiments (see, e.g., [START_REF] Jeon | In Vivo Anomalous Diffusion and Weak Ergodicity Breaking of Lipid Granules[END_REF][START_REF] Golding | Physical Nature of Bacterial Cytoplasm[END_REF][START_REF] Mosqueira | Cholesterol modulates acetylcholine receptor diffusion by tuning confinement sojourns and nanocluster stability[END_REF]), was often interpreted as an indication of non-ergodicity and thus attributed to CTRW as a basic non-ergodic model. Here we showed that the ergodic two-state model can lead to similar features.

A. Spread of TAMSD curves

In order to illustrate the spread of TAMSD curves, we perform Monte Carlo simulations. For a prescribed set of parameters D 1 , D 2 , k 12 and k 21 , we set the trajectory duration T = 1000 (i.e., the number of steps) and generate T random centered Gaussian increments with unit variance. We also generate a sequence of successive residence times in two alternating states according to the exponential laws with the rates k 12 and k 21 , that results in a random sequence of state variables i 1 , i 2 , . . . (each i n taking values 1 or 2). All the increments are rescaled by 2D in and then cumulatively summed up to produce a random trajectory, from which the TAMSD is computed for all lag times from 1 to T -1 by discretizing Eq. ( 1). This computation is repeated 10 000 times to obtain a reliable statistics. Fixing D 1 = 1, we are left with D 2 , k 12 and k 21 as the major parameters.

Figure 3 illustrates the spread of TAMSD curves for three sets of parameters. The first set with D 2 = 1 corresponds to normal diffusion (here, switching rates do not matter as D 1 = D 2 ). In this case, TAMSD curves are close to each other at small lag times and then getting more spread at larger lag times because the time average becomes less and less efficient (Fig. 3a). In the second set (with D 2 = 0.01, k 12 = k 21 = 0.1), the mean residence times 1/k 12 and 1/k 21 are chosen to be much smaller than the measurement time T . As switching is rapid enough, TAMSD curves remain close to each (Fig. 3b), as for normal diffusion. In the third set of parameters, we keep D 2 = 0.01 but decrease both switching rates: k 12 = 10 -2 and k 21 = 10 -3 . While 1/k 12 is still much smaller than T , 1/k 21 is equal to T that leads to a wide spread of TAMSD curves (Fig. 3c), in agreement with large values of the EB parameter in this case. Similar spreads were observed in single-particle tracking experiments in living cells (see, e.g., [START_REF] Jeon | In Vivo Anomalous Diffusion and Weak Ergodicity Breaking of Lipid Granules[END_REF][START_REF] Golding | Physical Nature of Bacterial Cytoplasm[END_REF][START_REF] Mosqueira | Cholesterol modulates acetylcholine receptor diffusion by tuning confinement sojourns and nanocluster stability[END_REF]). Another way of presenting the spread consists in plotting the distribution of TAMSD for a fixed lag time. The distribution of TAMSD for ergodic Brownian motion and other Gaussian processes was studied in [START_REF] Grebenkov | Probability distribution of the timeaveraged mean-square displacement of a Gaussian process[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF][START_REF] Sikora | Statistical properties of the anomalous scaling exponent estimator based on time averaged mean square displacement[END_REF][START_REF] Gajda | Large deviations of time-averaged statistics for Gaussian processes[END_REF], while the analysis of its asymptotic behavior for nonergodic CTRW was initiated in [START_REF] He | Random Time-Scale Invariant Diffusion and Transport Coefficients[END_REF][START_REF] Lubelski | Nonergodicity Mimics Inhomogeneity in Single Particle Tracking[END_REF] (see also reviews [START_REF] Metzler | Anomalous diffusion models and their properties: nonstationarity, non-ergodicity, and ageing at the centenary of single particle tracking[END_REF][START_REF] Burov | Single particle tracking in systems showing anomalous diffusion: the role of weak ergodicity breaking[END_REF]). As our computation for two-state switching diffusion is limited to the first two moments, we show in Fig. 4 the empirical distribution of TAMSD obtained from simulated trajectories for the same three sets of parameters. The empirical distributions are presented for three lag times: t = 1, t = 10 and t = 105. As expected, the distribution is getting larger with the lag time, reflecting less and less efficient time averaging. For the second set of parameters, the distributions are close to that for normal diffusion (compare Figs. 4a and4b). In contrast, the distribution for the third set of parameters is much broader and almost does not depend on the lag time. This is the characteristic feature of non-ergodic dynamics (e.g., the distribution of TAMSD for CTRW is also broad).

To complete this discussion, we outline a subtle difference in the behavior of the EB parameter between a continuous-time process and its discrete-time approximation. Let us first discuss Brownian motion, which is often approximated by random walks on a lattice or by a sequence of Gaussian jumps. While such approximations are known to converge to Brownian motion (see, e.g., [START_REF] Mörters | Brownian Motion (Cambridge Series in Statistical and Probabilistic Mathematics[END_REF][START_REF] Grebenkov | A multiscale guide to Brownian motion[END_REF]), some functionals involving integrals along sample trajectories of these processes can be different. So, the variance of the TAMSD for continuous-time Brownian motion, given by the first term in Eqs. [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF], is different from the variance of the TAMSD for a discretetime random walk with Gaussian jumps, first derived by Qian et al. [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF] and later analyzed in [START_REF] Grebenkov | Probability distribution of the timeaveraged mean-square displacement of a Gaussian process[END_REF][START_REF] Grebenkov | Optimal and suboptimal quadratic forms for noncentered Gaussian processes[END_REF]. In particular, when t ≪ T , the leading term of the EB parameter is 4t/(3T ) in the continuous case, and 2t/T in the discrete case. While the scaling form, t/T , is identical for both cases, the numerical prefactor is different. This difference is related to the fact that an erratic trajectory of Brownian motion between two successive discrete times is replaced by a jump in the discrete case.

There is a similar distinction between the continuoustime switching diffusion, studied in Sec. II, and discretetime Monte Carlo simulations presented here. In particular, the variance of the TAMSD, computed via the exact formulas [START_REF] Qian | Single particle tracking. Analysis of diffusion and flow in twodimensional systems[END_REF][START_REF] Andreanov | Time-averaged MSD of Brownian motion[END_REF], would differ from its numerical estimation by Monte Carlo simulations. In spite of this subtle difference between continuous-time and discretetime processes, the qualitative conclusion on distinct dif- fusivity states and insufficient measurement time as eventual causes of the apparent weak ergodicity breaking remains valid, as confirmed by numerical results of this section. The exact computation of the EB parameter for a discrete-time switching diffusion by adapting the general technique from [START_REF] Grebenkov | Probability distribution of the timeaveraged mean-square displacement of a Gaussian process[END_REF] presents an interesting perspective for future research.

B. Extension to multi-state models

The above computation can be formally extended to a switching model with J states, which is characterized by a set of diffusion coefficients D i and switching rates k ij . The propagators satisfy for each i 0 = 1, . . . , J:

∂ ∂t P i,i0 = D i ∂ 2 ∂x 2 P i,i0 + J j=1 k ji P j,i0 , (34)
where k ii ≡j =i k ij . The formula [START_REF] Weigel | Ergodic and nonergodic processes coexist in the plasma membrane as observed by single-molecule tracking[END_REF] for the marginal propagator for multiple successive points in the Fourier space remains practically unchanged:

P n (q n , t n -t n-1 ; . . . ; q 2 , t 2 -t 1 ; q 1 , t 1 ) (35) =    1 1 ... 1    † P(q n , t n -t n-1 ) . . . P(q 1 , t 1 )    p 1 p 2 ... p J    ,
where p i is the probability of starting in the state i, P(q, t) = exp(-M q t), and M q is now the J × J matrix of the form

M q = q 2 D -K † , with (D) ij = δ ij D i and (K) ij = k ij .
As the eigenvalues γ q of the matrix M q are determined as the zeros of the polynomial det(γI

-M q ) = 0, (36)
there is no explicit formula in general. However, as we are interested in evaluating the limit q → 0 to get the mean square displacement and other related quantities, one can apply the standard perturbation theory by treating q 2 D as a perturbation of the matrix M 0 = -K † . As the matrix M 0 is neither symmetric, nor invertible, the analysis is more involved. While analytical calculations for a general multi-state model seem challenging, numerical computations of the propagator and related quantities (such as the mean square displacement) are efficient due to the matrix form [START_REF] Lubelski | Nonergodicity Mimics Inhomogeneity in Single Particle Tracking[END_REF] when the number of states is not too high (say, below a thousand). It is also worth noting that multi-state switching diffusions can be seen as a discrete version of diffusing diffusivity models [START_REF] Chubynsky | Diffusing Diffusivity: A Model for Anomalous, yet Brownian, Diffusion[END_REF][START_REF] Jain | Diffusion in a Crowded,Rearranging Environment[END_REF][START_REF] Chechkin | Brownian yet Non-Gaussian Diffusion: From Superstatistics to Subordination of Diffusing Diffusivities[END_REF][START_REF] Lanoiselée | A model of non-Gaussian diffusion in heterogeneous media[END_REF][START_REF] Lanoiselée | Diffusion-limited reactions in dynamic heterogeneous media[END_REF][START_REF] Sposini | Random diffusivity from stochastic equations: comparison of two models for Brownian yet non-Gaussian diffusion[END_REF][START_REF] Sposini | First passage statistics for diffusing diffusivity[END_REF], in which the diffusivity D t of a tracer changes continuously in time (see [START_REF] Grebenkov | A unifying approach to first-passage time distributions in diffusing diffusivity and switching diffusion models[END_REF] for details).

C. Extensions to switching diffusion in bounded domains

While we focused on one-dimensional diffusion, the above computation can be carried on for more sophisticated processes.

Here, we briefly mention switching diffusion in an arbitrary bounded domain Ω ⊂ R d , for which the propagator for multiple successive points, P (x n , i n , t n ; x n-1 , i n-1 , t n-1 ; . . . ; x 0 , i 0 , 0), and its average over the state variables i k , P (x n , t n ; x n-1 , t n-1 ; . . . ; x 0 , 0), can be expressed in analogy to Eqs. [START_REF] Wong | Anomalous Diffusion Probes Microstructure Dynamics of Entangled F-Actin Networks[END_REF][START_REF] Jeon | In Vivo Anomalous Diffusion and Weak Ergodicity Breaking of Lipid Granules[END_REF]. For instance, Eq. ( 11) becomes

P (x n , t n ; x n-1 , t n-1 ; . . . ; x 0 , 0) = k1,...,kn u kn (x n ) u * kn (x n-1 )u kn-1 (x n-1 )u * kn-1 (x n-2 ) . . . u k1 (x 1 )u * k1 (x 0 ) ×P n ( λ kn , t n -t n-1 ; . . . ; λ k1 , t 1 ), (37)
where P n is given by Eq. ( 12) (or Eq. ( 35) in a multistate case), while λ n and u k (x) are the eigenvalues and L 2 -normalized eigenfunctions of the Laplace operator ∆, satisfying ∆u k (x) + λ k u k (x) = 0 in Ω. Boundary conditions determine the properties of the boundary in a standard way [START_REF] Redner | A Guide to First Passage Processes (Cam-bridge[END_REF]: Neumann condition incorporates passive impermeable walls whereas Dirichlet or Robin conditions allow one to describe diffusion-controlled reactions and the related first-passage time statistics, see [START_REF] Lanoiselée | Diffusion-limited reactions in dynamic heterogeneous media[END_REF] (note that P 1 (q, t) corresponds to Υ(t; q 2 ) in [START_REF] Lanoiselée | Diffusion-limited reactions in dynamic heterogeneous media[END_REF]). For instance, one can investigate diffusion-controlled reactions in presence of buffers that can reversely bind the diffusing molecule and thus affect its mobility. We emphasize that the above derivation requires that switching modifies only the diffusivity of the particle but does not change other properties (e.g., reactivity). In this case, restricted diffusion in each state is characterized by the same set of Laplacian eigenmodes. In contrast, this extension is not applicable to other intermittent processes [START_REF] Bénichou | Intermittent search strategies[END_REF] such as, e.g., surface-mediated diffusion [START_REF] Bénichou | Optimal Reaction Time for Surface-Mediated Diffusion[END_REF][START_REF] Bénichou | Mean First-Passage Time of Surface-Mediated Diffusion in Spherical Domains[END_REF][START_REF] Rojo | Enhanced diffusion through surface excursion: A master-equation approach to the narrow-escape-time problem[END_REF][START_REF] Rupprecht | Kinetics of Active Surface-Mediated Diffusion in Spherically Symmetric Domains[END_REF][START_REF] Rupprecht | Exact mean exit time for surface-mediated diffusion[END_REF], two-channel diffusion with switching reactivity [START_REF] Godec | First passage time statistics for two-channel diffusion[END_REF] or switching boundary conditions [START_REF] Bressloff | Escape from a potential well with a switching boundary[END_REF][START_REF] Lawley | Boundary Value Problems for Statistics of Diffusion in a Randomly Switching Environment: PDE and SDE Perspectives[END_REF].

From the propagator, one can compute the mean square displacement, as well as the mean and the variance of the TAMSD. However, these formulas include multiple infinite sums over eigenvalues and thus remain formal, so we do not present these results.

IV. CONCLUSION

We revisited the two-state switching diffusion model from a single-particle tracking perspective. The exact formulas for the mean and the variance of the TAMSD were derived, in order to investigate the behavior of the ergodicity breaking parameter. When the mean residence times 1/k 12 and 1/k 21 are small as compared to the measurement time T , switching is fast enough for a particle to probe both states so that a single trajectory is representative of the ensemble. In particular, the distribution of the TAMSD is narrow to allow for an accurate estimation of the mean diffusivity D. Roughly speaking, the twostate switching diffusion looks as normal diffusion with diffusivity D. In contrast, when 1/k 12 and/or 1/k 21 are comparable to or larger than T , an individual trajectory is not representative of the ensemble, the EB parameter may exceed 1, and the distribution of the TAMSD is broad. As a consequence, TAMSD curves exhibit a significant spread, in agreement with experimental observations in living cells. While the two-state switching diffusion is ergodic (when k 12 and k 21 are strictly positive), the analysis of individual trajectories may indicate non-ergodic features. As this fictitious non-ergodicity is the mere consequence of insufficient measurement time T , we called it "apparent weak ergodicity breaking". A similar behavior was observed for CTRW with an exponential cut-off: when the measurement time is much smaller than the cut-off time, the trajectory "looks" as that of usual non-ergodic CTRW but when the measurement time grows and exceeds the cut-off time, the nonergodic features disappear [START_REF] Lanoiselée | Revealing nonergodic dynamics in living cells from a single particle trajectory[END_REF]. The relative simplicity of switching diffusions presents their advantage as potential models for describing some intracellular and onmembrane transport [START_REF] Sungkaworn | Single-molecule imaging reveals receptor-G protein interactions at cell surface hot spots[END_REF][START_REF] Weron | Ergodicity breaking on the neuronal surface emerges from random switching between diffusive states[END_REF]. Moreover, these are natural models for describing the effect of buffers that can reversely bind the diffusing particle and thus affect its mobility. More generally, our results illustrate that a broad spread of TAMSD curves, often observed in experiments, is not necessarily a signature of weak ergodicity breaking. It may just be a consequence of the insufficient duration of acquired trajectories. Statistical tests of ergodicity breaking need to be systematically performed in such situations [START_REF] Lanoiselée | Revealing nonergodic dynamics in living cells from a single particle trajectory[END_REF][START_REF] Magdziarz | Anomalous diffusion: Testing ergodicity breaking in experimental data[END_REF].

1 FIG. 1 :

 11 FIG. 1: (Color online) The EB parameter χ as a function of the residence time 1/k for the case k12 = k21 = k/2, with t = 1, T = 1000, and several values of D2/D1 as indicated in the plot. Arbitrary units are used. Vertical dashed line indicates the measurement time T . Dotted horizontal line corresponding to normal diffusion (D1 = D2) is close to 4t/(3T ).

FIG. 2 :

 2 FIG. 2: (Color online) The EB parameter χ as a function of mean residence times 1/k12 and 1/k21, for t = 1, T = 1000, D1 = 1, and D2 = 0 (a) and D2 = 0.01 (b). Arbitrary units are used.

FIG. 3 :

 3 FIG. 3: (Color online) Spread of TAMSD curves obtained from a thousand of simulated trajectories of length T = 1000, with D1 = 1, pi = p eq i , and three sets of parameters: (a) D2 = 1 (normal diffusion); (b) D2 = 0.01 and k12 = k21 = 10 -1 ; and (c) D2 = 0.01, k12 = 10 -2 and k21 = 10 -3 . Arbitrary units are used.

FIG. 4 :

 4 FIG. 4: (Color online) Empirical distribution of TAMSD, δ 2 (t, T ), normalized by its mean value δ 2 (t, T ) , obtained from 10 000 simulated trajectories of length T = 1000, with D1 = 1, pi = p eq i , and three sets of parameters: (a) D2 = 1 (normal diffusion); (b) D2 = 0.01 and k12 = k21 = 10 -1 ; and (c) D2 = 0.01, k12 = 10 -2 and k21 = 10 -3 . Arbitrary units are used.

Acknowledgments

DG acknowledges the financial support by French National Research Agency (ANR Project ANR-13-JSV5-0006-01).