Brownian yet non-Gaussian diffusion in heterogeneous media: from superstatistics to homogenization

We discuss the situations under which Brownian yet non-Gaussian (BnG) diffusion can be observed in the model of a particle’s motion in a random landscape of diffusion coefficients slowly varying in space (quenched disorder). Our conclusion is that such behavior is extremely unlikely in the situations when the particles, introduced into the system at random at t = 0, are observed from the preparation of the system on. However, it indeed may arise in the case when the diffusion (as described in Ito interpretation) is observed under equilibrated conditions. This paradigmatic situation can be translated into the model of the diffusion coefficient fluctuating in time along a trajectory, i.e. into a kind of the ‘diffusing diffusivity’ model.


Introduction
Many experiments and numerical simulations point onto a wide-spread (if not universal) behavior of displacements of tracers in different complex systems: the mean squared displacement (MSD) of the tracers grows linearly in time, (with D 0 being the diffusion coefficient and d being the dimension of space), like in the normal, Fickian diffusion; the probability density function (PDF) of the displacements is, however, strongly non-Gaussian. At first such non-Gaussian behavior accompanied by the MSD growth which is linear in time was observed in a broad class of materials close to glass and jamming transitions, such as binary Lennard-Jones mixture [1], dense colloidal hard sphere suspensions [2,3], silica melt [4], and bidisperse sheared granular materials [5].
Later it was found, that the PDF in most of these cases follows the exponential (Laplace) pattern with the parameter l(t) characterizing the width of the distribution [6]. In complex fluids this dependence was explained by the dynamic heterogeneity of ensemble of tracers which results in the intermittent nature of particles' trajectories [7], see reference [8] for the review. Later on, a very similar behavior was observed in a large amount of systems of a very different nature. Thus, in Wang et al [9,10], such a behavior was observed for the motion of colloidal beads on phospholipid tubes and in entangled actin suspensions, and later in the motion of tracer particles in mucin gels [11], the cases very different from systems discussed above. In [9] this type of behavior was termed Brownian yet non-Gaussian (BnG) diffusion. More experimental examples are mentioned in [12], for recent simulation results see [13].
In [10] the behavior was attributed to the heterogeneity of the medium. Each tracer moves in the environment with its own diffusivity D, i.e., the PDF of each tracer's displacement follows (in one dimension, or in the projection on the x-direction). The PDF of displacements in an ensemble of tracers is than given by where p(D) is the probabitity density of the distribution of the diffusivities D. If p(D) is exponential, than P(x, t) is a Laplace distribution. Different diffusivities for different members of the ensemble may be connected with the fact that the properties of the medium slowly change in space or in time.
The heterogeneity can be attributed to the tracers themselves. Heterogeneous ensembles of tracers were invoked in several works. Thus, analyzing the consequences of the fundamental observation in population biology that individuals of the same species are not identical, it was shown in [14] that e.g., the Maxwell type of the distribution of speeds of flying insects may lead to the exponential distribution of diffusivities. Analyzing the movement data of the parasitic nematodes, Hapca et al [15] showed that there is a significant variation of an effective diffusion coefficient within the population, and that the distribution of the diffusion coefficients follows a gamma distribution. This again leads to the exponential tails of the displacement distribution.
In mathematical statistics the procedure given by equation (4) is called compounding [16]. A slightly different procedure formulated in the language of random variables, is associated with the concept of generalized gray Brownian motion (ggBm) [17][18][19][20]. For Brownian diffusion characterized by the Gaussian displacement distribution, equation (3), this procedure leads to a similar scheme as simple compounding. In physics, such procedures fall under the notion of superstatistics [21][22][23].
All above means that the BnG diffusion appears in many different physical systems and mathematical models, and, in physics, this is probably not a monocausal phenomenon.
Similar effects are also seen in anomalous diffusion [24][25][26][27][28] which case however is not discussed in the present work.
In many cases the form of the PDF changes at longer times for the normal, Gaussian one, see e.g. [9,10]. Such change cannot be described within the tracers' heterogeneity model unless the properties of tracers change in time (which is not the case in the experiments using well-controlled tracers). In the models of heterogeneous media, the exponential-to-Gaussian transition in the PDF may stem from slow temporal fluctuations of the diffusion coefficient (a diffusing diffusivity model [12,[29][30][31][32]).
As mentioned in [12,15], such a transition may also stem from the quenched spacial heterogeneity of the medium (quenched disorder). At short times the motion of a tracer is confined to a domain with some diffusivity D, while at longer times it travels between the domains with different diffusivities, so that the values of the diffusion coefficients fluctuate along the trajectory. The main difference between the quenched disorder model and different variants of diffusing diffusivity models is therefore the fact that now, having returned to an already visited patch of the system, the walker perceives the same diffusion coefficient, and not a different one.
In the present work we consider the motion of tracers in a medium with quenched distribution of diffusion coefficients slowly varying in space. Thus, at short times different tracers see different diffusion coefficients while at long times (and correspondingly large scales) the homogenization sets in, and the motion of each tracer can be described as taking place in a homogeneous medium characterized by some effective diffusion coefficient D * . The BnG diffusion implies that the mean diffusion coefficient D 0 sampled at short times is equal to the effective one D * in the homogenized regime. The main topic of the present discussion is to find out when this is likely to be the case, and, if it is case, how does the transition between the short-time Laplace and the long-time Gaussian PDF shapes takes place. As we proceed to show, the possibility to observe the BnG diffusion strongly depends on the initial conditions. Thus, for the case when the particles tracked are introduced into the system at random at the beginning of the observation, the BnG diffusion is not observed. On the other hand, if the particles were present in the system, and their distribution could equilibrate before the observation started, the BnG diffusion is indeed observed for some specific model of the medium.

