Random walks in correlated diffusivity landscapes

In recent years, several experiments have highlighted a new type of diffusion anomaly, which was called Brownian yet non-Gaussian diffusion. In systems displaying this behavior, the mean squared displacement of the diffusing particles grows linearly in time, like in a normal diffusion, but the distribution of displacements is non-Gaussian. In situations when the convergence to Gaussian still takes place at longer times, the probability density of the displacements may show a persisting peak around the distribution’s mode, and the pathway of convergence to the Gaussian is unusual. One of the theoretical models showing such a behavior corresponds to a disordered system with local diffusion coefficients slowly varying in space. While the standard pathway to Gaussian, as proposed by the Central Limit Theorem, would assume that the peak, under the corresponding rescaling, smoothens and lowers in course of the time, in the model discussed, the peak, under rescaling, narrows and stays sharp. In the present work, we discuss the nature of this peak. On a coarse-grained level, the motion of the particles in the diffusivity landscape is described by continuous time random walks with correlations between waiting times and positions. The peak is due to strong spatiotemporal correlations along the trajectories of diffusing particles. Destroying these correlations while keeping the temporal structure of the process intact leads to the decay of the peak. We also note that the correlated CTRW model reproducing serial correlations between the waiting times along the trajectory fails to quantitatively reproduce the shape of the peak even for the decorrelated motion, while being quite accurate in the wings of the PDF. This shows the importance of high-order temporal correlations for the peak’s formation.


