Fractional Brownian motion with random diffusivity: emerging residual nonergodicity below the correlation time

Numerous examples for a priori unexpected non-Gaussian behaviour for normal and anomalous diffusion have recently been reported in single-particle tracking experiments. Here, we address the case of non-Gaussian anomalous diffusion in terms of a random-diffusivity mechanism in the presence of power-law correlated fractional Gaussian noise. We study the ergodic properties of this model via examining the ensemble- and time-averaged mean-squared displacements as well as the ergodicity breaking parameter EB quantifying the trajectory-to-trajectory fluctuations of the latter. For long measurement times, interesting crossover behaviour is found as function of the correlation time τ characterising the diffusivity dynamics. We unveil that at short lag times the EB parameter reaches a universal plateau. The corresponding residual value of EB is shown to depend only on τ and the trajectory length. The EB parameter at long lag times, however, follows the same power-law scaling as for fractional Brownian motion. We also determine a corresponding plateau at short lag times for the discrete representation of fractional Brownian motion, absent in the continuous-time formulation. These analytical predictions are in excellent agreement with results of computer simulations of the underlying stochastic processes. Our findings can help distinguishing and categorising certain nonergodic and non-Gaussian features of particle displacements, as observed in recent single-particle tracking experiments.

corresponding plateau at short lag times for the discrete representation of fractional Brownian motion, absent in the continuous-time formulation. These analytical predictions are in excellent agreement with results of computer simulations of the underlying stochastic processes. Our findings can help distinguishing and categorising certain nonergodic and non-Gaussian features of particle displacements, as observed in recent single-particle tracking experiments.
Keywords: stochastic processes, anomalous diffusion, fractional Brownian motion, diffusing diffusivity, weak ergodicity breaking (Some figures may appear in colour only in the online journal)
To allow for such a crossover, the concept of 'diffusing diffusivity' (DD) was introduced [67]. Introducing a fluctuating instantaneous D(t), in this model exponential PDFs at short times and Gaussian PDFs at long times emerge, while the process still features a linear MSD with time-independent effective diffusivity. The crossover is characterised by the correlation time inherent to the diffusivity dynamics. A similar concept of distributed diffusivities was previously developed in reference [68] (see also reference [69]. The DD approach helps mimicking and rationalising the impact of static and dynamical heterogeneities [70][71][72][73][74][75] on various statistical quantifiers of the particle-spreading dynamics.
Recently, various modifications and extensions of the DD model [67] were developed . Specifically, a minimal DD model [87] provided the general analytical subordinationbased framework for Brownian yet non-Gaussian processes. Persistent and antipersistent anomalous diffusion processes of FBM of generalised-master-equation family are normally also Gaussian. However, antipersistent non-Gaussian dynamics have been observed as well [99,100].
To accommodate this dynamics, recent advances of the DD model include the superstatistical FBM approach [99] describing the exponential PDFs as observed, e.g., for cytoplasmic RNA-protein diffusion in bacterial and eukaryotic cells [99]. A more general approach for the superstatistical generalised Langevin equation was developed in reference [101]. The link between the DD model and random-coefficient autoregressive model for non-Gaussian diffusion was established [102] and a DD model for generalised grey BM was developed [88].
Strongly non-Gaussian behaviours were reported for a number of complex systems, such as, e.g., molecular diffusion of lipid molecules or proteins embedded in protein-crowded lipid membranes [100,103,104], dynamics of polymers transiently adsorbed at solid-liquid interfaces [106,107], spreading dynamics of micron-size tracers in mucin-polymer gels [105,108], transiently superdiffusive spreading of amoeboid cells in heterogeneous populations [109], anomalous transport of tracers in amoeboid cells [110], dynamics of colloidal particles near a wall [111], in dense matrices of micropillars [71] and anisotropic liquid crystals [112], diffusion in narrow corrugated channels with fluctuating cross-sections [113], and the dynamics of acetylcholine receptors on live muscle-cell membranes [114].
The paper is organised as follows. In section 2 we introduce the physical observables used in the description and the basic equations solved in the text. In sections 3.1 and 3.2 the timeaveraged MSD (TAMSD) and the EB parameter of the DD-FBM model, respectively, are calculated analytically and supported by the results of computer simulations. We provide analytical expressions for EB for the BM case H = 1/2 and numerical results for the whole range H ∈ (0, 1). The discussion and conclusions are summarised in section 4.

Ensemble-versus time-averaging
Single-particle tracking (SPT) routinely measures the trajectories of submicron or even singlemolecular tracers in the form of time series of the particle position at unprecedented spatial and temporal resolution. SPT is by now an established powerful tool to study 'microscopic' diffusion in a broad spectrum of physical systems at different length-and time-scales. For a large number of SPT trajectories, the ensemble-based MSD is a well suited statistical measure. The dynamics is, however, often assessed from a limited number of long SPT trajectories in terms of the TAMSD. In accord with the ergodic hypothesis [19], the TAMSD for a given trajectory of a particle exploring the entire system for long times is equivalent to the MSD computed for a large ensemble of identical particles diffusing in the same system. In contrast, when the system features weak ergodicity breaking [115,116] the MSD and TAMSD cease to coincide, even for long measurement times T [16,19]. However, even for ergodic processes an ensemble of TAMSDs features a finite spread for finite trajectories. This spread is The position x(t) is driven by fractional Gaussian noise ζ H (t), whose amplitude is modulated by stochastic diffusivity D(t). The latter is taken to be the square of the Ornstein-Uhlenbeck process y(t).
The normal-diffusion DD model [87] and the DD-FBM model driven by power-law correlated noise ζ H (t) [38,39] are central for the current study, see figure 1. The ergodic properties of the DD-FBM process-defined below in the Boltzmann-Khinchin sense of the equivalence of the long-time limit of the MSD (2) and TAMSD (3) [16,19,121]-are investigated via analysing the realisation-to-realisation amplitude variation of individual TAMSDs quantified by EB, see equation (6) below. This is also a practical definition of ergodicity used, for instance, in statistical physics dealing with SPT data, as experimentalists often measure time averages. We do not talk about rigorous conditions of ergodicity, for instance, in the sense discussed in reference [122]. A detailed comparison of the dynamics and the (non-)ergodicity to the behaviours of pure FBM and the DD model is provided and the discrepancies are quantified below. In particular, we identify a fundamental time-scale below which EB features a plateau and we determine the residual EB value.

Definition of main observables
The MSD-the standard observable for a stochastic process [16][17][18][19]-is the average of the squared particle position with the PDF of its displacements at time t, Here and below we consider one-dimensional systems; component-wise extension to higher dimensions is possible. The TAMSD for a time series x i (t) of the ith particle is typically defined as where Δ is the lag time and T is the total measurement time. In contrast to the 'statistically averaged' MSD, each TAMSD is an inherently random quantity (even for BM) [16,19]. Therefore, averaging over N independent TAMSDs is often performed to compute the mean, For a fixed lag time Δ, a stochastic process is called ergodic if its MSD and TAMSD are identical in the limit of long observation times [16,19], i.e. when lim Δ/T→0 The randomness of each TAMSD realisation at a finite T gives rise to a certain amplitude scatter of δ 2 i (Δ, T) around the average (4). This scatter can be quantified by the EB parameter [16,19,39] where ξ(Δ) = δ 2 (Δ) δ 2 (Δ) . Similar to the TAMSD, the EB parameter is clearly also a function of the trajectory length, EB(Δ, T). Hereafter, however, we write EB as a function of its most relevant variable (often, the lag time Δ). For ergodic processes, EB approaches zero for long observation times T and the distribution of normalised TAMSDs [42,43,123,124] (named below φ(ξ)), approaches the Dirac δ-function in the asymptotic limit [16,19,44], φ(ξ) → δ(ξ − 1). In that limit, the results of all individual measurements coincide. For example, for paradigmatic BM (in the continuous-time limit) one gets [39,117,120] lim see also below. For nonergodic stochastic processes-such, e.g., as continuous-time random walks and heterogeneous diffusion processes-EB attains finite values at Δ/T → 0 [16,19,29,47,50]. We refer the reader to some examples of transiently nonergodic [43,124] and non-Gaussian [70,125] behaviour. Note also that the spectral content of single non-Brownian trajectories was recently investigated [22,126,127] and extended to random-diffusivity dynamics [128].

Main equations of the DD model
In what follows we employ the minimal DD model defined by the system of equations (following reference [87]) Here σ is a noise intensity, while ζ H (t) and η(t) is fractional Gaussian [38,39] and white Gaussian [16,19] noise, respectively. Both noises have zero means and correlation functions (for t 1 = t 2 ) and Here D 2H is the generalised diffusion coefficient. Extending the minimal DD model [87,88] the diffusivity D(t) in (8) is set to be the square of the Ornstein-Uhlenbeck process [129] with the correlation time τ . This guarantees non-negative diffusivities, D(t) = y 2 (t). The physical units of some model parameters and quantities are: We note here that the approach combining the features of both the FBM and DD models is pioneered in [130] and further developed here and, to the best of our knowledge, it has not been considered in the literature before.

Magnitude and distribution of the TAMSDs for the DD-FBM model
Here, we compute the mean TAMSD for the DD-FBM model for H ∈ (0, 1) and check the MSD-TAMSD equivalence. To be able to use the stationarity of the DD model, the initial condition for y(t) is chosen from the equilibrium distribution (see section 4.3 for nonequilibrium conditions). The mean TAMSD can be obtained via expanding (4) and calculating the position-correlation function (see appendix A for details), where the velocity autocorrelation function of the DD model is and For H = 1/2 the TAMSD (11) reduces to in agreement with the results of references [87,88]. The MSD for the DD-FBM model is obtained via integrating (8) with initial condition x(0) = 0, yielding From equations (11) and (15) follows that the MSD and mean TAMSD are identical in the entire range of (lag) times, As the DD-FBM process is self-averaging at T → ∞, the equivalence holds on the level of single TAMSD trajectories and the diffusion process is therefore ergodic. These theoretical predictions and results of computer simulations for the MSD and TAMSD in the DD-FBM model are presented in figure 2 for three values of the Hurst exponent H. Expression (15) for the MSD states that for persistent fluctuations (1 > H > 1/2) both for short and long times one obtains while for the antipersistent situation (0 < H < 1/2) a crossover from subdiffusive to Brownian behaviour is observed at long times, see figure 2. This behaviour observed in simulations is consistent with the analytical predictions stemming from the general MSD expression (15), as clarified in detail in reference [130]. Mathematically, repeating the arguments of reference [130], at short times (s 12 τ ) the diffusivity correlator in equation (15) approaches (1/2)σ 2 τ , see figure 7. This yields the anomalous MSD scaling at short times, both for the persistent and antipersistent situations. At long times we have to separate persistent and antipersistent motion. In the persistent case, 1 > H > 1/2, the leading contribution to the integral (15) at long times comes from large s 12 , owing to a slow decay of the noise-autocorrelation function. We thus have a monotonically decreasing correlator (18) with the limit K(s 12 τ ) → (1/π)σ 2 τ (see figure 7) that yields anomalous MSD growth [130] x In the antipersistent case, via splitting the integral (15) we arrive at the leading MSD contribution 4t ∞ 0 K(s 12 ) ζ H (s 1 )ζ H (s 2 ) ds 12 , where (due to convergence of the integrand) the upper limit was set to infinity, t → ∞ [130]. This leading MSD term at long times yields the linear growth, with the effective diffusivityD = lim 12 (where δ is the smoothening parameter of the correlation function [38,130]). Here, the crossover time from the short-time law (19) to the long-time linear diffusion scaling laws (20) and (21) is always the correlation time τ , independent on the actual value of the Hurst exponent H. Physically, the absence of the crossover and MSD scaling (19) in the entire rage of times for persistent motion and the crossover from the short-time behaviour (19) to the long-time linear MSD behaviour (21) for antipersistent noise is owing to the fact that antipersistent motion for FBM delicately depends on the exact vanishing of the cumulative correlation, in contrast, e.g., to an analogous process with cut-off noise correlator [46].
Specifically, from our simulations at Δ ≈ τ and H = 1/10 this crossover is distinct both for the MSD and mean TAMSD, see figure 2. We observe that the spread of individual TAMSDs for short lag times is larger for subdiffusion. From equation (A5) we also see that the mean TAMSD is independent of the total time T for all values of H. This is in contrast to ageing nonstationary processes [23,116,131,132]. Note that experimentally short-time plateaus for the MSD and mean TAMSD can emerge due to localisation errors of particle positions [133,134]. We also refer to the recent analysis of the effect of particle-localisation errors [135] on deviations of the short-time EB behaviour, above the BM asymptote (7).

Plateau values of EB
Here, we present the results for the EB parameter of the DD-FBM model in the limit of long traces, such that T τ . As for H = 1/2 the correlation function for fractional Gaussian noise, equation (9), reduces to the δ-function, we get exact analytical results for EB. For arbitrary H ∈ (0, 1) the complicated correlation function (9) hampers an exact analytical expression for EB in the DD-FBM model and, thus, we resort to simulations.

3.2.1.
Brownian case H = 1/2 for the DD-FBM model. The EB parameter of the DD-FBM model at H = 1/2 is calculated separately in the domains of lag times 0 < Δ < T/2 and T/2 < Δ < T, similarly as in references [117,120], and for long trajectories (T τ ). For 0 < Δ < T/2 we get (see appendix B) From equation (22), the leading order in 1/T yields and we get the respective asymptotes as We thus find that the standard result for BM [117,120] is reached at Δ τ , while a remarkable plateau is reached for EB DD+BM (Δ) at short lag times. As follows from (22) and (24), the correlation time τ emerges as the fundamental time-scale that controls the crossover behaviour of EB DD+BM (Δ). For For long enough T the first term in expression (25) gives At τ Δ and τ (T − Δ) this expression yields similar to EB of standard BM [117,120]. Towards the end of the trajectories, at Δ → T, the value is reached and from equation (27) one gets the first-order correction to this value as We note that equations (27) and (29) also hold for BM in the respective range of lag times [117].
In figure 3 the analytical and numerical results for the EB parameter of the DD-FBM model at H = 1/2 are presented. For the case τ T, EB DD+BM (Δ) starts from the plateau value 2τ/T (thin dashed lines in the plot, equation (24)) at short enough lag times, Δ τ . The BM EB asymptote (7) is approached for Δ τ , as equation (24) predicts. In figure 3 the results for varying correlation times are plotted. We find that, as for longer τ the EB plateau value increases (see equation (24)) and the region of lag times where EB DD+FBM (Δ) stays nearly constant becomes more extended. Concurrently, for larger τ values the region of lag times where EB DD+BM (Δ) follows the BM law (7) shifts towards larger Δ values, see the curve for τ = 1 in figure 3.

3.2.2.
General case of H ∈ (0, 1) for ordinary FBM: discreteness effects. We start by discussing the ergodic properties of free (unconstrained) ordinary or standard FBM. The expression for EB FBM (Δ) at short lag times Δ/T 1 was derived analytically in reference [39] in the continuous-time representation, namely where the coefficients are and The variation of C 1 (H ) is presented in figure 8. The main conclusion is that for short lag times EB cont FBM (Δ) scales linearly with Δ/T for H < 3/4, while the scaling of EB(Δ/T) is sublinear for 1 > H > 3/4 (H = 3/4 is a 'critical point' [39,44], see figure 5 below and the detailed behaviour in figure 9). In other words, the degree of reproducibility of individual TAMSD realisations increases linearly with the trace length for H < 3/4 and the statistical uncertainties decrease slower than linearly with T in the range of Hurst exponents 1 > H > 3/4 (it was erroneously predicted in reference [39] to diverge at H = 3/4, see below). The canonical continuous-time result for free BM, equation (7), follows from equation (30) at H = 1/2.
In the continuous-time formulation, the original EB results for FBM [39] were recently reexamined [44]. Specifically, within a more rigorous analytical framework for EB FBM (Δ) it was revealed that the 'critical point' at H = 3/4 disappears and the behaviour of EB is, in fact, continuous as function of H across this point, see figure 1 in reference [44]. This continuity agrees with the results of our FBM simulations presented in figure 5. Moreover, some unexpected results from computer simulations of FBM regarding the longer-tailed, non-Gaussian distributions φ(ξ) of individual TAMSDs-in particular, for progressively superdiffusive FBM and at longer lag times-were reported in figure 3 of reference [44].
The analytical predictions and the results of our computer simulations of FBM for EB FBM (Δ) are shown in figure 4. We find that for H > 3/4 the continuous-time theory [39] and our computer simulations coincide in a large range of the lag times studied. However, due to the innate limitations of the short-lag-time EB expansion (30), towards the end of the trajectory the simulations yield EB → 2 (28), while the (extrapolated) continuous-time prediction [39] would give EB(Δ → T) ≈ C 2 (H ), as follows from equation (32).
In contrast, for H < 3/4 at short lag times a plateau-like, saturation behaviour of EB is found, while the continuous-time theory [39] predicts a linear scaling with (Δ/T), see equation (30). This is the vital effect of a finite time-step used in computer simulations, δt. The region of EB saturation with Δ is particularly pronounced for small Hurst exponents, at which the EB plateau can occupy a considerable range of lag times (see, e.g., the curve for H = 1/100 in figure 4 and also the results of figure 10).
From the general discrete-time expression (C4) for EB disc FBM (Δ) we find that at Δ 1 = δt and H = 1/2 the residual value is approached, while for H → 0 one gets (see also reference [137]) Here N = T/δt is the number of elementary time-intervals in the trajectory. In the region we get the approximate, rather weak variation of the residual EB value with H (see appendix C), Note that (33) follows from this expression at H = 1/2. The lag time up to which this saturating EB behaviour is detected can be estimated as (see appendix C and figure 10 for details), where C 1 (H ) is given by equation (31). In the range from equation (C4) we get (in the leading order) that is identical to the result of continuous-time theory [39], see equation (30).
In figure 4 we plot the result for EB FBM (Δ) starting from the shortest lag time, Δ 1 = δt. For H > 3/4 the results for EB from computer simulations are in full agreement with the predictions of the continuous-time theory (30) and the discrete-time framework (39) for all lag times. At the end of the trajectories EB approaches the expected value EB = 2 (see equation (C7)). At short lag times EB FBM (Δ) in the discrete-time framework features a plateau for 0 < H < H pl ≈ 0.64 (as given by equation (C12)), as follows from the theoretical estimations (C9) and (C10). This plateau trend is most pronounced for the smallest Hurst exponents, extending towards longer lag times (as figure 10 explicitly quantifies). All these features are also consistent with the results of computer simulations, as illustrated in figure 4 for the  (7)) plotted versus the Hurst exponent H at the shortest lag time Δ 1 ≡ δt = 10 −3 . The results of our computer simulations are the filled symbols. The analytical results of the continuoustime theory [39], see equation (30), are the dashed coloured curves. The results of the discrete-time EB-derivation scheme-see equation (C4) and also reference [137]-is the solid black curve (shown for T = 10 2 only, not to cover the dashed coloured curves of reference [39] for H > 1/2). The thin dashed black lines are the discreteness-induced plateaus, equations (33) and (34). The thin vertical dotted line indicates the 'critical point' of the continuous-time EB theory for FBM [39]. The behaviour of EB near H = 3/4 is detailed in figure 9. Parameters: T = 10 0 , 10 1 , 10 2 for the respective colours, D 2H = 1/2. time-step δt = 10 −3 and in figure 11 for δt = 10 −6 . Note, however, that for a nearly ballistic Hurst exponent, at H = 0.99, in the region of extremely small lag times a rapid and unexpected reduction of EB for FBM is observed. This effect requires additional future consideration.
The results of computer simulations and the discrete-time-induced EB plateau (33) at Δ 1 are superimposing in the region 0 < H 1/2, see figure 5. The predictions of the continuoustime theory [39] deviate from the results of computer simulations in the range 0 < H 1/2, and, most pronouncedly, for very small values of the Hurst exponent.
In the region 0 < H < 3/4 the continuous-time EB(Δ) results for BM and FBM are linear in Δ/T, see equation (30). This linearity yields the universal, T-independent behaviour for the normalised quantity EB cont FBM (Δ)/EB cont BM (Δ) at short lag times in this region of H exponents, see the coloured dashed curves in figure 5. On the contrary, because of the sublinear scaling of EB cont FBM (Δ) with Δ/T in the region 3/4 < H < 1 the predictions for the ratio EB cont FBM (Δ)/EB cont BM (Δ) split up as T is being varied, see figure 5. We also find that the H-dependent residual values of the rescaled EB parameter, EB FBM (Δ 1 )/EB cont BM (Δ 1 ), at 0 < H 1/2, equation (36), are nearly independent on the trajectory length T, while in the range of Hurst exponents 1/2 H < 1 the ratio EB FBM (Δ 1 )/EB cont BM (Δ 1 ) is highly sensitive to T. The values of the ratio EB FBM (Δ 1 )/EB cont BM (Δ 1 ) decisively split up for different trajectory lengths in the range 1 > H 1/2, as demonstrated in figures 5 and 9. As all EB FBM (Δ 1 ) data in figure 5 are renormalised to the continuoustime classical result for the EB parameter of BM, EB cont BM (Δ 1 ), at H = 1/2 we find that the result of the discrete-time theory is a factor of 3/2 higher than that of the continuous-time EB calculation, while at H = 0 the EB value in the discrete model EB disc FBM (Δ 1 ) is a factor (3/2) 2 larger than EB cont BM (Δ 1 ) (see equations (7), (33) and (34)). Physically, the nonmonotonicity of the rescaled EB parameter EB FBM (Δ 1 )/EB cont BM (Δ 1 ) as a function of H presented in figure 4 is owing to the fact that for pure BM (at H = 1/2) the EB parameter attains the smallest value (natural for the most ergodic situation). The system becomes slightly less ergodic as H decreases from H = 1/2 towards H = 0 and the deviations from ergodicity turn much more dramatic as the Hurst exponent grows in the opposite direction from H = 1/2 towards H = 1 (the ultimate ballistic regime). This physically intuitive behaviour is consistent with all relevant limits clarified in figure 4, both within the continuous-time [39] and discrete-time (reference [137] and the current study) approaches. The overall dependence of EB FBM (Δ 1 )/EB cont BM (Δ 1 ) for FBM-for the results of continuoustime theory, the discrete-time analytical theory, and the innately discretised computer simulations-as a function of exponent H (see figure 4) as well as the universality of the ratio EB disc FBM (Δ 1 )/EB cont BM (Δ 1 ) for varying trajectory lengths in the region 0 < H 1/2 (see figure 5) are similar to those trends we observed for another Gaussian anomalous-diffusion process, namely so-called scaled BM (e.g., see figure 3(a) of our recent study [136]).
After the current study was finished, we became aware of reference [137] presenting detailed calculations of EB in the discrete-time scheme both for standard BM as well as for FBM. In figure 5 we show that the analytical results of equation (37) in reference [137]-which are identical to our EB derivation in equation (C4)-agree excellently with the results of our computer simulations. Note also that the emergence of the residual EB value for the discretetime simulations was already presented (but not rationalised) in figure 6 of the original study [39]. H ∈ (0, 1) for the DD-FBM model. We now consider the situation of FBM-driven DD motion. We present results from computer simulations for EB DD+FBM (Δ) at different H in figure 6 and compare them to the results for EB FBM (Δ) obtained in figure 4. For short lag times (Δ τ ) and for all values of H, EB is found to approach a plateau (subscript 'pl') with the residual value EB pl,DD+FBM ≈ 2τ/T. (40) In figure 12 we check this functional form via examining the plateau values for varying correlation time τ . The magnitude of EB at Δ τ approaches the plateau (33) (thin dashed line in figure 6) and it shows the FBM-like asymptotic law (30) at long lag times Δ τ (thick dashed lines in figure 6). As figures 6 and 12 show, EB for larger H approaches this plateau at progressively shorter lag times. This trend is similar to that of Δ disc pl for the discreteness-induced EB plateau for pure FBM, shown in figure 10. Performing computer simulations at varying correlation times τ for the relatively large Hurst exponent H = 9/10 we confirmed the universal EB plateau within the DD-FBM model, see figure 13, which is realised also (at even shorter lag times) for H = 9/10 in figure 6. The EB plateau for the DD-FBM model is reached at shorter lag times also for larger Hurst exponents. To make this region visible, in figure 14 we present results of simulations for shorter trajectories and shorter elementary lag-time-step used in simulations.

