Geometry controlled anomalous diffusion in random fractal geometries : looking beyond the infinite cluster

We investigate the ergodic properties of a random walker performing (anomalous) diffusion on a random fractal geometry. Extensive Monte Carlo simulations of the motion of tracer particles on an ensemble of realisations of percolation clusters are performed for a wide range of percolation densities. Single trajectories of the tracer motion are analysed to quantify the time averaged mean squared displacement (MSD) and to compare this with the ensemble averaged MSD of the particle motion. Other complementary physical observables associated with ergodicity are studied, as well. It turns out that the time averaged MSD of individual realisations exhibits non-vanishing fluctuations even in the limit of very long observation times as the percolation density approaches the critical value. This apparent non-ergodic behaviour concurs with the ergodic behaviour on the ensemble averaged level. We demonstrate how the non-vanishing fluctuations in single particle trajectories are analytically expressed in terms of the fractal dimension and the cluster size distribution of the random geometry, thus being of purely geometrical origin. Moreover, we reveal that the convergence scaling law to ergodicity, which is known to be inversely proportional to the observation time T for ergodic diffusion processes, follows a power-law BT h with h o 1 due to the fractal structure of the accessible space. These results provide useful measures for differentiating the subdiffusion on random fractals from an otherwise closely related process, namely, fractional Brownian motion. Implications of our results on the analysis of single particle tracking experiments are provided.


Introduction
Roughly a century after Jean Perrin's groundbreaking experiments on Brownian motion using an elaborate single particle tracking method in 1908, 1 nowadays sophisticated single particle tracking of submicron tracer particles and even single molecules have become a routine tool in the study of passive and active transport dynamics in living biological cells as well as in various crowded fluids in vitro. 2 Distinguished from traditional ensemble averaging experiments single particle tracking enables one to directly access the diffusive properties of the tracer particles without any loss of information due to ensemble averaging.Based on the single particle tracking technique it was revealed that the tracer motion in 'superdense' 3 biological or complex environments often exhibits anomalous diffusion such that its ensemble averaged mean squared displacement (MSD) grows non-linearly with time in the form with the anomalous diffusion exponent a.We distinguish subdiffusion for 0 o a o 1 and superdiffusion with a 4 1.
][10][11][12] As examples from in vitro systems we mention the anomalous diffusion of microbeads in polymer networks or gels [13][14][15][16] and in artificially crowded media 17,18 as well as colloidal suspensions. 19,20Examples for superdiffusion due to active motion in living cells were provided in ref. 21-25.Effective superdiffusion on surfaces is also observed due to bulk mediation effects. 26,279][30][31][32] One process, that was identified as stochastic mechanism behind the motion of tracers in living biological cells and structured environments is the continuous time random walk. 33Continuous time random walks give rise to weakly non-ergodic subdiffusion due to multiple binding or caging events characterised by a long-tailed distribution of waiting times between successive motion events which features a diverging mean waiting time. 28,29,315][36][37][38][39][40] Such self-affine Gaussian processes have a fractal path and are ergodic, and they were shown to occur in viscoelastic environments. 40,4142,43 Fractals with their scale-free geometry were popularised by Mandelbrot in the 1980's to more closely resemble natural objects such as mountains or coastlines than do classical geometrical bodies. 441][52][53][54][55] Recently models of mixed origins of anomalous diffusion processes were also suggested. 29ollowing the advances in single particle tracking techniques it has become possible to garner sufficient evidence to diagnose the statistical characteristics of a stochastic process on the single trajectory level.7][58][59] Weakly non-ergodic processes in Bouchaud's sense exhibit the above-mentioned disparity between the ensemble averaged MSD (1) and the time averaged MSD defined below, even in the limit of long measurement times.Moreover, these processes exhibit ageing, the dependence of observables on the time span between initiation of the process and start of the measurement, and the amplitude of individual time averages fluctuates significantly.9][70] In contrast, stochastic motion driven by power-law correlated Gaussian noise and the random motion on fractals are the key ergodic anomalous diffusion processes. 29,713][74][75] It is known that the discrimination between them based on data analysis is challenging, as both models share the same asymptotic scaling behaviour of the velocity autocorrelation functions. 76Theoretical tests differentiating the two were recently proposed on the ground that the fractal dimension-or the number of visited sites-is smaller than the dimensionality of the embedding Euclidean space-or the total number of sites in space. 56,77We here explore in more detail the exact ergodic properties of de Gennes' random ant-in-thelabyrinth motion.Specifically, we analyse the single particle diffusion on two-dimensional percolation clusters at varying densities.Special emphasis is put on the ergodic properties of the motion revealed by the time averaged MSD.Remarkably, our simulations reveal that the fluctuations in the time averaged MSD have a unique statistical feature of geometrical origin which is significantly different from that of fractional Brownian motion.
This paper is organised as follows.In Section 2 we gather the information on the simulation procedure, the definitions of the averaging procedures used in our analysis, and the well known diffusion dynamics on percolation clusters in the framework of ensemble averages.In Section 3 we present our simulation results on the single trajectory level including the time averaged MSD, the amplitude scatter distribution of the time averaged MSD, and the ergodicity breaking parameter.In particular, we analyse the ro ˆle of clusters other than the infinite cluster.In Section 4 we discuss our results and conclude.