Position-dependent diffusion coefficient
Let us assume that the BnG diffusion phenomenon in some specific medium is fully due to the spacial heterogeneity of the local diffusion coefficient D(x) which is position-dependent. In dimensions higher than d = 1 the system will be assumed to be isotropic on the average, which considerably simplifies the further discussion. In order to observe heterogeneity at times which are not very short (much longer than experimental sampling time for the trajectory), the local diffusion coefficient D(x) has to vary only slowly in space. The correlation length λ of D(x) is therefore mesoscopic, and the characteristic time of the transition from the short time inhomogeneous (superstatistical) to the homogenized behavior is t H ∼ λ 2 /D 0 .
The knowledge of the diffusion coefficient D(x) as a function of the coordinates is not enough to uniquely define the properties of this diffusion. Systems with exactly the same D(x), the ones described by the Langevin equation equation (5), are described by different Fokker-Planck equations (FPEs), depending on what interpretation of the stochastic integral is assumed. The corresponding FPEs for a given realization with the interpretation parameter α taking the values in the interval 0 α 1. The typical interpretations, the Ito, Stratonovich and the Hänggi-Klimontovich (HK) ones, correspond to α = 0, 1/2 and 1, respectively. Each of these schemes appears in different physical situations (see e.g. reference [33]), i.e. the choice of the integration prescription in stochastic calculus implies different assumptions about our disordered medium. These differences get clear from the difference in numerical implementations of the master equations corresponding to the discretization of equation (6) with the constant C > 0 depending on the size of the system and on the number of particles therein. The HK case is the only one when this profile is flat. The Fokker-Planck equation for HK interpretation corresponds to the phenomenological second Fick's law The Ito interpretation leads to the Fokker-Planck equation which, in some cases, better reproduces the experimental results for heterogeneous diffusion and is connected with the continuous time random walk (CTRW) scheme [34].
In what follows we will provide rigorous inequalities for the long-time diffusion coefficient D * for different values of interpretation parameter α under different additional conditions, as well as estimates for this long-time diffusion coefficient as provided by the effective medium approximation (EMA). This will allow us to discuss the question, which model is the best candidate for reproducing the BnG behavior. We will moreover present simulation results for the time-dependence of the diffusion coefficient and of the corresponding PDF for systems corresponding to HK and Ito models, the ones with α = 1 and α = 0, respectively, which correspond to the barrier and to the trap model of the disordered medium. These simulations will confirm our theoretical assessment that, within the class of models with quenched disorder, the Ito model under equilibrium conditions is the only promising candidate for the description of the BnG diffusion.