General case
We observe three clear differences between the FBM and DD-FBM models. (i) At short lag times, in the DD-FBM model the DD-induced EB plateau (40) is different from the discreteness-induced residual value (33) for FBM. Moreover, while equation (40) figure 15), EB DD+FBM (Δ) grows nonmonotonically in Δ and shows a minimum at Δ ≈ τ . We quantify the position of this minimum in figure 15, also verifying the minimum via examining the width of the distributions of individual TAMSDs in figure 16. Near this minimum, at intermediate lag times, EB DD+FBM (Δ) also features a drop below the paradigmatic BM asymptote (7). (iii) For long lag times, Δ τ , EB DD+FBM (Δ) at H < 1/2 approaches the continuous BM result, while for 1/2 < H < 1 the FBM limit for the EB parameter (equation (30)) is obtained, see figure 6.
To quantify inaccuracies of the numerical computation of the means δ 2 (Δ) and δ 2 (Δ) 2 , in figure 17 we present the respective error bars versus lag time for the computations within the DD-FBM model. We find that, as expected, the error bars grow in magnitude at later lag times, due to worsening statistics of averaging. In the plateau-like region of EB DD+FBM (Δ) versus H, however, the error bars are small enough for us to be confident in the existence of the plateau region itself and of the dip in EB DD+FBM (Δ) at Δ ∼ τ for very small H values (see figure 6). These features are not artefacts of poor averaging. Moreover, the magnitude of the error bars increases for larger exponents H, in agreement with equation (39).