Simulation scheme
Our study is based on a two-dimensional square lattice with periodic boundary conditions, as commonly used in other studies. 50,77,78n this underlying structure we generate our random percolation environments.The system size varies depending on the load of the simulations.For instance, to simulate individual trajectories for analysis in terms of the time averaged MSD a lattice of 4096 Â 4096 is considered whereas in more computationally expensive simulations such as labelling clusters with different sizes, the lattice size is set to 1024 Â 1024.Each site is either filled with probability p or vacant with probability 1 À p.Here we follow the convention that p refers to the obstacle density in space. 50ercolation theory states that as the percolation density p approaches the critical value p c from above-that is for decreasing obstacle size-separated finite clusters merge and the correlation length x of clusters diverges as x B |p À p c | Àn with n = 5/36. 48t the percolation threshold p c and for lower obstacle densities p o p c there exists an infinite cluster of empty sites spanning the entire volume. 42,48The probability that a given site belongs to the infinite cluster at criticality scales as P N B ( p À p c ) b where b = 43/18. 42,48Concurrently as pp c the crossover length r c from anomalous to normal diffusion used below diverges as r c B | p À p c | Àn+b/2 .Thus at p = p c anomalous diffusion of the form (1) prevails at all times. 42,48,50Simultaneously the percolation geometry forms scale-free fractal objects over all length scales larger than the lattice constant. 42,48,50,79 For obstacle percolation considered here the critical percolation threshold is p c = 0.407256. 48 This critical case is depicted in Fig. 1.
This journal is © the Owner Societies 2015 In our simulations, after a random geometry with label z of given density p is generated a tracer particle starts to diffuse at its centre if it is vacant.Otherwise a randomly chosen vacant site close to the centre is chosen.In our Monte Carlo simulation, at every unit time step dt = 1 the tracer particle jumps to the nearest vacant site with probability 1/4 on our two-dimensional lattice.Uniform random numbers were generated by the ranðÞ function taken from Numerical Recipes. 81On an isolated unit cluster no move is allowed and the random walker becomes completely localised.In the limit of p -0 the diffusivity of the tracer particle is that of a free two-dimensional lattice, K ¼ a 2 4dt where a is the lattice constant.Our typical simulation time was 10 6 time steps, unless indicated otherwise.For the evaluation of physical observables N z random percolation geometries were generated at a given percolation density p, and N tracer trajectories for each percolation geometry were recorded.For the efficiency of the simulation study, the number N z was varied from 50 to 10 000 and N from 3 to 1000 depending on the observables.In each case sufficiently large samples were chosen to ensure a proper statistical convergence.Information on N z and N used in the evaluation of each observable are provided in the text and in the figure captions.

Ensemble averaged mean square displacement
Let us first address the potential ambiguity in defining ensemble averages in percolation systems. 42For a given fractal geometry z the MSD of an ensemble of particles starting at the origin at time t = 0 is given by where we choose the coordinate systems such that r i (0) = 0.The MSD hr(t) 2 i z will in general assume different values for different realisations z. ¶ To avoid the dependence on the specific realisation, an additional (disorder) average of the form is defined, where N z counts the different realisation z over which the average is taken.In what follows we refer to this quantity as the ensemble averaged MSD.We note that we did not average over initial positions for a given realisation z, however, we expect this averaging to be equivalent to the disorder averaging over different z.We also note that similar questions on the averaging arise in the analysis of diffusion in quenched energy landscapes. 82rior to our study of single trajectories below we here briefly summarise the ensemble averaged MSD for particles diffusing in percolation geometries, and we compare these results with our current Monte Carlo simulation.The ensemble averaged MSD has two distinct scaling regimes over time for percolation densities below the percolation threshold p c , namely 42,43,47,48,50 Thus, the diffusion is transiently anomalous below the crossover length r c ( p), corresponding to times shorter than the crossover time t c ' r dw c .Here the classical notation involves the walk dimension d w of the diffusing particle.Generally d w 4 2 implying subdiffusion with 0 o a = 2/d w o 1.In the simplest case when we have no obstacles ( p = 0) all sites are accessible to the random walker and we have normal diffusion at all times.In this case the crossover length r c vanishes.This case is shown in Fig. 2(a).At finite percolation densities below the percolation threshold, 0 o p o p c , the crossover length roughly corresponds to the average cluster size of locally fractal structures: local bottlenecks and cul-de-sacs generate anomalous diffusion, which eventually crosses over to normal diffusion when the typical distance travelled by the random walker exceeds r c .This case is shown in Fig. 2(b).We note that compared to the obstacle-free case p = 0, the effective diffusion coefficient in the case 0 o p o p c has a reduced value.At the percolation threshold p = p c E 0.407256 the percolation geometry is fractal on all available scales, and the crossover length r c diverges.Concurrently, the anomalous diffusion regime ranges over all time scales.8For the square lattice the anomalous diffusion exponent is 2/d w E 0.7, 50 which is again confirmed by our simulations, see Fig. 2. Indeed, Fig. 2 in computer generated percolation clusters milled into stacked plastic sheets. 83,84Above the percolation threshold the diffusion is more and more reduced.In particular, only finite, disjunct clusters remain such that the MSD eventually saturates to a plateau, see Fig. 2(d).

