2021

[1]Adrian Pacheco-Pozo

1]Institut für Physik, Humboldt-Universität zu Berlin, Newtonstraße 15,
D-12489, Berlin, Germany

2] IRIS Adlershof, Zum Großen Windkanal 2, D-12489, Berlin, Germany

# Random walks in correlated diffusivity landscapes

<EMAIL_ADDRESS>Igor M. Sokolov<EMAIL_ADDRESS>berlin.de [ [

###### Abstract

In recent years, several experiments 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.

###### keywords:

Disordered systems, Diffusion, Random walks, Correlations

##  1 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 Brown1828 . A. Einstein Einstein1905 was the
first to propose a mathematical description of this type of motion (see
Maiocchi1990 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 Balakrishnan2019 : (i)
The mean square displacement (MSD) of the particles from their initial
position grows linearly in time,

$\langle\mathbf{r}(t)^{2}\rangle=2dDt$ (1)

(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

$p(\mathbf{r},t)=\frac{1}{(4\pi
Dt)^{d/2}}\exp\left(-\frac{\mathbf{r}^{2}}{4Dt}\right).$ (2)

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 Perrin1916 . 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 Smoluchowski1906 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 Smoluchowski1916 . After Einstein and Smoluchowski, Langevin
Langevin1908 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 power-law time-dependence $\langle\mathbf{r}(t)^{2}\rangle\propto
t^{\gamma}$, with $0<\gamma<1$ (subdiffusion) or $1<\gamma<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 Metzler2000 ; Sokolov2012 ;
Hofling2013 ; Krapf2015 ; Oliveira2019 ; Wang2022 . 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 Wang2009 ; Leptos2009 ; Kurtuldu2011 ; Wang2012 ;
Skaug2013 ; Yu2013 ; He2013 ; Guan2014 ; Thorneywork2016 ; He2016 ;
Acharya2017 ; Wagner2017 ; Chakraborty2019 ; Kwon2019 ; Chakraborty2020 ;
Pastore2021 ; Pastore2022 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 Wang2012 . Some of the corresponding systems show a crossover
from the non-Gaussian distribution to a Gaussian one at long times Wang2009 ;
Pastore2021 . In several cases Wang2009 ; Wang2012 ; Skaug2013 ; He2016 ;
Wagner2017 ; Chakraborty2019 ; Kwon2019 ; Chakraborty2020 ; Pastore2021 ;
Pastore2022 , 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. Chubynsky2014 ; Jain2016 ; Chechkin2017 ; Lanoiselee2018_b ;
Slezak2018 . They assume slow random changes of the diffusion coefficient in
time. The particular variant of the model used in Chechkin2017 will be called
“the minimal model” of diffusing diffusivity in what follows. Another model
describing BnG diffusion is the diffusivity landscape model (DLM)
Postnikov2020 which considers that the diffusion coefficient varies slowly in
space.

The possible connection between the diffusing diffusivity and DLM was stated
in Ref. Chechkin2017 : 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 Luo2018 , 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 Luo2018 was immediately connected
to the correlated nature of disorder.

Recently, in Ref. Pacheco2021 , 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
$\mathbf{r}\to\mathbf{r}/\sqrt{t}$, $p\to 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.

##  2 Diffusivity landscape model

In what follows, we use the model proposed by Postnikov et al. Postnikov2020
which assumes the particles’ diffusion in a heterogeneous medium modeled by a
correlated diffusivity landscape $D(\mathbf{r})$. This motion is described by
the force-free Langevin equation with multiplicative noise

$\frac{d}{dt}\mathbf{r}=\sqrt{2D(\mathbf{r})}\;\boldsymbol{\xi}(t),$ (3)

with $\boldsymbol{\xi}(t)$ being a Gaussian white noise with
$\langle\boldsymbol{\xi}(t)\rangle=0$ and
$\langle\xi_{\mu}(t)\xi_{\nu}(t^{\prime})\rangle=\delta_{\mu\nu}\delta(t-t^{\prime})$
with $\mu,\nu$ representing Cartesian coordinates. This Langevin equation
corresponds to the Fokker-Planck equation

$\frac{\partial}{\partial t}p(\mathbf{r},t)=\nabla[(1-\alpha)\nabla
D(\mathbf{r})+D(\mathbf{r})\nabla]p(\mathbf{r},t)$ (4)

with $\alpha$ being the interpretation parameter taking values in the interval
${0\leq\alpha\leq 1}$ (see e.g. Arnoulx2023 for a comprehensive discussion).
The authors of Postnikov2020 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 ($\alpha=0$; then, of course, any other interpretation can be
used by introducing the corresponding deterministic force Arnoulx2023 ), 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 Wang2009 ; Leptos2009
; Kurtuldu2011 ; Wang2012 ; Skaug2013 ; Yu2013 ; He2013 ; Guan2014 ;
Thorneywork2016 ; He2016 ; Acharya2017 ; Wagner2017 ; Chakraborty2019 ;
Kwon2019 ; Chakraborty2020 ; Pastore2021 ; Pastore2022 , one can then show
that the single-point PDF of the diffusion coefficients in the corresponding
landscape should be given by a Gamma distribution:

$p(D)=\frac{\beta^{\beta}}{\Gamma(\beta)}\frac{1}{\overline{D}}\left(\frac{D}{\overline{D}}\right)^{\beta-1}\exp\left(-\beta\frac{D}{\overline{D}}\right),$
(5)

where $\Gamma(\cdot)$ is a Gamma function, and $\beta$ and $\overline{D}$ are
shape parameters dependent on the dimension of space.

In what follows, we will concentrate on the two-dimensional situation, for
which $\beta=5/2$ and $\overline{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
$\langle\mathbf{r}(t)^{2}\rangle=2dD_{0}t$ Postnikov2020 .

A finite-difference discretization of the Fokker-Planck equation, Eq. (4),
with $\alpha=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 $\alpha=0$, Eq. (4) can be rewritten in the form

$\frac{\partial}{\partial
t}p(\mathbf{r},t)=\Delta[D(\mathbf{r})p(\mathbf{r},t)],$

and its discrete version is

$\frac{d}{dt}p_{i}(t)=\sum_{k=1}^{4}\frac{D_{j_{k}}}{a^{2}}p_{j_{k}}(t)-\frac{4D_{i}}{a^{2}}p_{i}(t).$
(6)

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}\equiv
D(x_{i},y_{i})$, which are generated according to the following algorithm,
Ref. Postnikov2020 : 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 $\widehat{G}_{i}$ by applying the
Fourier filtering method Toral2014 to $G_{i}$. Like in Postnikov2020 , we take
the correlation function of the correlated field to follow

$\rho(\mathbf{r}_{ij})=\langle\widehat{G}_{i}\widehat{G}_{j}\rangle=\exp\left(-\frac{\mathbf{r}_{ij}^{2}}{2\lambda^{2}}\right),$
(7)

with $\lambda$ being the correlation length, and $\mathbf{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 checkerboard-like diffusivity
landscape.

Finally, the correlated Gaussian field $\widehat{G}_{i}$ is transformed into
the $\Gamma$-distributed diffusivity landscape $D_{i}$ by performing a
probability transformation:

$D_{i}=f(\widehat{G})=F^{-1}_{\beta}\left\\{\frac{1}{2}\left[1-\text{erf}\left(\frac{\widehat{G}}{\sqrt{2}}\right)\right]\right\\},$
(8)

where $\text{erf}(\cdot)$ is the error function and $F^{-1}_{\beta}(x)$ is the
inverse of the cumulative distribution function (CDF) $F_{\beta}(D)$ for the
PDF given by Eq. (5), which is given by

$F_{\beta}(D)=\int_{0}^{D}p(D^{\prime})dD^{\prime}=\frac{1}{\Gamma(\beta)}\gamma\left(\beta,\beta\frac{D}{\overline{D}}\right),$

with $\gamma(\cdot,\cdot)$ 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\times 256$ with $a=1$, $D_{0}=1$ and $\lambda=10$.

$0$
---
$100$
---
$x$
---
$200$
---
$0$
---
$100$
---

$y$
---

$200$
---
$1$
---
$2$
---
$3$
---
$4$
---
$5$
---
$6$
---

$D(x,y)$
---

Figure 1: A two-dimensional realization of the diffusivity landscape
$D(\mathbf{r})$ in the diffusivity landscape model. It corresponds to a
$256\times 256$ lattice with correlation length $\lambda=10$ and sampled
diffusion coefficient $D_{0}=1$.

Let us now return to our Eq. (6). Defining the transition rates as

$\omega_{i\to
j}=\frac{D_{i}}{a^{2}}=\frac{1}{4}\left(\frac{a^{2}}{4D_{i}}\right)^{-1}=\frac{1}{4}\frac{1}{\tau_{i}},$

with

$\tau_{i}=\frac{a^{2}}{4D_{i}}$ (9)

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

$\frac{d}{dt}p_{i}(t)=\sum_{k=1}^{4}\omega_{j_{k}\to
i}p_{j_{k}}(t)-\sum_{k=1}^{4}\omega_{i\to j_{k}}p_{i}(t),$

which can be rewritten as

$\frac{d}{dt}p_{i}(t)=\frac{1}{4}\sum_{k=1}^{4}\frac{1}{\tau_{j_{k}}}p_{j_{k}}(t)-\frac{1}{\tau_{i}}p_{i}(t).$
(10)

In Ref. Postnikov2020 , Eq. (10) was solved using the forward Euler method.
Here, we employ another approach. Like in Ref. Pacheco2021 , we use the fact
that the master equation (10) corresponds to a CTRW with exponential waiting
times distribution Klafter2011 . 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

$\psi(t\arrowvert\tau_{i})=\frac{1}{\tau_{i}}\exp\left(-\frac{t}{\tau_{i}}\right),$

with $\tau_{i}$ given by Eq. (9). Taking the lattice spacing to be the length
unit of the problem ($a=1$), we get

$\psi(t\arrowvert D_{i})=4D_{i}\exp\left(-4D_{i}t\right).$ (11)

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.

##  3 Space-time correlations

In the last decades, several correlated CTRW models were proposed which lead
to interesting behaviors Chechkin2009 ; Tejedor2010 ; Schulz2013 . 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 Pacheco2021 , which is not reproduced in such pre-averaged models
like the CTRW description of the DLM Pacheco2021 , or the minimal model of BnG
Chechkin2017 . 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
$\mathbf{r}\to\mathbf{r}/\sqrt{t}$, $p\to 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.

$(a)$
---
$(b)$
---
$\mathbf{r}_{0},t_{0}$
---
$\mathbf{r}_{0},t_{0}$
---
$\mathbf{r}_{1},t_{1}$
---
$\mathbf{r}_{2},t_{2}$
---
$\mathbf{r}_{3},t_{3}$
---
$\mathbf{r}_{4},t_{4}$
---
$\mathbf{r}_{1}^{\prime},t_{1}$
---
$\mathbf{r}_{2}^{\prime},t_{2}$
---
$\mathbf{r}_{3}^{\prime},t_{3}$
---
$\mathbf{r}_{4}^{\prime},t_{4}$
---
Figure 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.

Let us consider a particle whose initial position is $\mathbf{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 from the exponential distribution $\psi(t\arrowvert D_{0})$
given in Eq. (11), with $D_{0}=D(\mathbf{r}_{0})$. Next, the particle randomly
jumps to one of its neighboring lattice points whose position is
$\mathbf{r}_{1}$, and then waits for another time $t_{1}$ which is now drawn
from the exponential distributions $\psi(t\arrowvert D_{1})$, with
$D_{1}=D(\mathbf{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
$\\{\mathbf{r}_{0},\mathbf{r}_{1},\mathbf{r}_{2},\dots\\}$, which correspond
to a simple random walk, and a list $\\{t_{0},t_{1},t_{2},\dots\\}$ of the
corresponding waiting times between subsequent jumps which are drawn from the
exponential distributions $\psi(t\arrowvert D_{i})$ with
$D_{i}=D(\mathbf{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},\dots\\}$, 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.,
$\mathbf{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 $\mathbf{r}^{\prime}_{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
$\\{\mathbf{r}_{0},\mathbf{r}^{\prime}_{1},\mathbf{r}^{\prime}_{2},\dots\\}$
being sums of i.i.d. random steps, and the same list of waiting times
$\\{t_{0},t_{1},\dots\\}$ 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 $\mathbf{r}_{0}$ in the corresponding landscape.

$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$(a)\;t=10^{1}$
---
$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$(b)\;t=10^{2}$
---
$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$(c)\;t=10^{3}$
---
$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$(d)\;t=10^{4}$
---
Figure 3: A comparison of the one-dimensional cut of the PDF $q(\xi)=p(x,0)t$
of rescaled displacements $\xi=x/\sqrt{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\times 2048$ with $\lambda=10$
and $D_{0}=1$. Each realization contains $10^{5}$ particles. Plotted in Figure
3 is a cut of PDF $p(x,y\mathchar 24635\relax\;t)$ through the origin at
$y=0$. Moreover, following Pacheco2021 , we plot the PDF as a function of the
rescaled displacement $\xi=x/\sqrt{t}$. To keep the normalization of the PDF,
it has to be rescaled as $q(\xi)=t\cdot p(\xi)$. 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. Luo2018 and Pacheco2021 , 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.

##  4 Time-correlated continuous-time random walk

Ref. Pacheco2021 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

$\psi(t)=\int_{0}^{\infty}\psi(t\arrowvert
D)p(D)dD=\frac{5}{2}\left(\frac{3}{8}\right)^{5/2}\left(\frac{3}{8}+t\right)^{-7/2},$
(12)

for $D_{0}=1$. The corresponding waiting time density corresponds to a Pareto
type II distribution with mean waiting time $\langle t\rangle=1/4$ and second
moment $\langle t^{2}\rangle=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 Klafter2011

$\psi_{1}(t)=\frac{1}{\langle
t\rangle}\left[1-\int_{0}^{t}\psi(t^{\prime})dt^{\prime}\right]=\frac{3}{2}\left(\frac{3}{8}\right)^{3/2}\left(\frac{3}{8}+t\right)^{-5/2},$
(13)

which is also a Pareto type II distribution with a different exponent.

As we have shown in Pacheco2021 , 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 $\mathbf{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.

###  4.1 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(\mathbf{r})$ in terms of that of the correlated
Gaussian field $\widehat{G}$ defined in Eq. (7). The correlation function of
the diffusivity landscape is, by definition,

$C_{DD}(\mathbf{r})\equiv\xi(\mathbf{r})=\frac{\langle\delta
D(\mathbf{0})\delta D(\mathbf{r})\rangle}{\sigma_{D}^{2}}=\frac{\langle
D(\mathbf{0})D(\mathbf{r})\rangle-\overline{D}^{2}}{\sigma_{D}^{2}},$ (14)

with $\sigma_{D}^{2}$ the variance of the local diffusivity, and $\delta
D(\mathbf{r})=D(\mathbf{r})-\overline{D}$. In our case, $\overline{D}=5/3$ and
$\sigma_{D}^{2}=10/9$, for $D_{0}=1$ in two dimensions.

Let us turn our attention to the mean $\langle
D(\mathbf{0})D(\mathbf{r})\rangle$ in the last expression of Eq. (14). Just
for convenience, let us denote $D(\mathbf{0})$ and $D(\mathbf{r})$ as $D_{1}$
and $D_{2}$, respectively. By doing so, one can write

$\langle D(\mathbf{0})D(\mathbf{r})\rangle=\langle
D_{1}D_{2}\rangle=\int_{0}^{\infty}\int_{0}^{\infty}dD_{1}dD_{2}p(D_{1},D_{2}\mathchar
24635\relax\;\mathbf{r})D_{1}D_{2}.$ (15)

Now, we make use of the invariance of the probability measures,

$dD_{1}dD_{2}p_{D}(D_{1},D_{2},\mathbf{r})=d\widehat{G}_{1}d\widehat{G}_{2}p_{G}(\widehat{G}_{1},\widehat{G}_{2}\mathchar
24635\relax\;\mathbf{r}),$

with

$p_{G}(\widehat{G}_{1},\widehat{G}_{2}\mathchar
24635\relax\;\mathbf{r})=\frac{1}{2\pi\sqrt{1-\rho(\mathbf{r})^{2}}}\exp\left[-\frac{\widehat{G}_{1}^{2}+\widehat{G}_{2}^{2}-2\widehat{G}_{1}\widehat{G}_{2}\rho(\mathbf{r})}{2(1-\rho(\mathbf{r})^{2})}\right]$

being the bivariate distribution of the correlated Gaussian field used in the
first stage of construction of the diffusivity landscape, with
$\rho(\mathbf{r})$ being the correlation function of this field given by Eq.
(7). We note that the $\mathbf{r}$-dependence in this expression is fully due
to the one of the correlation function $\rho(\mathbf{r})$, and concentrate
only on this $\rho$-dependence, introducing the function
$g(\widehat{G}_{1},\widehat{G}_{2}\mathchar
24635\relax\;\rho)=p(\widehat{G}_{1},\widehat{G}_{2}\mathchar
24635\relax\;\mathbf{r})$. Now, we can write Eq. (15) as

$\langle
D_{1}D_{2}\rangle=\int_{-\infty}^{\infty}\int_{\infty}^{\infty}d\widehat{G}_{1}d\widehat{G}_{2}g(\widehat{G}_{1},\widehat{G}_{2}\mathchar
24635\relax\;\rho)f(\widehat{G}_{1})f(\widehat{G}_{2}),$ (16)

with $f(\widehat{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 $\langle D_{1}D_{2}\rangle$ is a
function of $\rho$ only, and therefore passing to the correlation function of
diffusivity landscape, which differs from $\langle
D(\mathbf{0})D(\mathbf{r})\rangle$ by shift and rescaling, we see that
$\xi(\mathbf{r})=\xi[\rho(\mathbf{r})]$, and the dependence $\xi(\rho)$ 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(\widehat{G})$ around zero up to the fourth
order:

$f(\widehat{G})\approx
a_{0}+a_{1}\widehat{G}+a_{2}\widehat{G}^{2}+a_{3}\widehat{G}^{3}+a_{4}\widehat{G}^{4}+O(\widehat{G}^{5}),$

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 $\rho$,
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 $\rho$ we find

$\langle D(\mathbf{0})D(\mathbf{r})\rangle\approx
b_{0}+b_{1}\rho+b_{2}\rho^{2}+b_{3}\rho^{3}+b_{4}\rho^{4}+O(\rho^{5}),$

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
$\overline{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}(\mathbf{r})\approx\xi[\rho(\mathbf{r})]$ with the function

$\xi(\rho)=\frac{b_{1}\rho+b_{2}\rho^{2}}{b_{1}+b_{2}}=c_{1}\rho+c_{2}\rho^{2}.$
(17)

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\leq\rho\leq 1$ as compared to the result of high-precision numerical
integration.

The transformation $\xi(\rho)$ 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

$\rho(\xi)=\sqrt{\left(\frac{c_{1}}{2c_{2}}\right)^{2}+\frac{\xi}{2c_{2}}}-\frac{c_{1}}{2c_{2}}.$
(18)

Substituting the expression for $\rho(\mathbf{r})$, Eq. (7), into Eq. (17) we
get the approximation for the correlation function of the diffusivity
landscape:

$C_{DD}(\mathbf{r})\approx
c_{1}\exp\left(-\frac{\mathbf{r}^{2}}{2\lambda^{2}}\right)+c_{2}\exp\left(-\frac{\mathbf{r}^{2}}{\lambda^{2}}\right).$
(19)

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(\mathbf{r}\arrowvert n)$ of the displacements given the
number of steps $n$:

$C_{DD}(n)=\int d\mathbf{r}\;C_{DD}(\mathbf{r})\;f(\mathbf{r}\arrowvert n).$
(20)

$0$
---
$2$
---
$4$
---
$6$
---
$8$
---
$10$
---
$n$
---
$\times 10^{2}$
---
$-2.0$
---
$-1.6$
---
$-1.2$
---
$-0.8$
---
$-0.4$
---
$0.0$
---

$\ln C_{DD}(n)$
---

Sim
---
Eq.(22)
---
Figure 4: Correlation function $C_{DD}(n)$ as a function of the number of
steps for the two-dimensional case with $\lambda=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(\mathbf{r}\arrowvert n)$ can be safely approximated by a two-
dimensional Gaussian distribution

$f(\mathbf{r}\arrowvert n)=\left(\frac{d}{2\pi
a^{2}n}\right)^{\frac{d}{2}}\exp\left(-\frac{d\,\mathbf{r}^{2}}{2a^{2}n}\right)=\frac{1}{2\pi\sigma^{2}n}\exp\left(-\frac{\mathbf{r}^{2}}{2\sigma^{2}n}\right),$
(21)

with $d=2$, $a=1$ and, respectively, $\sigma^{2}=1/2$. Within this
approximation, Eq. (20) takes the form

$C_{DD}(n)\approx
c_{1}\left(1+\frac{n}{2\lambda^{2}}\right)^{-1}+c_{2}\left(1+\frac{n}{\lambda^{2}}\right)^{-1}.$
(22)

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.

###  4.2 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 $\xi(n)=C_{DD}(n)$ given by
Eq. (22) we use Eq. (18) to obtain the correlation function $\rho(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 Toral2014 , we generate a correlated Gaussian vector
$\widehat{g}_{i}$ with correlation function $\rho(n)$ with $n=\arrowvert
i-j\arrowvert$. Using the probability transformation, Eq. (8), we transform
this correlated Gaussian vector into a one-dimensional array of diffusion
coefficients $\mathcal{D}_{i}$ with the desired correlation function
$C_{DD}(n)=\xi(n)$. The array of correlated diffusion coefficients
$\mathcal{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
$\psi(t\arrowvert\mathcal{D}_{i})=4\mathcal{D}_{i}\exp(-4\mathcal{D}_{i}t)$.
In each realization of the process, one repeats the procedure until obtaining
such a number $n^{\prime}$ of drawn elements that the sum of the first
$n^{\prime}$ elements does not exceed $t_{max}$ but the sum of the first
$n^{\prime}+1$ does. The number of elements $n^{\prime}$ 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(\mathbf{r},t_{max})=\langle f(\mathbf{r}\arrowvert
n^{\prime})\rangle_{n^{\prime}},$

with $f(\mathbf{r}\arrowvert n)$ being the PDF of displacements for a given
number of steps, Eq. (21), weighted with the waiting time of the first step.

$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$-1$
---
$0$
---
$1$
---
$10^{-1}$
---
$(a)\;t=10^{2}$
---

$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$-1$
---
$0$
---
$1$
---
$10^{-1}$
---
$(b)\;t=10^{3}$
---

Figure 5: A comparison of the one-dimensional cut of the PDF $q(\xi)=p(x,0)t$
of rescaled displacements $\xi=x/\sqrt{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\cdot 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\times 10^{6}$ different realizations of the correlated array of diffusion
coefficients $\mathcal{D}_{i}$, constructed with $\lambda=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\times 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\times 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.

##  5 The checkerboard model

$0$
---
$100$
---
$x$
---
$200$
---
$300$
---
$0$
---
$100$
---
$200$
---
$300$
---

$y$
---

$0$
---
$1$
---
$2$
---
$3$
---
$4$
---
$5$
---
6
---

$D(x,y)$
---

Figure 6: A two dimensional realization of the diffusivity landscape
$D(\mathbf{r})$ for the checkerboard model. It corresponds to a $300\times
300$ lattice where each cell has a size of $2\zeta$ with $\zeta=10$, and
sampled diffusion coefficient is $D_{0}=1$.

Let us now go a few steps back and consider how critical our assumption about
the shape of the correlation function $\rho(\mathbf{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\times N$
squares containing $2\zeta\times 2\zeta$ 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, $\zeta$ defines the correlation length of the landscape; to
compare to the results of the above DLM, we set $\zeta=\lambda$. Figure 6
shows one realization of the checkerboard-like diffusivity landscape on a
lattice of $300\times 300$ with $\zeta=10$ and $D_{0}=1$.

$-8$
---
$-4$
---
$0$
---
$4$
---
$8$
---
$\xi$
---
$10^{-4}$
---
$10^{-2}$
---
$10^{0}$
---

$q(\xi)$
---

$t=10^{1}$
---
$t=10^{2}$
---
$t=10^{3}$
---
$t=10^{4}$
---
$-1$
---
$0$
---
$1$
---
$10^{-1}$
---
Figure 7: A one dimensional cut of the PDF $q(\xi)=p(x,0)t$ of rescaled
displacements $\xi=x/\sqrt{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.

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\times 10^{4}$ different realizations of the
landscape, each one using $10^{4}$ particles. The landscape was constructed
with $N=109$ and $\zeta=10$, i.e., we consider a lattice of size $2180\times
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.

##  6 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 higher-order 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.

Acknowledgments

The work of A. P. P. was financially supported by Doctoral Programmes in
Germany funded by the Deutscher Akademischer Austauschdienst (DAAD) (Programme
ID 57440921).

## References

  *   * (1) Brown, R.: XXVII. A brief account of microscopical observations made in the months of June, July and August 1827, on the particles contained in the pollen of plants; and on the general existence of active molecules in organic and inorganic bodies. The Philosophical Magazine 4(21), 161–173 (1828). https://doi.org/10.1080/14786442808674769
  * (2) Einstein, A.: Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen. Annalen der Physik 322(8), 549–560 (1905). https://doi.org/10.1002/andp.19053220806
  * (3) Maiocchi, R.: The case of Brownian motion. The British Journal for the History of Science 23(3), 257–283 (1990)
  * (4) Balakrishnan, V.: Elements of Nonequilibrium Statistical Mechanics. Springer, Cham (2020)
  * (5) Perrin, J.B.: Atoms. Constable, London (1916)
  * (6) von Smoluchowski, M.: Zur kinetischen Theorie der Brownschen Molekularbewegung und der Suspensionen. Annalen der Physik 326(14), 756–780 (1906). https://doi.org/10.1002/andp.19063261405
  * (7) Smoluchowski, M.V.: Über Brownsche Molekularbewegung unter Einwirkung äußerer Kräfte und deren zusammenhang mit der verallgemeinerten Diffusionsgleichung. Annalen der Physik 353(24), 1103–1112 (1916) https://onlinelibrary.wiley.com/doi/pdf/10.1002/andp.19163532408. https://doi.org/10.1002/andp.19163532408
  * (8) Langevin, P.: Sur la théorie du mouvement Brownien. C. R. Acad. Sci. (Paris) (146), 530–533 (1908)
  * (9) Metzler, R., Klafter, J.: The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Physics Reports 339(1), 1–77 (2000). https://doi.org/10.1016/S0370-1573(00)00070-3
  * (10) Sokolov, I.M.: Models of anomalous diffusion in crowded environments. Soft Matter 8, 9043–9052 (2012). https://doi.org/10.1039/C2SM25701G
  * (11) Höfling, F., Franosch, T.: Anomalous transport in the crowded world of biological cells. Reports on Progress in Physics 76(4), 046602 (2013). https://doi.org/10.1088/0034-4885/76/4/046602
  * (12) Krapf, D.: Chapter five - mechanisms underlying anomalous diffusion in the plasma membrane. In: Kenworthy, A.K. (ed.) Lipid Domains. Current Topics in Membranes, vol. 75, pp. 167–207. Academic Press, New York (2015). https://doi.org/10.1016/bs.ctm.2015.03.002. https://www.sciencedirect.com/science/article/pii/S1063582315000034
  * (13) Oliveira, F.A., Ferreira, R.M.S., Lapas, L.C., Vainstein, M.H.: Anomalous diffusion: A basic mechanism for the evolution of inhomogeneous systems. Frontiers in Physics 7 (2019). https://doi.org/10.3389/fphy.2019.00018
  * (14) Wang, W., Metzler, R., Cherstvy, A.G.: Anomalous diffusion, aging, and nonergodicity of scaled Brownian motion with fractional Gaussian noise: overview of related experimental observations and models. Phys. Chem. Chem. Phys. 24, 18482–18504 (2022). https://doi.org/10.1039/D2CP01741E
  * (15) Wang, B., Anthony, S.M., Bae, S.C., Granick, S.: Anomalous yet Brownian. Proceedings of the National Academy of Sciences 106(36), 15160–15164 (2009) https://www.pnas.org/doi/pdf/10.1073/pnas.0903554106. https://doi.org/10.1073/pnas.0903554106
  * (16) Leptos, K.C., Guasto, J.S., Gollub, J.P., Pesci, A.I., Goldstein, R.E.: Dynamics of enhanced tracer diffusion in suspensions of swimming eukaryotic microorganisms. Phys. Rev. Lett. 103, 198103 (2009). https://doi.org/10.1103/PhysRevLett.103.198103
  * (17) Kurtuldu, H., Guasto, J.S., Johnson, K.A., Gollub, J.P.: Enhancement of biomixing by swimming algal cells in two-dimensional films. Proceedings of the National Academy of Sciences 108(26), 10391–10395 (2011). https://doi.org/10.1073/pnas.1107046108
  * (18) Wang, B., Kuo, J., Bae, S.C., Granick, S.: When Brownian diffusion is not Gaussian. Nature Materials 11(6), 481–485 (2012). https://doi.org/10.1038/nmat3308
  * (19) Skaug, M.J., Mabry, J., Schwartz, D.K.: Intermittent molecular hopping at the solid-liquid interface. Phys. Rev. Lett. 110, 256101 (2013). https://doi.org/10.1103/PhysRevLett.110.256101
  * (20) Yu, C., Guan, J., Chen, K., Bae, S.C., Granick, S.: Single-molecule observation of long jumps in polymer adsorption. ACS Nano 7(11), 9735–9742 (2013). https://doi.org/10.1021/nn4049039
  * (21) He, K., Babaye Khorasani, F., Retterer, S.T., Thomas, D.K., Conrad, J.C., Krishnamoorti, R.: Diffusive dynamics of nanoparticles in arrays of nanoposts. ACS Nano 7(6), 5122–5130 (2013). https://doi.org/10.1021/nn4007303
  * (22) Guan, J., Wang, B., Granick, S.: Even hard-sphere colloidal suspensions display Fickian yet non-Gaussian diffusion. ACS Nano 8(4), 3331–3336 (2014). https://doi.org/10.1021/nn405476t
  * (23) Thorneywork, A.L., Aarts, D.G.A.L., Horbach, J., Dullens, R.P.A.: On the Gaussian approximation in colloidal hard sphere fluids. Soft Matter 12, 4129–4134 (2016). https://doi.org/10.1039/C5SM03049H
  * (24) He, W., Song, H., Su, Y., Geng, L., Ackerson, B.J., Peng, H.B., Tong, P.: Dynamic heterogeneity and non-Gaussian statistics for acetylcholine receptors on live cell membrane. Nature Communications 7(1), 11701 (2016). https://doi.org/10.1038/ncomms11701
  * (25) Acharya, S., Nandi, U.K., Maitra Bhattacharyya, S.: Fickian yet non-Gaussian behaviour: A dominant role of the intermittent dynamics. The Journal of Chemical Physics 146(13) (2017). https://doi.org/10.1063/1.4979338. 134504
  * (26) Wagner, C.E., Turner, B.S., Rubinstein, M., McKinley, G.H., Ribbeck, K.: A rheological study of the association and dynamics of MUC5AC gels. Biomacromolecules 18(11), 3654–3664 (2017). https://doi.org/10.1021/acs.biomac.7b00809
  * (27) Chakraborty, I., Rahamim, G., Avinery, R., Roichman, Y., Beck, R.: Nanoparticle mobility over a surface as a probe for weak transient disordered peptide–peptide interactions. Nano Letters 19(9), 6524–6534 (2019). https://doi.org/10.1021/acs.nanolett.9b02764
  * (28) Kwon, T., Kwon, O.-S., Cha, H.-J., Sung, B.J.: Stochastic and heterogeneous cancer cell migration: Experiment and theory. Scientific Reports 9(1), 16297 (2019). https://doi.org/%****␣sn-article.tex␣Line␣1475␣****10.1038/s41598-019-52480-3
  * (29) Chakraborty, I., Roichman, Y.: Disorder-induced Fickian, yet non-Gaussian diffusion in heterogeneous media. Phys. Rev. Research 2, 022020 (2020). https://doi.org/10.1103/PhysRevResearch.2.022020
  * (30) Pastore, R., Ciarlo, A., Pesce, G., Greco, F., Sasso, A.: Rapid Fickian yet non-Gaussian diffusion after subdiffusion. Phys. Rev. Lett. 126, 158003 (2021). https://doi.org/10.1103/PhysRevLett.126.158003
  * (31) Pastore, R., Ciarlo, A., Pesce, G., Sasso, A., Greco, F.: A model-system of Fickian yet non-Gaussian diffusion: light patterns in place of complex matter. Soft Matter 18, 351–364 (2022). https://doi.org/%****␣sn-article.tex␣Line␣1525␣****10.1039/D1SM01133B
  * (32) Chubynsky, M.V., Slater, G.W.: Diffusing diffusivity: A model for anomalous, yet Brownian, diffusion. Phys. Rev. Lett. 113, 098302 (2014). https://doi.org/10.1103/PhysRevLett.113.098302
  * (33) Jain, R., Sebastian, K.L.: Diffusion in a crowded, rearranging environment. The Journal of Physical Chemistry B 120(16), 3988–3992 (2016). https://doi.org/10.1021/acs.jpcb.6b01527
  * (34) Chechkin, A.V., Seno, F., Metzler, R., Sokolov, I.M.: Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities. Phys. Rev. X 7, 021002 (2017). https://doi.org/10.1103/PhysRevX.7.021002
  * (35) Lanoiselée, Y., Grebenkov, D.S.: A model of non-Gaussian diffusion in heterogeneous media. Journal of Physics A: Mathematical and Theoretical 51(14), 145602 (2018). https://doi.org/10.1088/1751-8121/aab15f
  * (36) Slezak, J., Metzler, R., Magdziarz, M.: Superstatistical generalised langevin equation: non-Gaussian viscoelastic anomalous diffusion. New Journal of Physics 20(2), 023026 (2018). https://doi.org/10.1088/1367-2630/aaa3d4
  * (37) Postnikov, E.B., Chechkin, A., Sokolov, I.M.: Brownian yet non-Gaussian diffusion in heterogeneous media: from superstatistics to homogenization. New Journal of Physics 22(6), 063046 (2020). https://doi.org/10.1088/1367-2630/ab90da
  * (38) Luo, L., Yi, M.: Non-Gaussian diffusion in static disordered media. Phys. Rev. E 97, 042122 (2018). https://doi.org/10.1103/PhysRevE.97.042122
  * (39) Pacheco-Pozo, A., Sokolov, I.M.: Convergence to a Gaussian by narrowing of central peak in Brownian yet non-Gaussian diffusion in disordered environments. Phys. Rev. Lett. 127, 120601 (2021). https://doi.org/10.1103/PhysRevLett.127.120601
  * (40) de Pirey, T.A., Cugliandolo, L.F., Lecomte, V., van Wijland, F.: Path integrals and stochastic calculus. Advances in Physics 0(0), 1–85 (2023). https://doi.org/10.1080/00018732.2023.2199229
  * (41) Toral, R., Colet, P.: Stochastic Numerical Methods: An Introduction for Students and Scientists. Wiley-VCH, Weinheim (2014)
  * (42) Klafter, J., Sokolov, I.M.: First Steps in Random Walks: From Tools to Applications. Oxford University Press, Oxford (2011). https://books.google.de/books?id=y4BNxvImzTkC
  * (43) Chechkin, A.V., Hofmann, M., Sokolov, I.M.: Continuous-time random walk with correlated waiting times. Phys. Rev. E 80, 031112 (2009). https://doi.org/10.1103/PhysRevE.80.031112
  * (44) Tejedor, V., Metzler, R.: Anomalous diffusion in correlated continuous time random walks. Journal of Physics A: Mathematical and Theoretical 43(8), 082002 (2010). https://doi.org/10.1088/1751-8113/43/8/082002
  * (45) Schulz, J.H.P., Chechkin, A.V., Metzler, R.: Correlated continuous time random walks: combining scale-invariance with long-range memory for spatial and temporal dynamics. Journal of Physics A: Mathematical and Theoretical 46(47), 475001 (2013). https://doi.org/10.1088/1751-8113/46/47/475001