Sampled diffusion coefficients and local diffusivities
In heterogeneous media, the situations under equilibrium, and the non-equilibrium ones may show very different properties [35,36]. The first, equilibrium, situation corresponds to the case when the system (medium, already containing tracers) was prepared long before starting the observation. The position of the tracer is monitored from the beginning of observation (t = 0) on. If there are several tracers in the system, the ones to observe are chosen at random. In this case the probability density to find a tracer at position x at t = 0 is proportional to the equilibrium density n(x) of tracers at the corresponding position. If n(x) is not constant, the space is not sampled homogeneously. The second, opposite, situation corresponds to the case when the tracers were introduced at random exactly at the beginning of the observation, and the diffusivity values at initial particles' positions are sampled according to the distribution of D(x). We will call the first and the second situations the equilibrium sampling, and the homogeneous sampling, respectively. In [11], a moving time averaging over the long data acquisition time is used to get both the displacement's PDFs and the MSDs, which assumes that the system has enough time to equilibrate. In [9] the ensemble average was used, but time between preparation and the beginning of the observation was not specified. The situation will be equilibrated if this time is much larger than t H .

Distribution of sampled diffusion coefficients
At short times, i.e. in the superstatistical regime, the PDF of particles' displacements p(x, t) follows by averaging the Gaussian PDF of particles' displacements in a patch with local value of the diffusion coefficient D(x) over the distribution of these local diffusivities close to the particles' initial positions. Therefore, the PDF in the superstatistical regime (when each particle can be considered as moving with its own diffusion coefficient) is given by where p S (D) is the PDF of sampled diffusion coefficients. This PDF is defined by equation (10) and may or may not be equal to the one-point PDF p(D) of the random field D(x), depending on the type of sampling (equilibrium or homogeneous) implied in experiment. Let us take as a 'stylized fact' that the PDF of displacements at short times has the exponential form with l(t) defining the width of the distribution and A(t) being the normalization constant. Requiring that p(x, t)dx = 1 and that in d dimensions, we obtain the explicit form of the displacements' PDFs Passing to the Fourier transforms in the spacial coordinate in equation (10) and taking into account the central (rotational) symmetry of the PDF we get where k = |k|, and therefore see thatp(k, t) is the Laplace transform of p S (D) taken at the value of the Laplace variable s = tk 2 , and that p S (D) thus follows as the corresponding inverse Laplace transform. The values ofp(k, t) arep and therefore p S (D), following as the inverse Laplace transform of these equations in s = tk 2 , are: These distributions possess the form of a Gamma-distribution with the shape parameter β = (d + 1)/2 (with d being the dimension of space) and with the mean , which will be of use in what follows, diverges in d = 1 and is finite in higher dimensions:

Distribution of local diffusion coefficients
In the superstatistical regime corresponding to short times the particles move in different areas with different local diffusion coefficients D(x). Since the patches do not have to be sampled homogeneously, the distribution of sampled diffusion coefficients might differ from such of D(x). In situations corresponding to homogeneous sampling, the two PDFs, p(D) and p S (D) coincide, so that p(D) is given by equation (12). The same is true also for HK interpretation under equilibrium sampling, when the concentration profile given by equation (7) is flat. The situation for other interpretations under equilibrium sampling is different.
The probability density to find a domain with diffusivity D by picking up a particle at random is proportional to the particles' density in the corresponding region and to probability to find a region with density n among all regions with given diffusivity, i.e. to np(n, D), where p(n, D) is the joint probability density of n(x) and D(x). Thus where N is the normalization constant, and, in the second line, p(n|D) is a probability density of the particle concentration conditioned on the diffusion coefficient in the corresponding domain. The normalization constant N is trerefore given by where n|D in the last line is the first conditional moment of the particle concentration. This gives us a physical interpretation of the normalization constant N : it is the inverse of the concentration averaged over all possible landscapes, which coincides with a volume mean n(x) in the thermodynamical limit. Therefore, it is useful to introduce the normalized equilibrium concentration at a position x and put down equation (14) as (note that N n = ν and p(n|D)dn = p(ν|D)dν). According to equation (7) the connection between n(x) and D(x) is deterministic, and the conditional probability density p(ν|D) is given by Subsituting equation (18) into (16) we get with the additional requirement that Equations (19) and (20) are easily solved by noting that the first one gives the form of p(D) up to the normalization constant D α−1 : where in the second line we simply substitute the expression for p S (D) as following from equation (12).
Requiring the normalization of the lhs we get Substituting this expression into equation (21) we obtain and recognize on the rhs a Γ-distribution with the shape parameter β = β − α + 1 and with a mean Therefore the PDF of local diffusion coefficients is again given by a Gamma distribution. For future convenience we repeat the expressions for the parameters of the distribution of local diffusion coefficients as functions of the dimension of space for homogeneous sampling (p(D) = p S (D) with p S (D) given by equation (12)) and for equilibrium sampling, equation (24), in table 1.
The cumulative distribution function corresponding to the PDFs, equations (12) and (24), is with z = β or z = β respectively, and γ(z, x) being the lower incomplete Γ-function. This form will be used for generating diffusivity landscapes in simulations, as discussed section 5.
In the HK case and in the Ito case for homogeneous sampling the parameters are β = (d + 1)/2 and D = D 0 ; in the equilibrated Ito case they are β = (d + 3)/2 and D = d+3 dD for the case of homogeneous sampling we can put down the expression for ν(D) which will be useful when calculating the effective medium properties: For the case of equilibrium sampling D α−1 is given by equation (22), and