Summary of the main results
We considered the combination of DD dynamics [67,87] and canonical FBM. Our main focus was to quantify the TAMSD fluctuations naturally occurring in experiments and gauged by the EB parameter. In particular, we analysed the plateau-like residual EB behaviour. Specifically, assuming the stationarity of the diffusivity distribution, the MSD and TAMSD of the DD-FBM model were studied. We found that the MSD and mean TAMSD are equal for both normal and anomalous diffusion in the entire range of (lag) times. For H < 1/2 we described a crossover in the TAMSD from subdiffusion to normal diffusion at lag times of the order of the DD correlation time τ , figure 2. From this TAMSD behaviour, the correlation time of heterogeneous environments featuring DD properties can potentially be extracted for SPT data sets. We revealed an intricate nonergodic behaviour in this DD-FBM model. Specifically, considering long trajectories with T τ , a crossover behaviour was found: for short lag times, Δ τ , the EB parameter was shown to approach a plateau with a residual value EB pl,DD+FBM ∼ 2(τ/T), which scales linearly with the ratio of the DD correlation time τ and the total measurement time T. Conversely, for long lag times Δ τ , EB DD+FBM (Δ) behaved the same way as for ordinary or standard FBM, see figure 6. The residual value of EB DD+FBM (Δ) was shown to be universal for all values of the Hurst exponents H.
Moreover, we demonstrated that for small values of H the variation of EB DD+FBM (Δ) was nonmonotonic, featuring a clear systematic minimum at Δ ≈ τ , see figure 6. Towards the end of the trajectories, at Δ → T, we found the expected value EB = 2. When simulating standard FBM, we found a plateau-like behaviour for EB FBM (Δ) at H 1/2 that scaled as the ratio of the time-step to the total trace length and depended weakly on the Hurst exponent, see figure 4. The plateau-like behaviour of EB pl,DD+FBM and EB disc pl,FBM (both analytically and via computer simulations) is the key result of the current study. The correlation time τ in the DD and DD-FBM models is, therefore, a fundamental time-scale for the ergodic behaviour, similar to the single time-step in the free discrete-time dynamics.