(c) documents the anomalous diffusion at criticality. Anomalous diffusion in effectively two dimensional geometries was indeed verified by NMR measurements
3 Probing the ergodic behaviour of diffusion on percolation geometries

Time averaged mean squared displacement
To explore the ergodic properties of the diffusion process on percolation clusters in addition to the ensemble averaged MSD we now consider the time averaged MSD defined in terms of the moving average of a single trajectory r(t). 28,29,31Here D is the lag time corresponding to the width of the window slid along the time series, and T is the total observation (measurement) time.The time averaged MSD (5) provides information on the diffusive properties of a single particle and is routinely applied to evaluate experimental single particle tracking data.Typically diffusive processes, that are ergodic in the Boltzmann-Khinchin sense, have the property that for sufficiently long observation times T (formally, T -N), the time and ensemble averaged MSDs converge to each other, lim In agreement with previous numerical studies, 77 our simulations indeed demonstrate that this ergodicity relation is satisfied for our random walker on the percolation geometry for any percolation density if we average the time averaged MSD over both N individual simulated trajectories (index i) and over an ensemble of N z percolation geometries z, As seen in Fig. 2, this mean time averaged MSD-represented by the yellow solid line-coincides with the ensemble averaged MSD (black circles) for all cases.The equivalence between the two quantities for N z = 50 is not significantly improved when the sample size is increased tenfold to N z = 500 (not shown).
However, the equality between the ensemble averaged MSD (7) and individual long time averaged MSDs is not always fulfilled.To demonstrate this fact we show in Fig. 3    .This is due to the fact that at criticality, when the percolation geometry is fractal, the size of the accessible regions is scale-free and finite accessible clusters of all sizes are created in addition to the incipient infinite cluster.When seeded on such a finite cluster the time averaged MSD of the emanating trajectory will show a plateau depending on the very size of the visited cluster.We thus encounter a non-ergodic scenario in which the irreproducibility of trajectories stems from the quenched geometry.This non-ergodic behaviour is strengthened above the percolation threshold, p 4 p c , as shown in Fig. 3(d).Individual traces d 2 ðDÞ exhibit a pronounced scatter around the average g d 2 ðDÞ D E .In the sense of Bouchaud's definition 32 the fact that once seeded on a specific cluster the random walker cannot pass to an unconnected cluster is referred to as strong ergodicity breaking.It denotes the case when the phase space of a system consists of disjoint subvolumes and access to each subvolume is determined by the initial position of the particle.It is intrinsic for the percolation system in that it is reproducible when simulations are performed with the same number of runs for given trajectory length T.Moreover, a typical scatter of amplitudes of the plateaus always occurs at any trajectory length T. It is in fact noteworthy that albeit these features of individual tracer motion the ensemble average remains ergodic within the error of our simulations.Subdiffusive continuous time random walks 31,60,65,67,85,86 show ageing in the sense that the particle encounters longer and longer waiting times on its motion due to the scale free nature of the distribution of waiting times.The amplitude of the time averaged MSD given by the effective diffusivity therefore becomes a decaying function of the observation time T. Analogous ageing effects occur for scaled Brownian motion, 87 heterogeneous diffusion processes, 68 and their combination.Apart from transient ageing 88 the Gaussian processes with correlated increments of the fractional Brownian motion and fractional Langevin equation motion types do not exhibit such effects. 9In Fig. 4 we analyse the dependence of the time averaged MSD on the observation time T using the data from Fig. 3 for fixed lag time D = 100.As can be seen no significant ageing can be detected, the traces settle towards a constant function of T. The figure also illustrates the relatively small amplitude scatter below and at the percolation threshold, p r p c , as compared to the behaviour at the higher percolation density, p 4 p c .