Long-time behavior: homogenization
Now let us consider the long-time behavior of the diffusion coefficient corresponding to its homogenization. At long times particles sample large domains of the system and feel some effective diffusion coefficient D * which in general differs from the mean diffusion coefficient as obtained from the short time behavior. Note that this effective diffusion coefficient is understood as a long time limit of a diffusion coefficient averaged over the ensemble of different realizations of the random medium.
The problem of finding effective large-scale characteristics of an inhomogeneous medium is an old one. The simplest approaches correspond to static, time-independent problems, and the best-investigated situations are pertinent to the homogenization of the electric conductance, and to mathematically similar cases of homogenization of dielectric or magnetic successibility, see e.g. [37]. In the language of conductivity the situation is as follows: one considers a large piece of medium (inhomogeneous conductor), say in a form of a slab, with given boundary conditions for the potential (for example, with two opposite sides connected to a battery by highly conducting electrodes which are thus kept at constant potential difference Δφ) and measures the total current through the system. Locally the current density follows the Ohm's law j(x) = σ(x)E(x) giving a linear connection between a solenoidal field j(x) (for which divj = 0) and a potential field E(x) = grad φ(x) (for which in the static case rot E = 0 holds). Demanded is the connection between the volume mean current density j = 1 V V j(x)dx and mean electric field The mean values of the current density and of the electric field are immediately connected with the total current through the system and to the potential difference provided by a battery. A review of methods for calculating σ * for a variety of cases is given e.g. in reference [38].
Standard procedures of finding σ * can be applied to finding the effective long-time diffusion coefficient in the HK case, but in cases with inhomogeneous equilibrium. The reason is that the standard phenomenological Fick's equation corresponding to the HK case, equation (8), can be considered as a combination of a local continuity equation ∂ ∂t p(x) = − div j(x), with j(x) being now the probability or particles' flux, and the linear relation between j(x) and an obviously potential field grad p(x), namely the first Fick's law j(x) = D(x) grad p(x) serving as an analogue of the Ohm's law. At long times the probability distribution spreads, and the process of diffusion slows down, so that ∂ ∂t p(x) = − div j(x) → 0, and the methods used in the time-independent electrical case can be safely applied.
In the cases corresponding to inhomogeneous equilibrium the methods cannot be applied immediately. Although the corresponding diffusion equations can still be considered as a combination of a local continuity equation ∂ ∂t p(x) = − div j(x) with some linear response law, the last one does not have the form of the Ohm's law. For example, for the Ito case, equation (9), we have j(x) = grad D(x)p(x).
The homogenization of the diffusion coefficient in a heterogeneous medium is still similar to the one for the electric conductance [37,38], but with an important difference discussed below.
Let us assume that we are able to calculate the effective conductance σ * of an inhomogeneous system with local conductances σ(x) and denote this as a special type of an average, the homogenization mean, by σ * = σ(x) H . Then the effective diffusion coefficient in the homogenized regime for the general diffusive case is given by the same type of the average: where κ(x) = D(x)ν(x). This result follows by considering the stationary flow through a large piece of the medium with concentrations kept constant at its boundaries [39]. An alternative approach can follow the lines of [40] where some specific situations were considered. Here we give a simple explanation, not following the lines of the proofs but giving a physical intuition. The explanation relies on the property aσ(x) H = a σ(x) H which is evident from the initial definition of σ * as connecting the mean current density with the mean electric field. Let us assume that the diffusing particles carry a charge q, and move in electric field E under the force f = qE. Then the local current density is j(x) = q 2 n(x)μ(x)E(x), with μ(x) = D(x)/k B T being the local mobility, giving σ(x) = q 2 D(x)n(x)/k B T. Now we calculate the large scale conductance σ * = σ(x) H and assume that it is connected with the homogenized diffusion coefficient and the mean particles' density by the same Nernst-Einstein relation σ * = q 2 D * n(x) /k B T. Canceling all constant prefactors we get D(x)n(x) H = D * n(x) , which is our equation (28).
The homogenized conductance possesses the upper and the lower bound, from which the Wiener bounds, see e.g. [37], are the most general ones [41]. In our case they read The averages in the bounds in equation (29) can be considered either as volume means or averages over landscapes. In d = 1 the homogenization mean corresponds exactly to the lower bound which for conductance follows immediately from the Ohm's law. Using the distribution of the local diffusion coefficients D given by equations (12) and (23), and the connections between ν and D given by equations (26) and (27) for the homogeneous and equilibrium sampling, respectively, we can evaluate the upper and the lower Wiener bounds for the effective diffusion coefficient D * .
Thus, for homogeneous sampling we get so that the upper bound on D * is given by D 0 only for α = 1, i.e. in the HK-case. In one dimension the upper bound however cannot be realized by virtue of equation (30). For all α < 1 one has D * < D 0 , and this inequality is strict, so that no BnG diffusion can be observed under any circumstances. The lower Wiener bound is given by the expression  For equilibrium sampling the upper Wiener bound is universally equal to D 0 , and the lower bound is given by the expression .