Other DD-related models
The relative standard deviation of fluctuations of individual TAMSDs ( √ EB in our notations) was rationalised for a model of the Langevin equation with time-dependent and fluctuating diffusivity in reference [80]. This model of multiplicatively coupled Langevin equations enables one to assess EB via studying the relaxation behaviour of the noise-coefficient matrix (or the matrix of instantaneous diffusion coefficients). In this approach the process B(t) = √ 2D(t) × 1 in the (multidimensional) Langevin equation was assumed to be ergodic. The general expression for EB was derived [80] in the continuous limit for arbitrary two-time-point correlation functions of the diffusivity matrix. For the onedimensional case, in the limit of long observation times and short lag times-and, additionally, when the relaxation time of the diffusivity (denoted below τ 1 ) is much longer than the lag time, the diffusivity relaxation function ψ 1 (t) decays fast enough, and the relaxation time is much shorter than the trajectory length (i.e., when τ 1 Δ and τ 1 T)-the approximate EB expression was derived as (see equation (33) in reference [80]) For the simplest (and most common) situation of exponential relaxation, ψ 1 (t) = ψ 1 (0)e −t/τ 1 ∝ e −t/τ 1 , from (42) it follows that when the lag time is the shortest time-scale in the problem EB saturates at This value is identical to our predictions for the EB plateau in the DD-FBM model, equation (40).
As derived in reference [80], for the (Markovian) two-state diffusion model-the Kärger model [138] (see also references [139,140] for its recent applications)-with the diffusion coefficients D 1 and D 2 > D 1 a dependence similar to equation (43) can be obtained for the EB parameter in reference [80]. Namely, for the transition rates from the state with diffusivity D 1 to the state with D 2 being k 12 (and vice versa for k 21 ), the equilibrium probabilities of the respective diffusion states are p 1 = k 21 /(k 12 + k 21 ) and p 2 = k 12 /(k 12 + k 21 ). The characteristic relaxation time is and in the same limit (at τ 1,2 Δ and τ 1,2 T) EB has the same functional dependence on τ 1,2 /T, that is [80] where ψ 1 (0) = p 1 p 2 (D 2 − D 1 ) 2 /(p 1 D 1 + p 2 D 2 ) 2 (see equation (57) in reference [80]). These results were also confirmed by computer simulations in reference [80]. The generalisation of these EB calculations for such a dichotomic stochastic-diffusivity model for the situation when switching between the diffusion states is governed by a power-law distribution was developed in the same group [81,82,84]. The short-lag-time plateau of EB appears to be universal for the models of diffusing or fluctuating diffusivity, also in the presence of anomalous underlying dynamics, as we demonstrated above for the DD-FBM model. The model of temporally fluctuating diffusivity has recently been applied by the same group to rationalise the dynamic interactions between membranebinding proteins and lipids in model biomembranes [104]. Finally, we refer the reader also to the discussion of short measurement times and 'apparent' ergodicity breaking for the two-state switching diffusion, recently presented in reference [98].