Time averaged occupation probability
To gather more quantitative clues on the above observations of the time averaged MSD we determine the profile of individual time averaged MSDs and explore how it is related to the particles' random paths on a given percolation geometry z.To this end we introduce the time averaged occupation probability P z ðrÞ at the lattice point r for a given fractal geometry z in the form where t r is the accumulated residence time of the particle at the lattice site r averaged over the observation time T. Fig. 5 presents typical patterns of the time average P z ðrÞ obtained from simulations with N = 10 4 .For the fully accessible two-dimensional lattice ( p = 0) the time averaged occupation probability is smoothly spread out from the origin with a uniform radial distribution of bell shape.This is a typically expected result for Brownian motion as well as for other ergodic processes in which the ensemble averaged occupation probability hP z (r)i is identical with its time averaged counterpart.This case is shown in Fig. 5(a).For growing percolation densities below the percolation threshold we therefore observe similar patterns, albeit local fine structure will appear due to the existence of a finite correlation range of obstacles measured by the radius r c introduced above.In particular, growing p will increasingly limit the accessible range for the random walker, that is, the occupation probability at finite T will decay faster from the initial position.At the percolation threshold p = p c the pattern of the time averaged occupation probability P z ðrÞ exhibits a structural fingerprint of the specific underlying fractal geometry.Interestingly, we find that P z ðrÞ may assume two distinct patterns at the percolation threshold.Thus, while Fig. 5(b) reveals an emerging fractal spreading pattern, Fig. 5(c) shows a uniform distribution within a localised region-note the different scales of the panels.How can this come about?This phenomenon is again related to the fractal nature of the percolation geometry at criticality.In the former case, the tracer particle diffuses on the infinite incipient fractal cluster of vacant sites whose linear size is only limited by the size of the simulated lattice.Thus the corresponding time averaged MSD exhibits the anomalous scaling law (1) and the ergodic property lim T!1 d 2 ðDÞ ¼ r 2 ðDÞ is satisfied. 77The quantity P z ðrÞ also provides information on how the tracer particle diffuses away from its starting site.It is noteworthy that the spatial gradient of P z ðrÞ has no obvious discontinuous shape.This means that there are no sites at which the particle is significantly trapped for much longer times than at other neighbouring sites.In contrast to this, the latter case corresponding to Fig. 5(c) represents the case in which over the observation time T the particle visits repeatedly a small range of sites.These are either extremely poorly connected to the infinite cluster or completely separated.Accordingly the corresponding time averaged MSD saturates (outset).As seen in our previous analysis above the percentage of such complete confinement is relatively low compared to the situation above the percolation threshold.In the latter case p 4 p c confinement will always occur, as no infinite cluster of vacant sites remains.Patterns such as those shown in Fig. 5(c) will then emerge for clusters of varying size.On average clusters of decreasing size emerge when p is increased (not shown).