This corresponds to
in d = 2, and to in d = 3, respectively. The upper and the lower bounds coincide for α = 0 and the result D * = D 0 is in this case (i.e. for the Ito model) exact. Therefore, we know that for the equilibrated Ito model the BnG diffusion does definitely take place. Now we have to check whether the upper bound is realized in d = 2 and 3 for the HK model under homogeneous sampling, and for any of non-Ito models under equilibrium sampling. We first proceed to show that the last two possibilities are very unlikely, and then, by means of numerical simulations, that these are not realized, indeed.
Far from percolation transition σ * is typically well reproduced by the effective medium approximation (EMA), see [38] for the discussion. Within this approximation D * = κ H is given by the solution of the equation where the average again can be considered either as a volume average or as an average over the distribution of κ. We note that for d = 1 EMA reproduces the exact result, equation (30). Equation (32) is pertinent to quadratic and cubic lattices in d = 2 and 3 [42], and to continuous d-dimensional systems [38]. For known p(D) and ν(D) equation (32) takes the form and can be solved numerically by using the parameters of p(D) given in table 1 and ν(D) given by equations (26) and (27) for homogeneous and equilibrium sampling, respectively. The results of such numerical solution for d = 2 are shown in figure 1. We note that for the Ito cases, α = 0, equations (26) and (27)

Systems with homogeneous equilibrium
For systems with homogeneous equilibrium (the HK interpretation, or the random diffusivity model of reference [40]) ν = 1 and there is no difference between equilibrium and homogeneous sampling situations: the distributions of κ and of D coincide. The upper Wiener bound corresponds to D = D 0 . In d = 1 the first inverse moment ∞ 0 D −1 p(D)dD of the PDF equation (12) diverges, and the effective diffusion coefficient D * given by equations (28) and (30) vanishes, giving rise to anomalous diffusion [39]. In d = 2 and d = 3 the value of D −1 given by equation (13)  Independently on the quality of approximation given by EMA we note that the EMA result is realizable in the continuum case [43]: the ensemble of all 'disordered' configurations contains realizations with the effective conductance (diffusivity) equal to the one predicted by EMA. The lower Wiener bounds are also realizable [37]. Therefore situations leading to the BnG diffusion are highly unlikely. The ensemble of disordered systems should indeed contain the realizations with diffusivities close to the upper bound D 0 , but also the ones with considerably smaller diffusivities given by the EMA and by the lower bound. Therefore D * will typically be lower than D 0 .
In figure 2 we present the full time dependence of the diffusion coefficient in the HK interpretation as following from numerical simulations. The figure shows D(t) = 1 4 d dt |x| 2 (t) normalized on D 0 in d = 2. The details of our simulation approach are given in section 5. One readily infers that the diffusion coefficient decays with time, so that no BnG diffusion is observed. The value of the terminal diffusion coefficient D * obtained in simulations agrees well with the EMA prediction.

Systems with inhomogeneous equilibrium
For systems with inhomogeneous equilibrium the situations under homogeneous and equilibrium sampling are different. We start our discussion using the hints given by the EMA, as shown in figure 1. For homogeneous sampling, we see that the difference between the short time diffusion coefficient D 0 and the long-time one D * increases when α decreases from 1. Parallel to the discussion above, there is no reason to await the BnG behavior, which, for α < 1, is essentially forbidden by equation (31). The difference is maximal for the Ito case, when in d = 2 the terminal diffusion coefficient is exactly one third of the short-time one (in d = 3 it will make a half of an initial one).
The result of EMA for equilibrium sampling is shown in figure 1 as the upper curve. We see that the behavior is opposite to the one for homogeneous sampling, and that for α → 0 the difference between D * and D 0 vanishes, which result again does not depend on EMA. Moreover, this behavior persists even in d = 1. The simulation results for the Ito cases are shown in figure 2 along with the one for the HK.
The effective medium approximation only allows for comparing the initial and terminal values of the diffusion coefficients, but gives neither the full time-dependence of the MSD (and therefore of D(t)) nor the forms of the corresponding PDFs (essentially, no analytical method is known to reliably reproduce such PDFs in the intermediate time domain). Therefore here we have to rely on the results of numerical simulations. Figure 3 displays such PDFs for equilibrated Ito case at different times. These exhibit the