Effects of nonequilibrium initial conditions
Finally, the TAMSD and EB results obtained and discussed above involve the stationarity of the DD distribution. Nonequilibrium conditions can, however, also be relevant. For instance, recently the MSD of the DD model of normal diffusion with initial condition D(0) = 0 was discussed [88]. The MSD was demonstrated to be ballistic for t τ and linear for long times [88]. We found in the model (8) that the mean TAMSD for the initial condition D(0) = y(0) 2 = 0 and for H = 1/2 follows see figure 18. This expression converges to equation (14) as T → ∞ that is physically clear: the long measurement time eliminates effects of initial conditions in this system. The nonequilibrium initial DD conditions have no effect on EB, also in the general case H ∈ (0, 1).

Conclusions
Concluding, we here highlighted the role of the correlation time in the stochastic dynamics with random diffusivities. Similar to a finite elementary time-step in discrete-time diffusion models of free, unconstrained motion (BM, FBM), the correlation time represents the fundamental time-scale in the random-diffusivity dynamics: for lag times shorter than this correlation time the fluctuations of the TAMSDs reach a finite asymptotic spread, i.e., a finite residual value of the ergodicity breaking parameter EB. This effect is expected to be relevant for modern highresolution singe-particle-tracking experiments, to be considered in the data analysis. It will be interesting to analyse whether nonergodic anomalous-diffusion processes will exhibit similar features.  (3),