Amplitude scatter distribution
We now quantify the amplitude fluctuations of individual time averaged MSDs seen in Fig. 3 by measuring the normalised distribution f(x) as functions of the lag time D and the observation time T in terms of the dimensionless variable 28,29,60,89 x ¼ Fig. 6 shows the variation of the scatter distribution f(x) at two different values of D and for fixed T = 10 6 .As anticipated from the time averaged MSD traces in Fig. 3 the distribution for p o p c represented by Fig. 3(a) and (b), has a narrow, bell-shaped form.
Its centre is located at the ergodic value x = 1, and the width becomes broader at larger lag times.This is a typical behaviour for ergodic diffusion. 89At and above the percolation threshold p Z p c this behaviour changes drastically, and the scatter distribution reveals unique features that are not expected from an ergodic diffusion process, nor are these features known from other weakly non-ergodic processes.Consider first the longer lag time D = 10 3 at p = p c .We first observe an additional peak showing up at around x = 0 (Fig. 6(c)).This means that there is a finite contribution from traces in which the particles undergo severely confined diffusion on finite clusters of vacant sites.This peak corresponds to the localised and uniformly distributed time averaged occupation probabilities P z ðrÞ shown in Fig. 5(c).Note that this peak at x = 0 is not a statistical error due to an insufficient number of simulation runs.Our further investigation demonstrates that this peak is still observed in larger simulation sets.Second, the position of the main peak is shifted to a larger value above x = 1 while the bellshaped profile is preserved.For the shorter lag time D = 10 the shift is less severe, however, one can distinguish a fine structure of different peaks for x values closer to zero.These correspond to clusters of different size, which can be resolved here when the overall values of the MSD are smaller.
Remarkably, when the percolation density exceeds the critical value p c the profile of f(x) changes again significantly.In this overcrowded situation vacant sites exist in the form of finite clusters, forcing the tracer particles to undergo confined diffusion.At the shorter lag time, over which the tracer's motion is not seriously hampered by confinement corresponding to Fig. 6(d) at D = 10, the distribution f has a dominant contribution at x E 1. Concurrently, several peaks close to x = 0 reveal again a distribution of cluster sizes.However, at longer lag times the main peak disappears, and the distribution is monotonically decreasing with x: localised motion becomes dominant for the statistic.
In Fig. 7 we investigate the variation of the scatter distribution at the percolation threshold, p = p c with changing observation times T for fixed lag time.It shows that the position of the main peak remains largely unchanged while the peak gets increasingly sharper as T is increased.In contrast the height and width of the peak at x E 0 is quite insensitive to T. Naively speaking this population splitting of different trajectories into a distribution around ergodic motion x = 1 and an almost immobile fraction again results from the coexistence of the two distinct diffusion modes induced by the geometry: unrestricted motion on large (infinite) clusters and confined diffusion on small, finite clusters.The propensities of occurrence depend on the lag time D as well as the percolation density p.At criticality both modes are significant.We note that this behaviour is a geometry controlled analogue of the dynamic population splitting into mobile and immobile particles in subdiffusive continuous time random walks 85,86 and heterogeneous diffusion processes. 90et us spin this idea forward with some analytical considerations based on the cluster size (area) distribution (s).Imagine the situation in which a random walker moves on a finite cluster of size s and radius of gyration R s ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 1 2s 2 ð Þ Â P s i;j r i À r j 2 s . 48In the limit of sufficiently long trajectories (T -N) and for lag times D Z D s , where D s is implicitly defined by the equivalence the time averaged MSD and the squared gyration radius, the time averaged MSD is saturated to the value d 2 ðDÞ ' R s 2 .As the cluster appears fractal on length scales smaller than R s , it can be reasonably assumed that the size of the cluster is s ' R d f s close to the percolation threshold, where d f is the fractal dimension.We then relate the distribution of the saturation value of the time averaged MSD to the cluster size distribution P(s) in terms of The remainder N À N u then corresponds to saturated diffusion.As demonstrated in Fig. 7 the scatter distribution for the free diffusion part will be a d peak as T is increased to infinity.Thus in the long time limit this part corresponds to the contribution . For the complementary, confined fraction the scatter distribution will be proportional to Therefore the normalised scatter distribution in the long time limit can be written as in terms of the variable d 2 .Here is a normalisation constant.From eqn (12) the averaged value of the time averaged MSD is calculated to be Although the exact form of P is unknown, the fact that P(s) is a monotonically decaying distribution leads to the relation found in Fig. 6 and 7. Eqn (12) tells us that the fluctuations of the time averaged MSD given by f is fully determined by the cluster size P distribution as well as the lag time D. In percolation theory the form based on the scaling function f (x) is invoked where t and s are two scaling exponents. 48,91At the percolation threshold p = p c the cluster size distribution asymptotically has the power-law form P(s) C s Àt with 2 o t o 3. 48 Using this information we find that Thus N u slowly decreases with increasing D as t = 187/91 E 2.05 at p = p c on the square lattice.The sharp peak in f(x) is observed even for long lag times, as evidenced in the plot for p = p c in Fig. 6.For the distribution of finite size clusters, with the above we find the following simplified expression of eqn (10), The theoretical value of the scaling exponent of d 2 in this relation can be evaluated from the numerical estimation of d f and t in our simulations.From approximately 25 000 clusters of various size s extracted from 9480 percolation geometries z, the fractal dimension d f is estimated via the scaling law s ' R d f s of the gyration radius as d f E 1.865 with the confidence interval (1.824, 1.905) and t E 2.03 from eqn (15).**This measurement leads to the prediction that the exponent is given by À1.028 with the confidence interval (À1.027,À1.029).In the inset of Fig. 7 the decaying part of the distribution f(x) at D = 10 3 is plotted in on a logarithmic scale (black circles).The slope of the linear fit (blue line) is À1.063 with the confidence interval (À1.226,À0.899), agreeing well with the theoretically predicted value.
Above the obstacle percolation threshold, p 4 p c , when only finite, disjunct cluster remain, the cluster size distribution P(s) has the asymptotic form where 1/c is the characteristic cluster size. 48From this result we obtain that in terms of the incomplete Gamma function G(Á,Á).The last transformation follows from the fact that for p 4 p c the fractal dimension d f = 2 equals the embedding Euclidean dimension, that is d = 2 here, corresponding to the locally fully connected, finite clusters.Consequently in this situation we face the scaling d 2 ' s.At large lag times D satisfying the criterion and thus we obtain N u { N. Therefore in eqn ( 12) the second term involving the cluster size distribution P will dominate the scatter distribution f.This argument supports the monotonically decaying profiles of f at D = 10 3 and 10 5 in the plot for p = 0.5 displayed in Fig. 7 which is accordingly the very profile of P.
This also underlines the fact that at time scales for which d 2 D ) 1=c almost all tracer particles undergo confined diffusion.In the opposite case, at short lag times satisfying d 2 D ( 1=c, we see , which is the same as N u at p = p c with the replacement d f -2.This means that there are particles performing free diffusion over short lag times D. Thus f should have a peak around the ergodic value, as shown in the case of D = 10 for p = 0.5 in Fig. 6.

Ergodicity breaking parameter
We now study the functional behaviour of the fluctuations of the time averaged MSD as the observation time T is increased.
For this purpose, we use the ergodicity breaking (EB) parameter 28,29,31,60 EBðDÞ ¼ In this case the fluctuation of the time averaged MSD tend to zero as T -N: for this ergodic process the time averaged MSD becomes identical to the ensemble averaged MSD.This conventional convergence continues if the space is filled with obstacles with a concentration p well below the percolation threshold p c , as evidenced for the case p = 0.3 in Fig. 8(d).
However, the convergence of the EB parameter behaves quite differently when we approach the percolation threshold at p = p c .In this case EB does not converge to zero when T goes to infinity but, as demonstrated in Fig. 8(d), EB converges to the finite residual value shown by the dotted horizontal line.This value was numerically estimated from the scatter distribution in Fig. 6.The geometrically induced fluctuations in the time averaged MSD due to the existence of finite clusters do not vanish even in the limit of infinite observation time.
The EB curves shown in Fig. 8 converge towards EB N in good agreement with the algebraic form Here h and EB N Z 0 are, respectively, the associated scaling exponent and the limiting value of EB at T -N, and k is a proportionality constant.In Fig. 10 we show the functional relations of the exponent h and of EB N versus the percolation density p as estimated from the EB curves.In Fig. 8 (left) the solid lines depict the best fit to eqn (23).We mention the following noteworthy aspects of the results: (i) for percolation densities p sufficiently below the percolation threshold p c the EB parameter follows that of normal Brownian diffusion, that is, h = 1 and EB N = 0, as mentioned above.(ii) For the opposite regime of high obstacle densities, p 4 p c , we find h = 1 and a non-zero value EB N of the residual EB parameter that tends to increase with the percolation density p.In this regime the Brownian convergence speed EB B T À1 is derived from the long time confined Brownian motion of tracer particles.The value EB N a 0 is attributed to the heterogeneity of the saturation values of the time averaged MSD and is thus geometry controlled.This statement is also consistent with the fact that EB N depends on the lag time D in Fig. 8(g) and (h), due to the different resolution set by the value d 2 D .(iii) Close to the percolation threshold p c the behaviour of the EB parameter is distinguished from these two regimes.As p is increased towards the critical point p c , h consistently decreases from unity.Thus when the space becomes fractal on all length scales we find that h E 0.80 and EB N a 0.
We emphasise that the estimated value h o 1 at p c is a genuine convergence property due to the fractal structure of the explored space.To rule out the possibility that the value h E 0.8 o 1 is due to the small finite clusters responsible for EB N a 0, we plot in Fig. 9 the EB curves after excluding the contribution of the smallest clusters of size s = 1.The unit sized clusters are the most dominant contribution among the finite clusters, see the form of P(s), and results in the vanishing time averaged MSD, d 2 ¼ 0. Fig. 9 shows that after removing these unit size clusters the EB parameter always converges to EB N = 0 for all p values.However, the convergence exponent h remains at h E 0.8, compare also Fig. 10.We note that the convergence law of the EB parameter is not the same as that of fractional Brownian motion, although both models share the same anomalous diffusion scaling (1) and are ergodic (for the present  [92][93][94] Thus, compared to fractional Brownian motion the fractal geometry-induced subdiffusion has a slower convergence to ergodicity.The observed characters of the EB parameter also differ from those of other diffusion models such as scaled Brownian motion, heterogeneous diffusion processes, and continuous time random walk. 60,68,95herefore, the EB convergence law provides useful information for unveiling the physical origins of anomalous diffusion processes found in complex random media.

Conclusion
Based on extensive Monte Carlo simulations we studied the ergodic properties of single particles diffusing in random twodimensional fractal geometries modelled by a percolation geometry at varying percolation density of obstacles.While the asymptotic equality between the ensemble averaged MSD .Using the time averaged occupation probability P z ðrÞ we demonstrated that these outliers correspond to trajectories when the particles motion is restricted on finite clusters of gyration radius R s significantly smaller than the system size.Other particles moving on the infinite cluster at criticality, however, do show the convergence to the ensemble averaged MSD g r 2 ðDÞ h ion the single trajectory level.We thus observe ergodic motion for a fraction of particles conditioned to move on the infinite  incipient cluster, as reported earlier. 77Due to the mix of confined and freely diffusing trajectories, the scatter distribution of the amplitudes of individual time averaged MSDs shown in Fig. 6  As shown in Fig. 7, upon increase of the observation time T the latter part is preserved while the main peak around the ergodic value The profile of the scatter distribution f(x) is quite different from those for critical and lower than critical percolation densities, see Fig. 6.At sufficiently long lag times D the distribution decays almost monotonically with the variable x, due to the exclusive presence of finite clusters with an exponentially decaying distribution of cluster sizes.This is a purely geometrical effect, as revealed in the distribution of the time averaged MSD which turns out to be independent of the observation time T. In both Fig. 8 and 9 the EB parameter was shown to decay as T À1 , corresponding to the convergence to zero of a Brownian particle.This is due to the fact that in this overcrowded obstacle regime the particle explores small clusters with local Euclidean geometries.We note that while the specific values for the scaling exponents and percolation thresholds vary for different lattice types and dimensionality of the embedding dimension the generic features revealed here remain unchanged.This claim is supported by preliminary studies on cubic and hexagonal lattices (not shown).
Anomalous diffusion in a fractal geometry is a non-Gaussian process 42,43,48 and therefore different from the Gaussian fractional Brownian motion on a fundamental level.However, except for the Gaussianity the difference between these two processes has not been studied in detail.In particular, it has been said that both models share the same ergodic behaviour.However, as revealed in our study the ergodic properties of diffusion on fractals displays distinctly different behaviour due to the quenched nature of the underlying geometry.Only in certain cases (low obstacle concentration or conditional seeding of the particle exclusively on the infinite cluster close to criticality) we observe ergodic behaviour.While the statistical fluctuations in the time averaged MSD are homogeneous for fractional Brownian motion, the fractal geometry-induced anomalous diffusion is heterogeneous, and the strength and character of the heterogeneity depend on the obstacle density.Therefore the detailed analysis of the ergodic properties is indeed a useful measure to differentiate the type of ant-inthe-labyrinth motion from other models, along with recently developed theoretical tools estimating the fractal dimension d f . 56,77An advantage of the method developed herein studying ergodic properties is more feasible than estimating the fractal dimension d f .
As experimental single particle tracking studies become increasingly popular our results are expected to be helpful in analysing and interpreting experimental results for various problems of anomalous diffusion in complex environments.In the case of lateral diffusion in phospholipid membranes some recent simulation studies reported that lipid diffusion is a two dimensional fractional Brownian motion 9,37 while conventionally it was understood as diffusion on fractal lattices. 43,50,96In addition for more complex membranes the lateral diffusion exhibits non-ergodic continuous time random walk type motion 97 or was identified as the combination of motion on a fractal and continuous time random walks. 8Also fairly complex non-Gaussian ergodic motion types were reported. 98Such a stochastic variety in the lateral diffusion dynamics seems to be natural given the fact that biological membranes have a composition dependent, wide range of structural complexities.From the trajectory analyses presented in this work one can have additional insight about the lateral dynamics and static structural complexity of a membrane system under investigation.For instance, if the plot of EB versus T gives the scaling exponent h o 1 it gives a signature that the lateral anomalous diffusion is not of fractional Brownian motion type.Then the ageing test for the time averaged MSDs along with the moment ratio evaluation further differentiates the fractal induced subdiffusion from non-ergodic continuous time random walk process.The distribution of saturated TA MSDs and a non-vanishing EB N may be used to obtain information on the size distribution of confining domains in a membrane, if any.Our analysis can also be applied to nano-particle transport in porous media. 52,99From the profile of the distribution of the time averaged MSD one may obtain information on the pore size distribution as well as the porosity of a medium.Additionally the ageing test for the time averaged MSD shown in Fig. 4 appears to be informative to examining the nano-particle-pore interactions: if the motion exhibits features of ageing it is likely that there are nonspecific interactions between the particles and a porous medium which give rise to the temporal heterogeneity in the time averaged MSD.

Fig. 1
Fig. 1 Sample realisation of a random percolation geometry at the percolation threshold p = p c on a square lattice of size 1024 Â 1024.The copper coloured structure representing the vacant sites of an infinite cluster accessible to the random walker constitutes a statistically scale-free, fractal object.
individual time averaged MSD traces d 2 ðDÞ along with the ensemble average g d 2 ðDÞ D E (black solid line) from Fig. 2. As long as p o p c the individual d 2 ðDÞ follow nicely the mean g d 2 ðDÞ D E except for long lag times D E T when the statistic for the time averaging becomes insufficient (Fig. 3(a) and (b)).In this case the cluster turns from a fully accessible two-dimensional lattice at p = 0 to a geometry with growing but localised clusters of obstacles at 0 o p o p c .The connectivity remains high such that, independent of the initial position, eventually the random walker explores the entire structure, apart from inaccessible areas of sufficiently small measure.Right at the percolation threshold p = p c shown in Fig. 3(c) we already observe significant deviations from the mean g d 2 ðDÞ D E : a smaller number of

Fig. 2
Fig. 2 Ensemble averaged MSD g r 2 ðtÞ h i (black circles) and mean time averaged MSD g d 2 ðDÞ D E (yellow solid line) for four different percolation densities: (a) p = 0, (b) p o p c , (c) p E p c , and (d) p 4 p c .In each panel the inset depicts g r 2 ðtÞ h i .t versus t on a double logarithmic scale, the slope is a À 1.Note the terminal Brownian scaling for p o p c with a reduced diffusivity for the finite percolation density (b).The numbers of the random geometries and trajectories used for the plot are N z = 50 and N = 1000.

Fig. 3
Fig. 3 Individual time averaged MSD curves d 2 ðDÞ and their mean g d 2 ðDÞ D E (thick black curve, as already shown in Fig. 2).In each panel the variation of the individual time averaged MSD curves represents 10 different trajectories on each of the overall 50 different simulated percolation geometries.

Fig. 5
Fig. 5 Time averaged occupation probability P z ðrÞ, eqn (8), distributed over three given percolation geometries.This quantity is sampled from an ensemble of 10 4 particle traces on the same percolation geometry z starting from the same initial site.We show results for a two-dimensional lattice corresponding to (a) p = 0.1, (b) at the percolation threshold with p = p c with an emerging fractal pattern, and (c) at the percolation threshold with a localised pattern.The outsets show examples of the corresponding time averaged MSD curves for those particles.

Fig. 6
Fig. 6 Normalised scatter distribution f(x) as function of the dimensionless variable x ¼ d 2 ðDÞ g d 2 ðDÞ D E for four different percolation densities.In each panel the curves represent f(x) for the lag times D = 10 (black) and 10 3 (yellow), each obtained from 30 000 simulation runs corresponding to N z = 10 000 random geometries and N = 3 different trajectories performed on each geometry.Each single run is of length T = 10 6 .
) where we used that d 2 ' s 2=d f omitting a proportionality constant.Given this latter scaling relation between the saturation value of the time averaged MSD and the size of a cluster, we can use both quantities interchangeably.Now assume that there are N long (T -N) trajectories of tracer particles that performed a random walk on a randomly generated fractal cluster and have a value d 2 governed by the distribution P. At a given finite lag time D, a tracer particle will perform anomalous diffusion of the form d 2 D 4K a D a as long as the walker has not yet fully sampled the cluster of size s or, equivalently, the time averaged MSD d 2 D has not reached the saturation value d 2 .Among N sample trajectories the fraction showing unrestricted anomalous diffusion will be

Fig. 7
Fig. 7 Dependence of the scatter distribution f(x) on the observation time T at the percolation threshold p = p c and for the lag time D = 10 2 .Each distribution was obtained from 30 000 sample trajectories.The inset shows the double-logarithmic plot of the part of the scatter distribution at small x for T = 10 6 .The straight line is a fit with the power-law function x Àb , see text.

Fig. 8 $ 2 3 D
shows the EB parameter as function of the observation time T for four distinct cases.In each panel two EB curves are plotted at lag times D = 10 and 10 2 .At given percolation density p the EB curves are shown on linear (left column) and double logarithmic (right column) scales.At p = 0 when the accessible ** The box counting method gives the estimation d f 1.9191 AE 0.0095 from space is the full two-dimensional square lattice the EB parameter displays the theoretically expected scaling behaviour EBðDÞ T for two-dimensional Brownian motion shown by the dashed slope in the double logarithmic panel for p = 0, see Fig. 8(b).

Fig. 8
Fig. 8 Ergodicity breaking parameter EB on linear scales (left column) and on double logarithmic scales (right column) as function of the observation time T. The red curves correspond to D = 10 2 and the blue curves represent D = 10 3 .The black solid lines depict the best fits with eqn (23).The grey lines in the double logarithmic plot for p = 0 (b), correspond to the EB parameter for normal Brownian motion with D = 10 2 and 10 3 .The dashed lines in the double logarithmic plots show the asymptotic behaviour of EB.The results are from N z = 3000 percolation geometries and N = 3 trajectories.

g r 2
ðDÞ h i and the time averaged MSD g d 2 ðDÞ D E averaged over all individual trajectories and many percolation geometries is observed at any percolation density, individual time averaged MSDs do not always behave like this average and are thus nonergodic.As we showed here this non-ergodic behaviour is geometry controlled and thus corresponds to strong ergodicity breaking due to a topologically disconnected phase space.Thus, at low obstacle densities p { p c the single particle diffusion exhibits typical ergodic behaviour as seen in the scaling law of the time averaged MSDs, their scatter distribution, and the ergodicity breaking parameter.In this case, the ensemble averaged MSD shows no disparity with individual time averaged MSDs as long as the observation time is sufficiently long.Close to the percolation threshold p E p c , however, such a typical ergodic character is no longer observed.There exists a fraction of time averaged MSDs which significantly deviate from the averaged curve g d 2 ðDÞ D E

Fig. 9
Fig. 9 Ergodicity breaking parameter as a function of the observation time T when clusters with s = 1 are removed.The double logarithmic plot shows a linear behaviour of EB versus the observation time T, independent of the percolation density p.The red curve corresponds to D = 10 2 and the blue curve represents D = 10 3 .The black solid lines for the linear scales plots are the best fits to eqn (23), the dashed black lines in the double logarithmic plots show the slope of the EB curves.Same numbers for N z and N as in Fig. 8.

Fig. 10 (
Fig. 10 (a) Variation of the scaling exponent h and (b) the residual ergodicity breaking parameter EB N as function of the percolation density p. Values from fit of the scaling function (23) to the EB curves in Fig. 8 and 9.The results are reported for D = 10.Red diamonds: fit from Fig. 8 including the smallest clusters of size s = 1.Blue circles: fit from Fig. 9 neglecting the smallest clusters.

d 2 ðDÞ g d 2 ðDÞ D E ¼ 1
becomes a sharp, almost d function like peak.Relating the cluster size distribution to the distribution of the time averaged MSDs by eqn(12) we qualitatively explained the observed behaviour of the scatter distribution f.An interesting behaviour is also observed for the ergodicity breaking parameter defined in eqn(21).From numerical analysis we found that EB generally follows an algebraic decay [EB(D) À EB N ] B T Àh( p) towards the finite residual EB parameter EB N reached at T -N, and h is a scaling exponent.The exact value of EB N depends on the percolation density p.As p is increased to p c the decay of the EB parameter thus deviates from that for typical ergodic diffusion, for which h = 1 and EB N = 0. Importantly, on approximating the percolation threshold, p E p c , the EB parameter has a slower convergence with scaling exponent h E 0.8 o 1 as well as a nonzero value EB N 4 0 due to the contribution of confined particles.Above p c the particle diffusion always takes place on confined clusters.Large fluctuations from the average g d 2 ðDÞ D E are present in individual time averaged MSDs.
and 7 acquires an asymmetric bell shaped form around the main peak at