Simulations of the PDF and MSD
In our simulations we start from a discretized model and consider the situation described by the master equation where i and j number the sites of a square or cubic lattice with lattice constant a = 1. Only the transitions between neighboring sites are possible. This master equation describes a random walk scheme and can be considered as a spatial discretization scheme for the corresponding Fokker-Planck equations, Equation (8) or (9). The transition rates follow the distribution similar to the distribution of local diffusivities p(D) given by equation (24); the local diffusivity for the rates which vary slowly in space is simply D = a 2 w = w. For the HK case, equation (8), the rates w ij = w i←j satisfy the condition of the detailed balance, i.e., in the absence of the external force w ij = w ji as discussed in reference [33]. The rates follow from the PDF p(w) similar to p(D). To simulate the Ito situations, equation (9), we assume that the transition rates from each site to all its neighbors are the same.
The transition rate from the site j to any of its neighboring sites i is w ij = w j and this w j is distributed according to the corresponding p(w). The transitions now are asymmetric: w ij = w ji , which makes a difference. The first prescription corresponds essentially to a barrier model of a random medium, the second one-to the trap model [33].
To simulate the two situations we generate two or three dimensional arrays of correlated transition rates w ij . For the Ito cases all w ij = w j are defined on the sites j of a simple square or simple cubic lattice, corresponding to the lattice of sites at which the probabilities p j are defined. For the HK case the transition rates w ij = w ji are defined at the midpoints of bonds of the corresponding lattice of p j . In d = 2 this lattice of midpoints is a square lattice with the lattice constant equal to a/ √ 2 and with main axes rotated by π/4 with respect to the axes of the p j -lattice, but can also be considered as a quadratic lattice with lattice constant a with basis (with an additional site placed at a center of a square), which is shifted by a/2 with respect to the lattice of p j . For three-dimensional case the lattice of the midpoints of bonds of a simple cubic lattice is an octahedral lattice, again considered as a simple cubic lattice with basis. Note that the arrays used for simulation of HK and Ito cases are different in size: for example, for simulating a 2 d lattice with N = l × l sites we need l 2 different values of transition rates for the Ito cases and 2l(l − 1) (i.e. approximately twice as many) different rates for the HK case.
The correlated random variables w ij on the corresponding lattices of transition rates can be easily obtained by a probability transformation. Let us call F −1 β (y) the function inverse to F β (D), as given by equation (25). Then the probability transformation transforms the Gaussian variable z with zero mean and unit variance into a Γ-distributed w with shape parameter β and unit mean. The corresponding function F β (x) can be easily inverted (there exists a standard MATLAB implementation for this inverse), and therefore the corresponding fields can be easily simulated for any given two-point correlation function. In our simulations we use independent Gaussian variables for the uncorrelated case, or a correlated Gaussian landscape with Gaussian correlation function z i z j = exp(−r 2 ij /2λ 2 ) with r ij being the distance between the sites i and j on the corresponding lattice, with λ being the correlation length. Such a correlated Gaussian array is easily obtained by filtering of the initial array of independent Gaussian variables with their subsequent renormalization necessary to keep z 2 = 1 (note that initial arrays of independent Gaussian variables must be sufficiently larger than its 'internal' part used in simulations). The corresponding example of the diffusivity landscape as generated by this method is shown in figure 4.
The simulations of the situations under homogeneous sampling follow by numerical solution of the master equation for a particle starting at the origin. In 2 d the system is of the size (2L + 1) × (2L + 1) (i.e. one has −L k, l L) where L = 256 is used in simulations. The master equation is solved by forward Euler integration scheme to get p i (t) for each site i characterized by coordinates r i = (k, l). Figure 5 shows the exemplary distributions of p i (t) for the HK and Ito situations, which allow to grasp the differences between the cases. Note that the corresponding figures use different realizations of the landscape.
The probabilities p j (t) = p (k,l) (t) are then used for plotting the corresponding PDFs. The MSD for a particle starting at the origin (k = l = 0) is given by x 2 (t) = L k,l=−L (k 2 + l 2 )p (k,l) (t). The PDFs and MSDs are then averaged over different realizations of landscapes w ij (typically 1500 realizations). In the Ito case under equilibrium sampling the corresponding probabilities and MSDs are weighted with the inverse transition rate from the origin w −1 0 , which is proportional to D −1 (0) and therefore to the equilibrium concentration at the origin. Note that since the mean local diffusivity in this case is not equal to D 0 , additional normalization is applied to keep D 0 = 1. The PDFs p(x, t) shown in figures 3, 6 and 7 depict such PDFs for r i = (k, 0). The approach in 3 d is exactly the same (except for the fact that r i = (k, l, m)) but the size of the system is smaller: L = 64.
The results for time-dependent diffusivity D(t) are obtained by numerical differentiation of the corresponding MSD: D(t) = 1 2d d dt |x| 2 (t) . The results for HK case are shown in figure 8. Similar results for the Ito case are presented in figure 9. The results shown in figure 2 in the previous section are the curve for λ = 10 from figure 8, and the two corresponding curves for the homogeneous and eqilibrium sampling from figure 9. To show that the behavior in d = 3 is similar we present in figure 10 the results for the time dependent diffusion coefficient in d = 3 for λ = 3. The EMA prediction for the HK case corresponds to D * /D 0 ≈ 0.852, and for the Ito case under homogeneous sampling the analytical prediction is D * /D 0 = 0.5.
The results obtained via the direct solution are confirmed by stochastic simulations of continuous time random walks as described by the corresponding master equations using the Gillespy approach [49].