Introduction
The erratic motion of particles diffusing in a fluid medium (Brownian motion) has drawn considerable attention of scientists since Robert Brown first systematically investigated it [1].A. Einstein [2] was the first to propose a mathematical description of this type of motion (see [3] for a detailed historical account).Einstein, who essentially did not know about Brownian motion, found out, that such a phenomenon is an inavoidable consequence of the kinetic theory of heat, and closely connected it to diffusion.In this picture of what we now call normal diffusion, the particles' motion possesses two important properties [4]: (i) The mean square displacement (MSD) of the particles from their initial position grows linearly in time, (with d being the dimension of space, and D being the diffusion coefficient), and (ii) The probability density function (PDF) of the particles' displacements at a given time follows a Gaussian distribution The properties (1) and (2) were tested in many experiments, and their confirmation laid a solid foundation to our understanding of the atomistic structure of matter [5].The random walk approach used by Einstein assumed that one can approximate the particle's motion by a sequence of independent steps in random directions under the condition that the times necessary to make a step are the same, and the displacement in a single step has a finite second moment.This approach was closely mirrored in many early experiments using stroboscopic measurements.Independently of Einstein, Smoluchowski [6] presented a more formal mathematical description of the Brownian motion which lead to the same results as Einstein's, and set the ground to a new branch of probability theory concerning the diffusion processes [7].After Einstein and Smoluchowski, Langevin [8] proposed a new mathematical tool for the description of the particle's motion, the stochastic differential equation.The standard picture corresponds to the tracer's motion in a homogeneous, quiescent fluid.In the course of time, many deviations from this kind of behavior were found for other media.Numerous experiments on transport in complex media (disordered solids, rocks, biological media, etc.) showed that, instead of a linear time dependence as given by Eq. ( 1), the MSD often follows a powerlaw time-dependence ⟨r(t) 2 ⟩ ∝ t γ , with 0 < γ < 1 (subdiffusion) or 1 < γ < 2 (superdiffusion).A system whose MSD shows such a time dependence is said to exhibit anomalous diffusion.Depending on the specific case, different mathematical models have been proposed to describe this anomalous behavior by focusing on different aspects of the motion [9][10][11][12][13][14].Some classical models are: the uncorrelated continuous time random walks (CTRW) with power-law waiting time distributions, the fractional Brownian motion, and Lévy walks and Lévy flights.The PDF in these models may or may not be Gaussian.
Several recent experiments [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31] reported a new type of diffusion in which the MSD grows linearly in time, like in the normal diffusion, yet the PDF of displacements shows considerable deviations from the Gaussian shape.Usually, the PDF of displacements is well-described by a Laplace (two-sided exponential) distribution.This behavior was called Brownian yet non-Gaussian (BnG) diffusion [18].Some of the corresponding systems show a crossover from the non-Gaussian distribution to a Gaussian one at long times [15,30].In several cases [15,18,19,24,[26][27][28][29][30][31], for times at which the crossover takes place, the PDF presents a peak close to its mode.This peak resembles a part of the initial Laplace distribution, while the parts of the distribution further from its mode have already a more or less Gaussian shape.
Many of the systems in which the BnG diffusion is observed are pertinent to soft matter, and almost all of the experimental systems with BnG diffusion may show a great deal of spatial and temporal inhomogeneity, or disorder.Thus, the medium in which the particle moves may be spatially heterogeneous, or change in time.The properties of the diffusing tracer may change in time as well.
Different assumptions about the heterogeneity involved lead to different classes of models which were proposed for the description of BnG diffusion.The most popular class corresponds to the diffusing diffusivity (DD) models, see e.g.[32][33][34][35][36].They assume slow random changes of the diffusion coefficient in time.The particular variant of the model used in [34] will be called "the minimal model" of diffusing diffusivity in what follows.Another model describing BnG diffusion is the diffusivity landscape model (DLM) [37] which considers that the diffusion coefficient varies slowly in space.
The possible connection between the diffusing diffusivity and DLM was stated in Ref. [34]: the temporal randomness of the diffusion coefficient can be considered as stemming from its spatial change along the trajectory of a diffusing particle, so that the "minimal model" is a kind of a mean-field approximation for the case of spatial changes.Even if the DD and DL models are gauged in such a way that they reproduce the main features of the phenomenon, their predictions differ in some details.Looking particularly into these details may deliver valuable experimental insights into the kind of disorder involved.Thus, the DLM (and other models with correlated spatial disorder like the one discussed in [38], not necessarily exhibiting the BnG behavior) show a pronounced central peak at the mode of the PDF of particles' displacements.This central peak is, however, absent in the minimal model.The existence of this central peak in [38] was immediately connected to the correlated nature of disorder.
Recently, in Ref. [39], we concentrated on the behavior of the PDF of displacements close to its mode and showed that the PDF of displacements in several classical strongly disordered systems displays such a peak at its center.The behavior of this central peak is quite peculiar, since its presence shows that the convergence to a Gaussian (i.e.normal) behavior under homogenization may follow a different pathway than the one commonly known from the Central Limit Theorem (CLT) applied to sums of many independent, identically distributed (i.i.d.) random variables following some continuous distribution (in our case this should be the short-time Laplace one).This standard situation suggests that the initially sharp peak would smoothen and lower.However, under homogenization, the central peak in the considered classical strongly disordered systems gets narrower under the rescaling r → r/ √ t, p → t d/2 p implied by the CLT, while approximately keeping the height.Passing from the spatially disordered systems to their mean-filed counterparts (like the corresponding CTRWs, or the minimal model) restores the standard convergence pathway like the one predicted by the CLT.
The differences in the convergence pathways have to do with the fact that some important local information about the system is erased when passing to the pre-averaged (mean-field) description.Now, one could ask, what is the important information erased?In the present work, we try to answer this question by simulating the particles' trajectories in DLM (described as a continuous-time random walk of particles on a lattice with position-dependent waiting times) and erase the correlations between the waiting times and positions, while fully preserving the temporal structure of the walk.The result of the discussion shows that the existence of the persistent peak is connected to spatiotemporal correlations, and destroying them (while fully preserving the temporal structure of the problem) leads to a different kind of behavior.We note that the answer to this question may apply in other similar situations in strongly disordered systems.
The article is structured as follows: In Section 2, we revisit the diffusivity landscape model being the base of our investigation.Section 3 explores the idea that the DLM presents strong spario-temporal correlations which ultimately leads to the PDF exhibiting a central peak.We show that destroying these spatiotemporal correlations while fully preserving the temporal structure of steps reproduces the PDF in DLM at short times, but leads to lowering and disappearing of the peak at long ones.Section 4 provides a CTRW model with correlated waiting times which partially reproduces the behavior found in this decorrelated DLM model, but fails to fully describe the situation.In Section 5, we discuss the role of the particular shape of the correlation function of diffusivities assumed in DLM by considering a slightly different model.Finally, Section 6 presents concluding remarks.

Diffusivity landscape model
In what follows, we use the model proposed by Postnikov et al. [37] which assumes the particles' diffusion in a heterogeneous medium modeled by a correlated diffusivity landscape D(r).This motion is described by the force-free Langevin equation with multiplicative noise with ξ ξ ξ(t) being a Gaussian white noise with ⟨ξ ξ ξ(t)⟩ = 0 and ⟨ξ µ (t)ξ ν (t ′ )⟩ = δ µν δ(t − t ′ ) with µ, ν representing Cartesian coordinates.This Langevin equation corresponds to the Fokker-Planck equation with α being the interpretation parameter taking values in the interval 0 ≤ α ≤ 1 (see e.g.[40] for a comprehensive discussion).The authors of [37] asked, under which condition would Eq. ( 4) describe the BnG diffusion, and found out that the two following conditions should be met: First, Eq (3) without external potential must be interpreted in the Ito sense (α = 0; then, of course, any other interpretation can be used by introducing the corresponding deterministic force [40]), and second, initial positions of diffusing particles must be sampled from the equilibrium distribution.Taking as a "stylized fact" that the PDF at short times has been observed to follow a Laplace distribution [15][16][17][18][19][20][21][22][23][24][25][26][27][28][29][30][31], one can then show that the single-point PDF of the diffusion coefficients in the corresponding landscape should be given by a Gamma distribution: where Γ(•) is a Gamma function, and β and D are shape parameters dependent on the dimension of space.
In what follows, we will concentrate on the two-dimensional situation, for which β = 5/2 and D = 5D 0 /3, with D 0 being the sampled diffusion coefficient, i.e., the one defining the slope of the "experimental" MSD assumed to strictly follow the linear dependence ⟨r(t) 2 ⟩ = 2dD 0 t [37].
A finite-difference discretization of the Fokker-Planck equation, Eq. ( 4), with α = 0 on a square lattice with lattice constant a leads to a master equation (see Eq. ( 6) below) which, in its turn, defines a random walk scheme.The corresponding random walks are exactly what will be simulated in what follows.
For α = 0, Eq. ( 4) can be rewritten in the form and its discrete version is Here, the discretization point i corresponds to coordinates (x i , y i ) on a rectangular grid with the lattice constant a, and points j k are the four nearest neighbors of the lattice point i.
Under the above discretization, the random diffusivity field at each lattice point translates into correlated values of local parameters D i ≡ D(x i , y i ), which are generated according to the following algorithm, Ref. [37]: One begins by constructing an array of independent Gaussian random variables G i with zero mean and unit variance.Then one generates a correlated Gaussian field G i by applying the Fourier filtering method [41] to G i .Like in [37], we take the correlation function of the correlated field to follow with λ being the correlation length, and r ij the Euclidean distance between lattice points i and j.We note that the choice of Eq. ( 7) is not dictated by any physical reasons but by the ease of numerical implementation and further calculations.In Section 5, we will explore the consequences of changing the correlation function of the diffusivity landscape by considering a checkerboardlike diffusivity landscape.Finally, the correlated Gaussian field G i is transformed into the Γdistributed diffusivity landscape D i by performing a probability transformation: where erf(•) is the error function and F −1 β (x) is the inverse of the cumulative distribution function (CDF) F β (D) for the PDF given by Eq. ( 5), which is given by with γ(•, •) being the lower incomplete Gamma function.The procedure above generates a diffusivity landscape D i whose correlation function follows from that of the correlated Gaussian field, Eq. ( 7), by a transformation which will be discussed in Sec.4.1.
Figure 1 shows a realization of the diffusivity landscape D i for a lattice of 256 × 256 with a = 1, D 0 = 1 and λ = 10.Let us now return to our Eq.( 6).Defining the transition rates as with being the mean waiting time at a site, and 1/4 corresponding to the probability to choose one of the four neighbors to jump to.We put Eq. ( 6) into a standard form of a master equation which can be rewritten as In Ref. [37], Eq. ( 10) was solved using the forward Euler method.Here, we employ another approach.Like in Ref. [39], we use the fact that the master equation ( 10) corresponds to a CTRW with exponential waiting times distribution [42].Thus, to solve the master equation, i.e., to obtain the evolution of the PDF, we generate random walk trajectories whose waiting times follow the exponential waiting time density with τ i given by Eq. ( 9).Taking the lattice spacing to be the length unit of the problem (a = 1), we get Note that the single-step displacements in our CTRW are i.i.d.random variables (each step has a unit length and arbitrary, random direction), while waiting times are not independent since D i at neighboring points are correlated.An illustration of the procedure to generate random trajectories can be seen in panel (a) of Figure 2. As we shall see, this alternative method allows us to study the role of space-time correlations in the DLM, which would be impossible to do by solving the ordinary differential equations (ODEs).Moreover, generating random trajectories is considerably less computationally expensive than solving ODEs, and allows us to have much better statistics of the desired quantities.

Space-time correlations
In the last decades, several correlated CTRW models were proposed which lead to interesting behaviors [43][44][45].However, all these correlated models focus only on either the temporal part or the spatial part separately or simultaneously, but leave the spatiotemporal correlations out of the picture.The reason for this is that dealing with such cross-correlations is, in general, a very complex task, even from a computational point of view.All (semi-)analytical results usually come from applying mean-field techniques which partially or completely ignore the fine-scale structure of the system, so that some interesting features of the spatially disordered systems are not reproduced.This is the case, e.g., for the behavior of the central peak seen in the DLM [39], which is not reproduced in such pre-averaged models like the CTRW description of the DLM [39], or the minimal model of BnG [34].It is in this regard that we seek to know to what extent the spatiotemporal correlations are responsible for the persistence of the central peak and the unusual art of convergence to a Gaussian distribution by narrowing of the central peak under rescaling r → r/ √ t, p → t d/2 p, instead of its lowering.To assess the effects of spatiotemporal correlations in the DLM, we remove them by randomization of step directions and see what changes by comparing the PDF in the decorrelated motion with that in the correlated one.
Let us consider a particle whose initial position is r 0 , and follow its true motion as given by a random walk scheme corresponding to the master equation (6).At that position, the particle waits for a time t 0 which is drawn Springer Nature 2021 L A T E X template Random walks in correlated diffusivity landscapes 9 (a) (b) Fig. 2 Schematics of the procedure to decouple space and time.The real particle follows the blue trajectory, whereas the decoupled particle follows the red trajectory.Notice that the waiting times for both particles are the same, and they depend only on the positions of the real particle.
from the exponential distribution ψ(t|D 0 ) given in Eq. ( 11), with D 0 = D(r 0 ).Next, the particle randomly jumps to one of its neighboring lattice points whose position is r 1 , and then waits for another time t 1 which is now drawn from the exponential distributions ψ(t|D 1 ), with D 1 = D(r 1 ).This process of jumping and waiting is repeated until the maximal simulation time t max is exceeded by the sum of waiting times.At the end, the trajectory of our particle (which we will call real particle to distinguish its motion from its randomized counterparts) is given by a list of positions {r 0 , r 1 , r 2 , . . .}, which correspond to a simple random walk, and a list {t 0 , t 1 , t 2 , . . .} of the corresponding waiting times between subsequent jumps which are drawn from the exponential distributions ψ(t|D i ) with D i = D(r i ), and which are therefore dependent on the particles' positions.This dependence of the exponential distribution on the value of the local diffusion coefficient generates the spatiotemporal correlations in the DLM.The procedure to obtain the trajectories of a real particle is sketched in panel (a) of Figure 2. We now use the above trajectories of true motion of particles, which we refer to as "real trajectories" in what follows, to generate new trajectories in which space and time are uncorrelated.Let us start by taking the temporal part of a real trajectory, i.e., the list of waiting times {t 0 , t 1 , . . .}, and discard the spatial part.Then we proceed as follows: Let us consider a new particle, which we will call the decoupled particle, whose initial position is the same as for the real one, i.e., r 0 .The decoupled particle then waits a time equal to the first waiting time of the real particle, namely t 0 , and makes a jump to one of the neighboring lattice points with position r ′ 1 .At this position, the decoupled particle waits a time equal to the second waiting time of the real particle, namely t 1 , to make the next jump in a random direction.This process is repeated until the same number of steps as for a real particle is done, and the maximal time t max is exceeded.Hence, one ends up with a trajectory for the decoupled particle consisting of a list of the positions {r 0 , r ′ 1 , r ′ 2 , . . .} being sums of i.i.d.random steps, and the same list of waiting times {t 0 , t 1 , . . .} as for the real one, which are however decoupled from the corresponding particle's positions.This process is depicted in panel (b) of Figure 2. The trajectories of real and decoupled particles are then used for obtaining the PDFs of displacements in a given realization of the landscape.Similar PDFs are obtained for different realizations of the diffusivity landscapes and then weighted-averaged under the equilibrium condition: the corresponding weight is proportional to the waiting time t 0 at r 0 in the corresponding landscape.q(ξ) Fig. 3 A comparison of the one-dimensional cut of the PDF q(ξ) = p(x, 0)t of rescaled displacements ξ = x/ √ t for the real particle (black) and decoupled particle (red), see text for details.Each panel represents a particular time.The flattening of the the central peak in the PDF of positions of the decoupled particle is evident at longer times.
The resulting PDFs for real and decoupled particles are displayed in Figure 3.The four panels of the plot present four different maximal times: t max = 10 1 , 10 2 , 10 3 , and 10 4 .Each panel presents a comparison of the PDF for the real (black dots) and decoupled (red dots) particles.Each PDF is the average over 10 4 realizations of the diffusivity landscape over a lattice of 2048 × 2048 with λ = 10 and D 0 = 1.Each realization contains 10 5 particles.Plotted in Figure 3 is a cut of PDF p(x, y; t) through the origin at y = 0.Moreover, following [39], we plot the PDF as a function of the rescaled displacement ξ = x/ √ t.To keep the normalization of the PDF, it has to be rescaled as q(ξ) = t•p(ξ).Figure 3 shows that the decoupling of the spatial and temporal aspects of the motion changes the art of convergence to the Gaussian from the unusual one, by narrowing of the central peak, to the CLT-like convergence, by lowering and smoothening the peak.In Refs.[38] and [39], the existence of the central peak was connected with the set of particles which started their motion in a patch with a very low local diffusivity, so that they could hardly leave the patch until very long times.The randomization results show that this is only a partial explanation, since at the beginning of its motion a decoupled particle experinces the same, very long waiting times as the real one provided it started in such a patch.The trajectory of a decoupled particle is simply a different realization of a random walk with the same starting point associated with the same list of waiting times, so that only kind of correlations which are destroyed by our procedure correspond to what happens if the particle returns to a close vicinity of its initial position after making an excursion to the outside of the patch (a real particle will again experience long waiting times, while for a decoupled one these new waiting times are not necessarily long, and new long waiting periods occure at different positions).Thus, it is a behaviour after an excursion that makes the peak persistent.
Figure 3, however, shows that the central peak for a decoupled motion is still present at a times as long as 10 3 .The reason for its existence may only be connected with correlations between waiting times along the trajectory, which are not destroyed by decoupling.Therefore, our next step will be to include the temporal correlations into a space-time-decoupled CTRW model.

Time-correlated continuous-time random walk
Ref. [39] presented a mean-field description of the DLM.This mean-field description is constructed as an uncorrelated CTRW model whose waiting time distribution is found by averaging the mean waiting time distribution at a site (Eq.( 11)) over the distribution of the diffusion coefficients, which is given by the one-point distribution of the diffusivity landscape (Eq.( 5)).This mean-field waiting time distribution is given by for D 0 = 1.The corresponding waiting time density corresponds to a Pareto type II distribution with mean waiting time ⟨t⟩ = 1/4 and second moment ⟨t 2 ⟩ = 3/8.The fact that the initial state of the system in the DLM must be at equilibrium should also be included in its mean-field description.This is done by taking the first waiting time to follow the PDF [42] ψ which is also a Pareto type II distribution with a different exponent.
As we have shown in [39], this mean-field description, being a pre-averaged model (cf.Eq. ( 12)) neglecting all correlations, does not show any peak at the center of the distribution, except for a decaying remains of the initial condition at r(0) = 0.The PDF of the particles' displacements in this model is shown in Figure 5 to be compared with the results for a decoupled particle and for a CTRW model reproducing the serial correlations along the trajectory discussed below.

Correlation function of the diffusion coefficient along the trajectories
We would like to know to what extent the PDF for the decoupled particle can be replicated if temporal correlations are included.To do so, let us first determine the correlations of the diffusion coefficients along the trajectories of the random walk.Let us start by finding an approximation for the correlation function of the diffusivity landscape D(r) in terms of that of the correlated Gaussian field G defined in Eq. ( 7).The correlation function of the diffusivity landscape is, by definition, with σ 2 D the variance of the local diffusivity, and δD(r) = D(r) − D. In our case, D = 5/3 and σ 2 D = 10/9, for D 0 = 1 in two dimensions.Let us turn our attention to the mean ⟨D(0)D(r)⟩ in the last expression of Eq. ( 14).Just for convenience, let us denote D(0) and D(r) as D 1 and D 2 , respectively.By doing so, one can write Now, we make use of the invariance of the probability measures, being the bivariate distribution of the correlated Gaussian field used in the first stage of construction of the diffusivity landscape, with ρ(r) being the correlation function of this field given by Eq. (7).We note that the rdependence in this expression is fully due to the one of the correlation function ρ(r), and concentrate only on this ρ-dependence, introducing the function g( G 1 , G 2 ; ρ) = p( G 1 , G 2 ; r).Now, we can write Eq. ( 15) as with f ( G) being the function that transforms the correlated Gaussian field into the diffusivity landscape, Eq. ( 8).Note that according to Eqs. ( 16) and ( 8) the value of the function ⟨D 1 D 2 ⟩ is a function of ρ only, and therefore passing to the correlation function of diffusivity landscape, which differs from ⟨D(0)D(r)⟩ by shift and rescaling, we see that ξ(r) = ξ[ρ(r)], and the dependence ξ(ρ) is not influenced by a particular shape of the correlation function of the Gaussian landscape.Thus, the transformation from the Gaussian field to a Gamma-distributed landscape corresponds to a pointwise transformation of their correlation functions.This property will be used several times.The integration in Eq. ( 16) can only be performed numerically.However, one can still find an analytical approximation to this integral.We begin by Taylor expanding the function f ( G) around zero up to the fourth order: with a 0 = 1.4505, a 1 = 0.9704, a 2 = 0.2194, a 3 = 0.0130 and a 4 = 0.0011 for the values of parameters used.The coefficients correspond to the numerical evaluation of the analytical expressions of the corresponding derivatives of f which is easily done with Mathematica.This last expression can now be used to compute the integral in Eq. ( 16) as a function of ρ, since the corresponding integral reduces to the sum of moments of a bivariate Gaussian weighted with different prefactors.Keeping contributions up to the fourth order in ρ we find with c 1 = 0.918465 and c 2 = 0.081535.This simple quadratic approximation has the relative accuracy better than 0.0005 in the whole domain 0 ≤ ρ ≤ 1 as compared to the result of high-precision numerical integration.The transformation ξ(ρ) is invertible, and gives therefore the possibility to construct a Gaussian filed whose probability transformation would produce a Gamma-field with a desired two-point correlation function.We will use this possibility in what follows, when considering the correlated CTRW scheme in Sec.4.2.The inverse transformation is given by the solution of the quadratic equation, giving the inverse function Substituting the expression for ρ(r), Eq. ( 7), into Eq.( 17) we get the approximation for the correlation function of the diffusivity landscape: Now we can use this approximation to find the correlation function of the diffusivity landscape along the trajectories, or in other words, as a function of the number of steps C DD (n).To do so, Eq. ( 19) has to be averaged using the PDF f (r|n) of the displacements given the number of steps n: Fig. 4 Correlation function C DD (n) as a function of the number of steps for the twodimensional case with λ = 10 and D 0 = 1.We compare the numerical results obtained in simulations of particle diffusion in the diffusivity landscapes (green line), with the approximation Eq. ( 22) (red line).The standard errors of the mean (SEM) are represented by the light green area.Excellent agreement is observed in the whole range of the steps' numbers.
Given that the spatial part of the motion is a two-dimensional simple random walk, the PDF f (r|n) can be safely approximated by a two-dimensional Gaussian distribution with d = 2, a = 1 and, respectively, σ 2 = 1/2.Within this approximation, Eq. ( 20) takes the form Figure 4 shows a comparison between this approximate expression and the result from simulations of particle diffusion on the diffusivity landscape.Excellent agreement is observed in the whole range of steps' numbers.Note that the correlation function of diffusion coefficients is extremely long-ranged.

Correlated CTRW
Let us now use the correlation function of the diffusivity values along the trajectories, Eq. ( 22), to construct a time-correlated CTRW scheme.The process of generating correlated waiting times is similar to the one used for generating the landscape.Starting from values of ξ(n) = C DD (n) given by Eq. ( 22) we use Eq. ( 18) to obtain the correlation function ρ(n) of a Gaussian vector which then will be transformed to the one of diffusivity values and finally into waiting times along the trajectory.We proceed by generating a one-dimensional uncorrelated Gaussian vector g i by assigning to each entry of the vector a random number drawn from a Gaussian distribution with zero mean and unit variance.Then, using the Fourier filtering method [41], we generate a correlated Gaussian vector g i with correlation function ρ(n) with n = |i − j|.Using the probability transformation, Eq. ( 8), we transform this correlated Gaussian vector into a one-dimensional array of diffusion coefficients D i with the desired correlation function C DD (n) = ξ(n).The array of correlated diffusion coefficients D i is then used to generate waiting times of our CTRW scheme by drawing random numbers t i from an exponential waiting time distribution ψ(t|D i ) = 4D i exp(−4D i t).In each realization of the process, one repeats the procedure until obtaining such a number n ′ of drawn elements that the sum of the first n ′ elements does not exceed t max but the sum of the first n ′ + 1 does.The number of elements n ′ is then the number of steps performed by a walker until t max .The PDF of displacements for this correlated CTRW can be estimated by the average p(r, t max ) = ⟨f (r|n ′ )⟩ n ′ , with f (r|n) being the PDF of displacements for a given number of steps, Eq. ( 21), weighted with the waiting time of the first step.Fig. 5 A comparison of the one-dimensional cut of the PDF q(ξ) = p(x, 0)t of rescaled displacements ξ = x/ √ t for the decoupled particle (red) and the correlated CTRW (black) for t = 10 2 and 10 3 when the differences between the behaviors are considerable.The data for the decoupled particles are the same as in Figure 3, the result for coupled CTRW correspond to 8 • 10 6 independent realizations, see text for details.The green dots show the results for an uncorrelated CTRW model as given by Eqs. ( 12) and (13).
Figure 5 shows the resulting PDF for two different times t max = 10 2 , and 10 3 , one time per panel.Each panel presents a comparison between the PDF of the decoupled particle and that in the correlated CTRW.The PDFs for decoupled particles are the same as the ones in Figure 3 for the corresponding times.For the correlated CTRW, each PDF corresponds to the average over 8 × 10 6 different realizations of the correlated array of diffusion coefficients D i , constructed with λ = 10.As one can see, both PDFs are indistinguishable in their wings, and both present a central peak which, instead of narrowing, flattens out.However, the shapes of the peaks in both cases are significantly different.We note that the uncorrelated CTRW model shows a very different behavior in the wing (its convergence to a Gaussian is much faster) and does not show any peak except for some remains of initial condition at a shorter time.
It is worth mentioning that to generate the PDFs of the correlated CTRW we have used an extremely high number of realizations, namely 8 × 10 6 ; in our simulations this number was subdivided into five independent runs, and the results were both considered separately, and pooled for the plot in Figure 5.The analysis of the subsets shows that the height of the peak in different sets of 1.6 × 10 6 realizations still fluctuates considerably, so that this height is dominated by rare events, while both in the initial model and in the decoupled variant thereof the behavior in the peak may be considered as much more typical.
Since the approximations used to construct the correlated CTRW are quite accurate and the number of realizations is high enough to guarantee sufficiently good statistics, the differences suggest that our correlated model fails to capture important details of temporal correlations.Since serial correlations along the trajectories are reproduced correctly, one can conclude that these are the higher-oder correlations that play a key role in the development of the central peak but are of minor importance in the wings.Let us now go a few steps back and consider how critical our assumption about the shape of the correlation function ρ(r), Eq. ( 7), is, i.e., what happens if this function is chosen differently.To do so, we consider a DLM with a checkerboard-like diffusivity landscape.On a lattice, a checkerboard-like diffusivity landscape consists of an array of N × N squares containing 2ζ × 2ζ lattice points.A constant diffusion coefficient D (i) is assigned to each square.These diffusion coefficients are drawn from the distribution given by Eq. ( 5), the condition needed for the diffusion to be BnG.This choice of diffusivity landscape strongly changes the shape of the correlation function.Moreover, the changes in diffusion coefficients on the borders of the squares are now discontinuous, while the previous diffusivity landscape was assumed to model a smooth situation.In this model, ζ defines the correlation length of the landscape; to compare to the results of the above DLM, we set ζ = λ.As in the case of the DLM, we perform random walk simulations of particles diffusing on an ensemble of checkerboard-like landscapes, from which the PDF of displacements can then be constructed.Figure 7 shows the time evolution of the PDF averaged over 2 × 10 4 different realizations of the landscape, each one using 10 4 particles.The landscape was constructed with N = 109 and ζ = 10, i.e., we consider a lattice of size 2180 × 2180.One can see that the central peak is preserved.Moreover, the transition to the Gaussian limits follows the same type of convergence by its narrowing.This suggests that the form of the correlation function of the diffusivity landscape does not change the overall behavior.A closer look, though, reveals the presence of some discontinuities near the center of the distribution, which are expected from the fact that the diffusivity landscape itself is discontinuous.

Conclusions
In this work, we study the diffusivity landscape model (DLM) characterized by a diffusion coefficient slowly varying in space.Under specific conditions, this model leads to a Brownian yet non-Gaussian diffusion, that is, the MSD is linear in time, but the shape of the PDF changes from a Laplace distribution at short times to a Gaussian distribution at long ones.The art of convergence to the Gaussian is quite a peculiar one since the PDF at all times displays a central peak that does not decay with time, but narrows under rescaling.We show that the persistence of the peak is due to strong spatiotemporal correlations introduced by correlations of local diffusion coefficients in space.Destroying the spatiotemporal correlations on the level of single trajectories (by considering a different relaization of steps' directions while keeping the same list of waiting times as for the real motion) lets the peak to lower and to disappear at longer times.This kind of behavior is qualitatively reproduced by a correlated CTRW model with serial correlations of waiting times along the trajectory mimicking the ones observed in simulations.The model however fails to quantitatively reproduce the PDF for the decoupled case, showing a considerably lower peak.We attribute this fact to an important role of higherorder correlations which are not reproduced by the model.By considering a different variant of correlated disorder (the checkerboard model) we moreover show that the existence of the peak is insensitive to the exact shape of the correlation function of local diffusivities, and its shape is hardly sensitive to it.

Fig. 1 A
Fig. 1 A two-dimensional realization of the diffusivity landscape D(r) in the diffusivity landscape model.It corresponds to a 256 × 256 lattice with correlation length λ = 10 and sampled diffusion coefficient D 0 = 1.
), with b 0 = 2.77717, b 1 = 1.01867, b 2 = 0.09043, b 3 = 0.00101 and b 4 = 0.00003.The first coefficient (b 0 ) is equal to D 2 , therefore it vanishes when plugging back into Eq.(14).Moreover, since the coefficients b 3 and b 4 are small compared to b 1 and b 2 , they can be neglected.Under this approximation we get C DD (r) ≈ ξ[ρ(r)] with the function

Fig. 6 A
Fig. 6 A two dimensional realization of the diffusivity landscape D(r) for the checkerboard model.It corresponds to a 300 × 300 lattice where each cell has a size of 2ζ with ζ = 10, and sampled diffusion coefficient is D 0 = 1.

Figure 6
Random walks in correlated diffusivity landscapesshows one realization of the checkerboard-like diffusivity landscape on a lattice of 300 × 300 with ζ = 10 and D 0 = 1.

Fig. 7 A
Fig.7A one dimensional cut of the PDF q(ξ) = p(x, 0)t of rescaled displacements ξ = x/ √ t for the diffusion of particles in the checkerboard model.The straight line corresponds to the Laplace distribution, whereas the dotted line corresponds to the Gaussian distribution.The inset shows a close-up of the central part of the distribution exposing the peak.