Acknowledgments
and using the MSD (15), we have Next, we compute the position autocorrelation function with and G(s 12 ) defined in equation (12)  (A5)

Appendix B. EB for the DD-FBM model at H=1/2
For H = 1/2 the 'second moment' of the TAMSD (after splitting the integrals into two parts), is the most challenging quantity to compute. Using the Isserlis-Wick theorem for Gaussian processes with zero mean, we expand all higher-order correlators via the pair correlators to get where To evaluate the integral (B1), we split the consideration into two cases, 0 < Δ < T/2 and T > Δ > T/2, that yields and Combining equations (B4), (B7) and (B8) with equation (6), we straightforwardly obtain the EB expressions of equations (23) and (25) in the main text.

Appendix C. EB for ordinary FBM
Here, we analyse the discrepancy between the theory and simulations of the EB parameter for FBM at short lag times, Δ T, arising from a discrete-time scheme employed in our simulations. Specifically, the TAMSD at the discrete points Using the Isserlis-Wick theorem (B2), from (C2) we get-as obtained initially in equation (A2) of reference [39]-that where t i − t j = (i − j) × δt. From equations (6) and (C2), the EB parameter for this discretetime FBM scheme can be expressed as At k n, using the Taylor expansion both for n = 1 or n 1, that coincides with equation (27), see also figure 4. At Δ 1 = δt and at H = 0 from (C4) we obtain expression (34) in the main text, while at H = 1 one gets EB FBM (Δ) = 2. (C7) For the Hurst exponents 0 < H 1/2 the first term in expression (C4) dominates and at Δ 1 = δt one gets the approximate expression (36). We checked these theoretical EB predictions at Δ 1 = δt versus FBM-based computer simulations in figure 5. From equation (C4) we also find that at Δ 1 = δt for H = 1/2 the EB parameter is For the case 0 < H < 3/4 in the limit n 1 via approximating the sum in equation (C4) by a continuous integral we obtain For H < 3/4 at n = 1 the sum is approximated by the leading term (at H < 1/2 only one term is enough). This ansatz works well for small Hurst exponents, while as H → 3/4 a progressively larger number of terms is to be accounted in the sum for an adequate approximation for the EB parameter. The condition of equality of the first and second term in equation (C9) provides a rough, H-independent estimate for the lag time, see equation (37) in the main text and equation (C10) below, below which the plateau-like, 'saturation' behaviour of EB disc FBM (Δ) is expected to occur. The analytical threshold for this lag time given by equation (37) is plotted in figure 10 versus H and shows that longer saturating regimes of EB disc FBM (Δ) emerge for smaller H values, consistent with the results of our simulations, as presented in figure 4.
In reality, however, the situation is more involved. If the plateau of EB persists in simulations for a long lag time, the first term in equation (C9) is typically much larger than the second one. This scenario is realised in the region 0 < H < 1/2 for Δ = Δ 1 . When the Hurst exponent increases and as H → 3/4 the two terms in equation (C9) become comparable in magnitude. Technically, based on equation (14) and figure 1 of reference [39] both illustrating the behaviour of the coefficient C 1 (H ) versus H, we find that equating the two terms in expression (C9) yields in this discrete-time scheme the following condition for the lag time This gives a universal, δt-independent condition for the critical lag-time value below which the EB plateau is realised. The value of the Hurst exponent via C 1 (H ) fully controls the EB plateau existence and the lag-time-range over which it persists. Therefore, if a plateau of EB disc FBM (Δ) is expected for lag times shorter than the elementary time step, at Δ < δt. This effect is thus 'undetectable' in our simulations (see figure 4). The condition (C11) is satisfied for large Hurst exponents 1 > H > 3/4 and also for in the region of small H (with H = 3/4 being the transition point between 'large' and 'small' H values). The shaded region in figure 10 demarcates the plateau-containing region of the EB behaviour for canonical FBM. Therefore, in the framework of this discrete-time FBM scheme we conclude using (C10) that at short lag times no plateau of EB disc FBM (Δ) is expected for the Hurst exponents in the range 1 > H H pl . This theoretical expectation is supported by the results of our computer simulations presented in figure 4.

Appendix D.
Below, we include additional figures supporting the claims presented in the main text.  (18), as predicted theoretically [132] and calculated from simulations, computed for the same parameters as in Figure 2.   (11) of reference [44] are the solid curves, the approximate analytical results of equation (12) of reference [44] are the asterisks and the results of the continuous-time infinite-trajectory-length theory of EB for FBM [39] are the dashed lines. All symbols and lines have the respective colours for varying lengths of the trajectory (see the legend). Other parameters are: Δ 1 = 10 −3 and D H = 1/2.         (6)) as obtained from our computer simulations. The bars are symmetric about the means (asymmetric in log-log scale). For error bars larger than the mean only the values above the mean are shown in logarithmic scale. Parameters are the same as in figure 6 and H=1/100, 1/2, and 9/10 (for the graphs from top to bottom, respectively).