Discussion
The properties of diffusion in random diffusivity landscapes strongly depend on the model adopted, but still share some similarities. The PDFs averaged over the realizations of diffusivity landscapes show similar features in all cases: the gradual transition from exponential to a Gaussian form which happens at the wings of the distribution while all distributions still show a cusp at the origin at intermediate times. These findings are similar to what is seen in experiments [9][10][11]44], and in theoretical models of disordered systems like the trap model of references [45,46] close in spirit to our Ito model with homogeneous  sampling, or a barrier model of reference [26], close in spirit to the HK one. The behavior of the MSD in different systems however differs, i.e. may show the BnG diffusion, the crossover between different types of diffusive behavior, and even anomalous diffusion. As we have already seen, within the model class adopted (spatially inhomogeneous systems with slowly varying diffusion coefficient), the equilibrated Ito model is the only promising candidate for a model showing the BnG diffusion, i.e. the diffusion coefficient staying constant over the time.
The Ito interpretation relies on the martingale property, which, in the Gaussian case, means that the increments of the process during small time intervals are symmetric [47]. A random walk interpretation of this process can be a continuous time random walk with locally symmetric steps in space, in which the spatial change of the diffusivity is attributed to coordinate-dependent waiting times. Such a random walk scheme corresponds to a trap model [33] which is thus the most prominent candidate for modeling BnG diffusion. This random walk scheme with exponentially distributed waiting times at a site, with the mean proportional to the inverse of w j , indeed, exactly corresponds to the master equation we solve numerically in the Ito case. In higher dimensions (in d = 3 and in d = 2 approximately, up to logarithmic corrections which are hard to detect) random trap models may be mapped to uncorrelated CTRW under disorder averaging [48]. CTRW is a process subordinated to a simple random walk (in a continuous limit-a process subordinated to Brownian motion). In our case, the waiting times in our CTRW would however be correlated, which is different from the standard CTRW schemes. The diffusing diffusivity model [12] is also a representative of the class of models subordinated to Brownian motion, and shares some properties with the corresponding correlated CTRWs. We note however that this diffusing diffusivity model shows a different kind of transition from exponential to Gaussian PDF which does not lead to a cusp at the origin. This difference may be connected with the intrinsically mean-field character of the diffusing diffusivity approach which models the effects of quenched disorder by assuming the system to be spatially  homogeneous but possessing the corresponding distribution of waiting times (i.e. assuming transition rates to be random along the walker's trajectory). The fact that the transition rates are the same, when the particle returns to the same site of the lattice, is fully disregarded within this model.