# From diffusion in compartmentalized media to non-Gaussian random walks

Jakub Ślęzak  Corresponding author<EMAIL_ADDRESS>Stanislav Burov
<EMAIL_ADDRESS>Physics Department, Bar-Ilan University, Ramat Gan 5290002,
Israel

###### Abstract

ABSTRACT. In this work we establish a link between two different phenomena
that were studied in a large and growing number of biological, composite and
soft media: the diffusion in compartmentalized environment and the non-
Gaussian diffusion that exhibits linear or power-law growth of the mean square
displacement joined by the exponential shape of the positional probability
density. We explore a microscopic model that gives rise to transient
confinement, similar to the one observed for hop-diffusion on top of a
cellular membrane. The compartmentalization of the media is achieved by
introducing randomly placed, identical barriers. Using this model of a
heterogeneous medium we derive a general class of random walks with simple
jump rules that are dictated by the geometry of the compartments. Exponential
decay of positional probability density is observed and we also quantify the
significant decrease of the long time diffusion constant. Our results suggest
that the observed exponential decay is a general feature of the transient
regime in compartmentalized media.

##  I Introduction

The Brownian motion is ubiquitous in applications, be it the microscopic
motion of molecules Feynman (1964), search patterns of animals Humphries  _et
al._ (2010), or the prices of financial options Merton (1971). Following the
classical argument of Einstein Einstein (1905), if there exists a time scale
in which the changes of the investigated variable can be treated as a result
of additive, independent and homogenous fluctuations with finite second
moment, then in larger time scales the observed process is the Brownian
motion. This insight was later formalised as the functional central limit
theorem Billingsley (1999).

The Brownian motion became a natural start for theoretical and experimental
investigations of more complicated stochastic models. For example, lowering
the requirement of finite moments led to the rich theory of Lévy flights
Shlesinger  _et al._ (1999); lowering the assumption of independence was one
of the cornerstones of the anomalous diffusion modelling Metzler and Klafter
(2000). The Langevin theory of diffusion investigates the time scales lower
than those required by the central limit theorem Coffey  _et al._ (1996) which
leads to the motions in which the deviations from Brownian motion appear in
the memory structure; the probability density of motion is still Gaussian,
only with a different scale.

In recent years the influx of experimental data proved the existence of a
robust class of systems exhibiting non-Gaussian, in particular exponential,
tails of the probability density Wang  _et al._ (2012); Metzler (2017);
Bhattacharya  _et al._ (2013), together with normal or anomalous mean square
displacement Hapca  _et al._ (2009); Lampo  _et al._ (2017); Wang  _et al._
(2009). This common phenomenon is called _Brownian, yet non-Gaussian
diffusion_ Wang  _et al._ (2012); Metzler (2017). Its raising prominence and
importance in understanding the biochemical nature of the transport stimulated
various attempts to provide some – at least effective – description. These
include works using variants of the Langevin equation that use
superstatististical Beck and Cohen (2003); Beck (2006, 2011); Ślęzak  _et al._
(2018) or diffusing diffusivity Chubynsky and Slater (2014); Lanoiselée  _et
al._ (2018); Chechkin  _et al._ (2017) approach.

It is commonly suspected that the true source of non-Gaussianity lays in the
heterogeneity of the medium Chechkin  _et al._ (2017); Ediger (2000); He  _et
al._ (2016); Wang  _et al._ (2012). It is also hypothesised that the correct
description may be related to random walks with traps Sokolov (2010a);
Postnikov  _et al._ . Yet, formulating and solving suitable models is still a
challenge far from completion. Here, we attempt to overcome it for systems in
which the diffusion is locally impeded by barriers, providing a general
mechanism for the appearance of the non-Gaussian diffusion; for a simplified
illustration of the type of medium considered see Fig. 1.

Such restrictions of the molecular motions are ubiquitous in nature Novikov
_et al._ (2011), especially in composite or porous materials Song  _et al._
(2000); Mair  _et al._ (1999). For example, according to the fences and
pickets model the cellular membrane is compartmentalized (likely by actin-
based membrane skeleton and various transmembrane proteins Fujiwara  _et al._
(2016); Sadegh  _et al._ (2017)), a fact essential in understanding its
structure and functions Chein  _et al._ (2019); Kusumi  _et al._ (2005). It is
also widely agreed that this phenomenon is crucial in determining the material
properties of the biological media, the transport of proteins, lipids and
their functions Syková and Nicholson (2008); Cory and Garroway (1990);
Yablonskiy  _et al._ (2002); Chein  _et al._ (2019).

Figure 1: A schematic representation of the diffusion in a two dimensional
compartmentalized medium. A particle (black trajectory) diffuses freely inside
the domains (white areas) but only rarely crosses the thick divisions between
them (grey area) which are thin compared to the average size of the domains
and act as barriers.

The transient confinement caused by these obstacles should result in a motion
similar to a random walk, with a particle jumping between adjacent domains.
This phenomenon was experimentally observed and termed hop diffusion Krapf
(2018); Fujiwara  _et al._ (2002). Its characteristic property is that the
long-time macroscopic effective diffusion coefficient is greatly reduced
compared to the short-time in-compartment diffusion coefficient with their
ratio providing important insight into the molecular properties of the system.
This ratio was found to be between 10 and 100 in various experiments, see the
list in Murase  _et al._ (2004). Similar results were found for porous
materials Latour  _et al._ (1994); Sen (2003); Wang  _et al._ (2019) and for
colloidal particles near the glass transition; in the last case what they call
the cage effect was observed together with the exponential tails of the
probability density Weeks and Weitz (2002); Chaudhuri  _et al._ (2007); Weeks
_et al._ (2000).

In what follows, we first make a universal observation on how the confinement
results in the static form of non-Gaussianity (Section II.). The argument is
quite general and shows the time scales at which the increments of the
diffusing particle’s position exhibit exponential tails. It does not specify
the dynamics, but provides a single jump’s distribution, which is a basic
building block of the random walk models which we derive. We study them using
one dimensional model of the system. We start from reducing the local
diffusivity model to the diffusion equation with interface conditions (Section
III.), which we solve domain-wise in order to derive the transition times of
the jumps between domains. Consequently, we manage to describe the diffusion
using environment-dependent random walk (Section IV.). It can be effectively
approximated using much simpler model of continuous time random walk (CTRW),
yielding the formula which links the effective diffusion coefficient to the
barriers’ permeability (Section V.). Next, we decouple the correlations of
this process and reveal its broad-tail non-Gaussian nature using the kurtosis
and the logarithm of characteristic function (Section VI.). Finally, we
consider a more reductive but mathematically elegant approximating CTRW for
which we calculate the exact probability density; it exhibits exponential
tails characteristic for the Brownian, yet non-Gaussian diffusion. The
derivations presented in Sections V. and VI. show the mathematical origin of
the analytical predictions presented in Figs 8-11, but for general
understanding of the results the form of the reduced model (14) and the
leading exponential-tail behaviour (36) suffice. The overview of the subject
and the results is given in Section VII. We also provide a list of commonly
used notation. Our simulation code is attached in the supplemental material.

##  II  Non-Gaussianity stemming from confinement

A possible approach for understanding the prevalence of the Gaussian
distribution is to look at it through the notion of entropy. For all different
variables $X$ in free space with a given variance, characterised by the
probability density function (PDF) $p_{X}$, the ones with Gaussian
distributions (with different means) maximise the entropy
$-\expectationvalue{\ln p_{X}(X)}$. Thus, if they are no other constrains
(such as potential field), any coordinate coupled to a heat reservoir will
converge to a Gaussian state Ebeling and Sokolov (2005).

This naturally opens up the question in what conditions the Gaussianity is
_not_ to be expected. Following the entropic argument, the simplest such case
is any bounded domain in which the entropy is maximized by the uniform, not
Gaussian distribution. Physically speaking, it corresponds to particles being
confined in some finite area by a reflecting barrier. Such a perfect local
confinement excludes the diffusion in the macroscopic scale, however in a more
detailed setting one can think about a system with dynamics dominated by two
time scales: the relaxation time $t_{r}$ of reaching the domain-wise
stationary uniform state, and the average escape time $t_{e}\gg t_{r}$ of
passing through the barrier. In this setting, for $t\ll t_{r}$ the scale of
the dynamics is too small to be affected by the confinement. For $t\gg t_{e}$
as long as there are no “hard” traps it is expected that the system will be
again homogenised and Brownian (even for the cases when the intra-domain
dynamics was not, e.g. was subdiffusive). For the time range in-between, the
motion is dominated by the confinement in one ($t\approx t_{r}$) or few
($t\approx t_{e}$) domains and always non-Gaussian. This time scale will be in
the centre of our discussion.

For $t<t_{e}$ the situation is simple: for a given particle its distribution
will be a uniform one with the shape reflecting the domain it was found in.
For the whole ensemble the distribution will be a mixture of uniform
distributions reflecting the domain shapes all over the medium. As we can see
we are in the situation in which the observed ensemble distribution is non-
Gaussian and completely determined by the random geometry of the system, to
some extent independently form the details of the dynamics.

This geometry is clearly hard to determine for many real systems and may be
very complex, e.g. porous media are often fractal-like. However, in many media
such as cellar membranes it is reasonable to assume that the domains are
mostly spherical and they vary mostly in size. Furthermore, the simplest size
distribution to consider is the exponential one. The argument for this is the
lack of dependence. We consider one dimensional line in the medium and ask if
in any given $\differential{x}$ it crosses the barrier or not. If the crossing
in $\differential{x}$ does not affect the distribution of other crossings
(they do not “see” each other) they must be located according to the Poisson
point measure Kallenberg (2017) and the distances between any two subsequent
barriers on this line are exponential.

We are now ready to calculate the stationary PDF which will be observed for
the whole ensemble. Let $L$ denote the diameter of a given domain. By choosing
proper displacement units we can always assume that $\expectationvalue{L}=1$
and its distribution of the sizes is a simple $p_{L}(l)=\exp(-l)$. However,
this is not the proper PDF to average over as we need to account for the
average number of particles in each domain. If the domains are not correlated
to any kind of a trap, this number will be proportional to the volume of the
domain. This leads to the equilibrium PDF Cox (1962) $p_{L}^{\text{eq}}(l)=\
c_{d}l^{d}\exp(-l),c_{d}=1,1/2,1/6$ in dimensions $d=1,2,3$ respectively
(these are gamma distributions $\mathcal{G}am(d+1,1)$, see the Notation
section).

For a given domain with a fixed $L$ the distribution is uniform, for one
dimension it is a simple $p_{X}(x|L)=1/L,|x|<L/2$, which yields

$p_{X}(x)=\int_{|x|<l/2}\\!\\!\\!\differential{l}\frac{1}{l}p^{\text{eq}}_{L}(l)=\int_{2|x|}^{\infty}\\!\\!\differential{l}\mathrm{e}^{-l}=\mathrm{e}^{-2|x|}.$
(1)

This simple argument shows how the geometry of the system determines the
distribution of $L$ and is reflected in the stationary PDF of particles’
displacements, which is found to be the Laplace (also called “two sided
exponential”) distribution with mean 0 and variance $1/2$, symbolically
$\mathcal{L}ap(0,1/2)$. For two and more dimensions the procedure is slightly
more complicated: for each coordinate $x,y$ or $z$ we will observe only the
projection of the total multidimensional probability mass, see Fig. 2. For
spheres these are semicircular law for two dimensions, and parabolic law for
three dimensions (as for each $\differential{x}$ we squeeze onto it a circle
with surface $\pi((1/2)^{2}-x^{2})$). The projections and averaged PDFs are
shown in the first and second row of Table 1.

Figure 2: Projection procedure applied to an exemplary domain. All the probability mass in two dimensions (the oval shape up) is projected onto a given axis, i.e. moved down leaving lengths of the depicted arrows intact, resulting in one dimensional mass (schematically shown as grey area below); then it is renormalized. | 1D | 2D | 3D
---|---|---|---
$p_{X}(x|L)$ | $\frac{1}{L}$ | $\frac{4}{\pi L}\sqrt{1-\mathopen{}\mathclose{{}\left(\frac{2x}{L}}\right)^{2}}$ | $\frac{3}{2L}\mathopen{}\mathclose{{}\left(1-\mathopen{}\mathclose{{}\left(\frac{2x}{L}}\right)^{2}}\right)$
$p_{X}(x)$ | $\mathrm{e}^{-2|x|}$ | $\frac{4|x|}{\pi}K_{1}(2|x|)$ | $\frac{1}{4}(1+2|x|)\mathrm{e}^{-2|x|}$
$p_{\Delta X}(x)$ | $\sim\mathrm{e}^{-|x|}$ | $\sim\frac{8}{3\pi}\mathrm{e}^{-|x|}$ | $\sim\frac{6}{|x|}\mathrm{e}^{-|x|}$
Table 1: Confined one-dimensional PDFs for the particles trapped in 1, 2, 3D
spheres. The conditional (i.e. with fixed trap diameter $L$) PDFs $p_{X}(x|L)$
in the top row are Wigner $n$-sphere distributions. Middle and bottom rows are
the PDFs averaged over the trap sizes, of the displacement $X$ and the
increments $\Delta X\mathrel{\mathop{:}}=X_{t+\Delta t}-X_{t},\Delta t\gg
t_{r}$ respectively. By $K_{1}$ we denote the modified Bessel function of the
second kind (see (DLMF, , Eq. 10.32.8)) which has tails
$K_{1}(2|x|)\sim(\pi/x)^{1/2}\exp(-2|x|)/2$. In the bottom row we show only
the asymptotics, but in the case of 1D and 3D the exact values of $p_{\Delta
X}$ can be calculated by direct integration, the formulas are complicated but
may be expressed using special function Ei and polynomials.

The important thing to notice there is that they all decay like $\exp(-2|x|)$
multiplied by some power-law factor. Integrating by parts shows that this a
more universal behaviour, the requirements for its occurrence are only that
the PDF of the sizes has a tail $\propto l^{\alpha}\exp(-l)$ and the
stationary PDF within a fixed domain has a cut-off at $\pm L/2$. However,
there is a simplifying assumption in (1) which we did not immediately discuss:
the centre of the domain is located at $x=0$, whereas in typical experimental
data the position of the particle at the start of the measurement is taken as
0. This can be fixed with a litte effort and the corrected PDF still has
exponential tails, see the footnote 111 In this more realistic scenario the
left end of the initial domain lays at $\Theta L$ and the right one at
$(\Theta-1)L$; variable $\Theta$ has the uniform distribution
$\mathcal{U}nif(0,1)$. Calculation analogical to (1) yields incomplete gamma
PDF $p_{X}(x)=\Gamma(0,|x|)/(2|x|)$ which has tails $\sim\exp(-|x|)/(2|x|)$,
they are even thicker than $\exp(-2|x|)$. On the other hand this distribution
has logarithmic singularity at $x=0$ which make the motion more constrained at
short distances..

The limitation of the practical usefulness of the insights made above is that
the stationary distributions are visible in experiments only if many data
points within the range $t_{r}<t<t_{e}$ can be measured. However, we can
extend our approach if the sampling rate of the measurements is sufficiently
smaller than the escape time, $t_{r}\ll\Delta t<t_{e}$. Because $\Delta
t<t_{e}$ the subsequent values $X_{t}$ and $X_{t+\Delta t}$ will most often
still be inside the same domain. Consequently, the series of increments
$\Delta X_{t}\mathrel{\mathop{:}}=X_{t+\Delta t}-X_{t}$ will consist of long
intervals of the in-domain differences only rarely interrupted by the jumps to
the adjoining domains; these rare spikes can be then neglected. Because
$t_{r}\ll\Delta t$, the values $X_{t}$ and $X_{t+\Delta t}$ will be
independent (as the particle will transverse the domain multiple times during
$\Delta t$), so its PDF will determined by the difference of two independent
variables with $p_{X}(x|L)$ distribution: this is triangular distribution in
one dimension, complicated but elementary polynomial distribution in three
dimensions and a rather involved one in two dimensions (this is because the
convolution of two $\sqrt{1-x^{2}}$ functions is not elementary). In any case
their tails can be derived using integration by parts, see the third row of
Table 1. Using the ergodicity of $\Delta X_{t}$, this stationary distribution
can be estimated using even a single long trajectory, no ensemble averaging is
required. For all the cases this PDF also exhibits exponentially decaying
tails.

The limit of this line of arguments is the time scale $t>t_{e}$ in which the
dynamics becomes important. We need to know how long it takes for the particle
to jump from one domain to another and for this a more detailed model of the
barriers is required. This is the subject of the following sections.

##  III Transiently confined diffusion with locally Brownian dynamics

We imagine the heterogeneous medium as the collection of domains with regular
shapes inside which the particle moves relatively freely, separated by narrow
structures composed of a thick material (e.g. actin meshwork in the case of
plasma membranes Fujiwara  _et al._ (2016); Sadegh  _et al._ (2017); Chein
_et al._ (2019); Kusumi  _et al._ (2005)) which impedes the diffusion within
(these are grey areas in Fig. 1). We model this behaviour by making the
diffusion coefficient locally small, which results in the stochastic
differential equation

$\differential{X}_{t}=\sqrt{2D(X_{t})}\differential{B_{t}}.$ (2)

governed by the Brownian increments $\differential{B_{t}}$. The local
diffusivity function $D(x)$ appearing here is a combination of two microscopic
kinematic parameters: mean free path and correlation time. As such, the
equation by itself is physically ambiguous and requires interpretation. It was
previously established that if $D(x)$ is to describe an environment with
barriers (and no traps) present, the proper choice is kinetic
Hänggi–Klimontovich interpretation in which the PDF of the diffusion solves
the Fokker-Planck equation Sokolov (2010b)

$\partialderivative{t}p_{X}(x;t)=\nabla\cdot\big{(}D(x)\nabla
p_{X}(x;t)\big{)}.$ (3)

In a bounded domain with volume $V$ the constant PDF $p_{X}(x)=1/V$ is the
unique stationary solution ($\partial p_{X}/\partial t=0$) of this equation,
so this, and only this, interpretation agrees with our assumption of “no other
constraints” in the medium and consequently leads to the maximal entropy
distribution being uniform.

Up to this moment our considerations were quite general, from now on we will
limit ourselves to the one dimensional system, which will allow us to simplify
the geometry immensely and consequently obtain quite straightforward
description of the dynamics. In one dimension the barriers are thin $\Delta x$
intervals with a small local diffusivity $D_{b}\ll D$. Let a barrier be a
layer starting at $x_{k}$ and ending at $x_{k}+\Delta x$. Preservation of the
number of particles forces the flux $D(x)\partial p_{X}/\partial x$ to be a
continuous function, in particular at the borders of the barrier

$\displaystyle D\partialderivative{x}p_{X}(x_{k}^{-};t)$
$\displaystyle=D_{b}\partialderivative{x}p_{X}(x_{k}^{+};t),$ (4)
$\displaystyle D_{b}\partialderivative{x}p_{X}(x_{k}+\Delta x^{-};t)$
$\displaystyle=D\partialderivative{x}p_{X}(x_{k}+\Delta x^{+};t).$

When $D_{b}\ll D$ these conditions cause $p_{X}$ to change rapidly inside the
barrier, though this function must still be continuous. For the flux to be
also continuous the derivative $\partial p_{X}/\partial x$ must be
discontinuous at $x$ and $x+\Delta x$, but with $p_{X}$ still being smooth
inside the barrier. Now, we may express $p_{X}$ on the outside edges of
barrier using the values inside, then make a Taylor approximation and apply
(4)

$\displaystyle p_{X}(x_{k}+\Delta x^{+};t)-p_{X}(x_{k}^{-};t)$
$\displaystyle=\ $ $\displaystyle p_{X}(x_{k}+\Delta
x^{-};t)-p_{X}(x_{k}^{+};t)$ $\displaystyle=\ $ $\displaystyle\Delta
x\partialderivative{x}p_{X}(x_{k}^{+};t)+\order{\Delta x^{2}}$
$\displaystyle=\ $ $\displaystyle D\frac{\Delta
x}{D_{b}}\partialderivative{x}p_{X}(x_{k}^{+};t)+\order{\Delta x^{2}}.$ (5)

If we go to the limit $\Delta x\to 0$ and $D_{b}\to 0$ in such a way that
$D_{b}/\Delta x\to\kappa D$ we end up with a diffusion determined by the heat
equation within the domains

$\partialderivative{t}p_{X}(x;t)=D\partialderivative[2]{x}p_{X}(x;t),\quad
x\neq x_{k}$ (6)

which breaks at the locations of barriers where two interface conditions are
linking the values of the PDF on their left and right side. The first one is
again the flux continuity, the second one is the reduced form of (III)

$\displaystyle\partialderivative{x}p_{X}(x_{k}^{-};t)$
$\displaystyle=\partialderivative{x}p_{X}(x_{k}^{+};t),$ (7)
$\displaystyle\partialderivative{x}p_{X}(x_{k};t)$
$\displaystyle=\kappa\big{(}p_{X}(x_{k}^{+};t)-p_{X}(x_{k}^{-};t)\big{)}.$

The parameter $\kappa$ here is the barrier permeability; mind that some
authors denote $\kappa^{\prime}=\kappa D$ as permeability instead. This
argument can be generalized to two or three dimensional domains in a
straightforward manner, see Sánchez-Palencia (1980). It is also worth adding
that the same interface conditions can be derived using the barriers modelled
by the bumps of potential, but the derivation is quite technical Mandrekar and
Pilipenko (2016). Mathematically speaking, this diffusion has a rather
peculiar PDF: by making barriers’ thickness negligible we caused it to exhibit
finite jump discontinuities at each $x_{k}$ but at the same time still has
continuous directional derivatives everywhere. It is smooth only in the
stationary state, which (in a finite space $V$) is uniform as required,
$p_{X}(x;t\to\infty)=1/V$.

Conditions (7) also appear in chemistry Erban and Chapman (2007); Andrews
(2009) and should come as no surprise: it is nothing but a stochastic version
of the _Newton’s law of cooling_. In the most typical form it states that the
heat flux at a boundary is proportional to the temperature difference,
$\partial Q/\partial t\propto-(T_{\text{in}}-T_{\text{out}})$. Combining it
with the Fourier’s law $\partial Q/\partial t\propto-\partial
T_{\text{in}}/\partial x$ we end up exactly with (7). For this reason its
commonly called “Newton” or “convection” boundary condition. However, our case
is far less typical since we look for the solution on both sides of each
barrier. This is in contrast to the thermal problems in which $T_{\text{out}}$
is most often taken to be a fixed state of the environment. For this reason
the mathematical difficulty in solving the stochastic variant increases
significantly.

If only one barrier is present, let it be at $x_{0}=0$, the solution is still
straightforward to obtain. One only needs to split the initial condition into
symmetric and antisymmetric terms,
$p_{X}(x;0)=p_{X}^{S}(x;0)+p_{X}^{A}(x;0)=(p_{X}(x;0)+p_{X}(-x;0))/2+(p_{X}(x;0)-p_{X}(-x;0))/2$.
Then the solution itself can be split into symmetric and antisymmetric terms
which evolve independently with reflective (Neumann), $\partial
p_{X}^{S}(0^{+},t)/\partial x=0$, and radiation (Robin), $\partial
p_{X}^{A}(0^{+},t)/\partial x=2\kappa p_{X}^{A}(0^{+},t)$ boundary conditions.
Both can be solved using multiple standard methods. This argument also helps
to understand physical meaning of the parameter $\kappa$. For $\kappa\to 0$
the boundary conditions converge to the purely reflective one and the particle
stops being able to transverse the barrier. For $\kappa\to\infty$ they reduce
to the condition $p_{X}^{A}(0^{+};t)=0$. But as $p_{X}^{A}$ is antisymmetric
by definition, it only forces $p_{X}^{A}$ to be a continuous function; the
interface conditions and the barrier disappear altogether. The in-between case
of a finite $\kappa$ is a partially reflective barrier or a semi-permeable
barrier; look at Fig. 3 for an illustration how the resulting trajectories
look like.

Figure 3: Exemplary simulated trajectory for $\kappa=1/50,D=1$ and
discretization step $\Delta t=0.01$. Barriers were located according to the
standard Poisson random measure.

Alas, for any system with more than one barrier we loose this convenient
symmetry 222In the overall valuable work of Powles _et al._ Powles  _et al._
(1992) a system with the regularly placed barriers at
$x_{k}\in\\{\ldots,-1,0,1,2,\ldots\\}$ is considered. They claim they obtained
the exact PDF for this particular case but no proof is provided and the result
may be doubted.. The heat equation and interface conditions (7) are linear,
but the dependence on the environment is non-linear, as expected for a
heterogeneous system. For any finite number of barriers using the Fourier or
Laplace transform in the position space reduces the problem to solving a
system of algebraic equations, but the solutions are quite complicated and
inverting the transform seems to be beyond reach for more than two barriers
present.

Still, there exists a convenient stochastic representation of this diffusion
which provides another layer of physical meaning and facilitates numerical
simulations. In the mathematically impressive series of works Antoine Lejay
derived such a representation for a system with one barrier Lejay (2006, 2016,
2018) which extends easily to the general case. If the space would be discrete
the obvious and correct representation would be a random walk for which at
each visit to the barrier there is a chance of passing or being reflected, see
e.g. Nagylaki (1976); Barton (2008). For continuous space this approach does
not work directly because Brownian trajectories are so irregular they pass any
threshold infinitely many times. The representation becomes valid if we
restate this model as “each time the particle spends a unit of time at a
barrier there is a chance of being reflected”. Now, as we decrease the
discretization mesh the occupation times converge to the _local time_ random
field
$\ell_{t}(x)\mathrel{\mathop{:}}=\int_{0}^{t}\differential{\tau}\boldsymbol{1}_{X_{\tau}}(x)$
Borodin (2017) (this is basically a histogram calculated from the sample
trajectory $\\{X_{t}\\}$). At the same time the discrete escape times which
have geometric distribution converge to the smooth exponential variables.

Indeed, Lejay showed that if we take a series of independent and identically
distributed (i.i.d.) random times
$\tau_{1},\tau_{2},\ldots\overset{d}{=}\mathcal{E}xp(\kappa)$ the process
which is the reflected Brownian motion on the one side of the barrier until
$\ell_{t}(x_{0})>\tau_{1}$, then the reflected Brownian motion on the other
side until $\ell_{t}(x_{0})>\tau_{1}+\tau_{2}$ and so on, has a PDF which
solves (6) and (7). This makes sense if we come back to the derivation: for a
barrier with a finite thickness $\Delta x$ the particle must explore the
inside of the barrier for the time long enough it will reach the other side.
For small $\Delta x$ it is no surprise the passing is a Markovian event, thus
the waiting time must be exponential.

The joint distribution of the reflected Brownian motion and its local time is
known, so the condition $\ell_{t}(x)>\tau_{1}$ can be directly implemented in
a stochastic simulation Lejay (2016). For multiple barriers system only the
one on the left and the one on the right side of the current domain are
important. Because the dynamics is Markovian and local, at each step it is
only required to look for the closest barrier and check if the particle
escaped through it in a given $\Delta t$ or not. This is the method used to
simulate the trajectory shown in Fig. 3.

##  IV From the transient confinement to a random walk

It is natural to suspect that the diffusion dominated by the transient
confinement is a type of random walk, which is even suggested by the term “hop
diffusion”. We will derive the random walk which corresponds to the model of
the medium described in Section II. Similar approach was used as a possible
explanation for Lévy flights Barthelemy  _et al._ (2008) in quenched
disordered media Burioni  _et al._ (2010, 2014). However, in contrast to our
system, for these phenomena the short scale motion is assumed to be ballistic
and the distances between barriers have a power law distribution.

In our case, the process can be imagined as a particle which at each given
time is localised according to a uniform distribution in the domain it is in,
but after waiting a random escape time it transitions to a uniform
distribution in one of the adjacent domains with the cycle starting again.

Let us consider one such domain with ends at $\pm L/2$. The escape time is
determined by the solution of the diffusion equation if we remove returns to
the domain, that is, after the escape event $\ell_{t}(\pm L/2)>\tau_{1}^{\pm}$
took place the particle becomes absorbed and cannot return. This is equivalent
of putting absorbing (Dirichlet) boundary conditions just outside the domain
walls. By instantly drawing out all the probability mass outside the domain
this procedure reduces the interface conditions (7) to the radiation boundary
conditions

$\partialderivative{x}p_{X}((\pm L/2)^{\mp};t)=\mp\kappa p_{X}((\pm
L/2)^{\mp};t).$ (8)

The resulting process is a type of _partially reflected Brownian motion_
Grebenkov (2006).

Using again the equivalence to thermal problems, the dynamics can be described
by two dimensionless parameters: the Fourier number $Dt/L^{2}$ regulating the
ratio of diffusive to transport motion and the Biot number $\kappa_{B}=\kappa
L$ describing the ratio of heat resistance inside to that of the surface.
Systems with small Biot number have uniform temperature, which is precisely
our assumption of the dominance of local equilibria, which can now be
specified to be the requirement that $\kappa_{B}\ll 1$. There is a nuance here
as $L$ and thus $\kappa_{B}$ are domain-dependent, but as the distribution of
$L$ has short tails, extreme values of $\kappa_{B}$ are improbable and we may
just require that the average Biot number is small,
$\expectationvalue{\kappa_{B}}\ll 1$.

By rescaling $t$ and $x$ the equation can be reduced to the dimensionless

$\displaystyle\partialderivative{t}u(x;t)$
$\displaystyle=\partialderivative[2]{x}u(x;t),$
$\displaystyle\partialderivative{x}u((\pm 1/2)^{\mp};t)$
$\displaystyle=\mp\kappa_{B}u((\pm 1/2)^{\mp};t)$ (9)

which we expand into a series of eigenfunctions

$u=\sum_{n=0}^{\infty}a_{n}\mathrm{e}^{-\lambda_{n}^{2}t}\cos(\lambda_{n}x)+\sum_{n=0}^{\infty}b_{n}\mathrm{e}^{-\beta_{n}^{2}t}\sin(\beta_{n}x).$
(10)

The eigenvalues solve $\tan(\lambda_{n}/2)=\kappa_{B}/\lambda_{n}$ and
$\cot(\beta_{n}/2)=-\kappa_{B}/\beta_{n}$. For small $\kappa_{B}$ we can
expand the trigonometric functions around their zeros and obtain

$\displaystyle\lambda_{n}$
$\displaystyle=2n\pi+\frac{2\kappa_{B}}{\sqrt{n^{2}\pi^{2}+2\kappa_{B}}+n\pi}+\order{\kappa_{B}^{2}},$
$\displaystyle\beta_{n}$ $\displaystyle=\lambda_{n}+\pi,\qquad
n\in\\{0,1,2,\ldots\\}.$ (11)

Due to our assumptions we expect the distribution to be mostly uniform, and it
truly is: there is a spectral gap between the smallest eigenvalue
$\lambda_{0}=\sqrt{2\kappa_{B}}+\order{\kappa_{B}^{2}}$ and all the rest, thus
the corresponding eigenfunctions decay much faster. After a brief relaxation
only the mode $\cos(\lambda_{0}x)=1-\order{\kappa_{B}}$ remains and then it
also decays with the slower rate $\exp(-2\kappa_{B}t)$. This is the PDF of the
escape times for small $\kappa_{B}$. After coming back to the general case
with all the physical constants explicitly present this distribution is found
to be $\mathcal{E}xp(2\kappa D/L)$. Observe that the average time spent in the
domain during a visit is proportional to its length $L$, therefore this
approximation preserves the globally uniform distribution of particles. We
illustrate how the resulting approximation looks like for multiple domains in
Fig. 4.

Figure 4: Schematic comparison of the original PDF (black lines) and the base
mode approximation (red lines) with gradients exaggerated; for small $\kappa$
the functions are close to a step function with jumps at the barriers (green
lines). The approximation has the shape $A_{k}(t)\cos(\sqrt{2\kappa/L_{k}}x)$
inside each domain, centred for the middle. The cosines are then approximated
by ones, so the PDF is given by the amplitudes $A_{k}(t)$. The random walk
approximation is equivalent to stating that $A_{k}$ evolve according to the
master equation $\differential A_{k}/\differential t=2\kappa
DL_{k}^{-1}(-A_{k}+A_{k-1}/2+A_{k+1}/2)$; one can also think about it as a
trap model Bouchaud and Georges (1990), though with traps having varying sizes
and exponential waiting times.

We are now able specify the random walk model. First, we need to fix the
locations of the domains $x_{k}$. In Section II. we noted that if they are
independent from each other hey should be drawn from the stationary (i.e.
translationally invariant) Poisson point measure.

Figure 5: Comparison between the Euler scheme simulation ($\Delta t=0.0005$)
of the diffusion process (6), (7) with $D=1,\kappa=1/50$, the random walk
approximation $X_{2\kappa Dt}$ in the quenched environment, and the annealed
version with the correlated jumps and waiting times. (The PDFs are
symmetrical, we only show the right half.) In the regime $t=10$ the mass
contained in the initial domain is dominating; for $t=50$ the majority of the
PDF corresponds to few jumps, $t=100$ is the start of the Gaussian regime (in
the semi-log space they are nearly parabolas).

For simulations one can just use the fact that the standard Poisson measure
(corresponding to the first barrier put at point $x_{0}=0$) converges to the
stationary one rather quickly. It suffices to put barriers one after another
until sufficiently large $x$ are reached (for domain size 1 few hundred is
more than enough) and then move the beginning of the coordinate frame there.
More elegant approach is to correct the initial domain. The proper choice was
described in the footnote Note (1), we draw
$L_{0}\overset{d}{=}\mathcal{G}am(2,1),\Theta\overset{d}{=}\mathcal{U}nif(0,1)$
and put the left end at $(\Theta-1)L_{0}$, the right end at $\Theta L_{0}$.
The rest of the domains is unaffected and they have lengths $L_{k}$ with
i.i.d. $\mathcal{E}xp(1)$ distributions.

The random walk starts at domain $k=0$, after
$T_{1}\overset{d}{=}\mathcal{E}xp(2\kappa D/L_{0})$ it jumps with probability
1/2 to $k=1$ domain or to $k=-1$ one. The jumping with the analogical domain-
dependent waiting times repeats until we reach the final domain $k_{f}$ for
which the sum of waiting times exceeds $t$. The position of the particle is
then drawn from the uniform distribution
$\mathcal{U}nif(x_{k_{f}},x_{k_{f}+1})$ where $x_{k}$ denotes the left end of
the $k$th domain.

Conveniently, this approximation decouples parameters $\kappa$ and $D$ from
the dynamics, now they only determine the timescale. We may standardise this
random walk by saying $X_{t}$ has $\mathcal{E}xp(1/L_{k})$ waiting times, the
physical units are can then be returned considering the process with rescaled
time $X_{2\kappa Dt}$.

This random walk depends on the local barriers’ placement (that is quenched)
which makes the analytical analysis difficult. Both the domain sizes and
waiting times are short-tailed distribution, so in this model the particle
explores the space relatively freely and it may be expected that the annealing
procedure should yield a good approximation. Actually, we will make another
simplification at the same time by removing the correlation between the
subsequent jumps and waiting times. For the quenched random walk each
transition moves the probability mass into a uniform distribution inside some
domain and is therefore correlated with the (domain size dependent) subsequent
escape time. If we think about the quenched random walk as jumping from the
middle to the middle of the domain, each transition has length
$L_{k}/2+L_{k\pm 1}/2$. One can just consider annealed random walk with this
property, but it can be considerately simplified if we just divide all the
transition lengths in half.

The annealed model is thus as follows: we draw the initial domain size as
$L_{0}\overset{d}{=}\mathcal{G}am(2,1)$ and the subsequent ones as i.i.d.
$L_{k}\overset{d}{=}\mathcal{E}xp(1)$. The corresponding waiting times are
drawn as $T_{k}\overset{d}{=}\mathcal{E}xp(2/L_{k})$ or equivalently
$T_{k}=E_{k}L_{k}/2,E_{k}\overset{d}{=}\mathcal{E}xp(1)$. Using them we
construct a CTRW, i.e. the random sum $\sum_{k=0}^{N_{t}}J_{k}$ generated by
jumps $J_{k}\mathrel{\mathop{:}}=\pm L_{k}/2$ and the counting process
$N_{t}\mathrel{\mathop{:}}=\\#\\{k\colon T_{1}+\ldots+T_{k}\leq t\\}$. This
process should approximate the centre of the domain the particle is in; to
obtain the final position we account for the last local equilibrium by adding
$\mathcal{U}nif((\Theta-1)L_{0},\Theta L_{0})$ if $N_{t}=0$ or
$\mathcal{U}nif(-L_{N_{t}}/2,L_{N_{t}}/2)$ if $N_{t}\geq 1$.

We show the positional PDFs of the simulated processes in Fig. 5. There we
have chosen $D=1$ which determines the timescale and
$\expectationvalue{L_{k}}=1$ which analogically fixes the position scale. The
comparison proves that the random walk models are quite successful at
capturing the bulk of the probability. There is a noticeable (but not huge)
difference in the rate of exponential decay of the annealed process, but it is
to be expected in this type of approximation.

The simulations also suggest that all the processes have similar Gaussian
limit. For the annealed CTRW the Brownian limit can be derived in a formal way
using the standard approach. We consider the rescaled process
$X_{ct}/\sqrt{c}$ and then push the parameter $c$ to infinity. In this limit
the correction at the last site gets squeezed to zero, so it may be ignored.
The rescaled counting process $N_{ct}/c$ converges to
$t/\expectationvalue{T_{k}}=2t$ (this is the renewal theorem) and the
underlying classical random walk
$S_{n}\mathrel{\mathop{:}}=\sum_{k=0}^{n}J_{k}$ converges to the Brownian
motion,
$S_{cn}/\sqrt{c}\to\sqrt{\expectationvalue{J_{k}^{2}}}B_{n}=B_{n}/\sqrt{2}$.
The CTRW $X_{t}$ is given by subordination of $S_{n}$ by $N_{t}$, so it
converges to $B_{2t}/\sqrt{2}\overset{d}{=}B_{t}$ Magdziarz  _et al._ (2015).
Important thing to notice here is because the counting process $N_{t}$
collapses to a deterministic function, the dependence between waiting times
$T_{k}$ and jumps $\pm L_{k}/2$ becomes irrelevant in the long time limit.

As a consequence, after accounting for all the physical constants, the long
time effective diffusion coefficient
$D_{\text{eff}}\mathrel{\mathop{:}}=\lim_{t\to\infty}\delta_{X}^{2}(t)/(2t)$
of the annealed process is

$D_{\text{eff}}=\expectationvalue{\kappa_{B}}D.$ (12)

This agrees with the former results which were established for the systems
with periodically placed barriers Tanner (1978); Powles  _et al._ (1992);
Crick (1970); Latour  _et al._ (1994); Grebenkov  _et al._ (2014). A
comparison of $D_{\text{eff}}$ observed in the simulations of the transiently
confined diffusion and the theoretical value (12) is shown in Fig. 6.

Figure 6: Time dependent diffusivity defined as
$D(t)\mathrel{\mathop{:}}=\delta_{X}^{2}(t)/(2t)$ calculated from the
numerical simulation of the transiently confined process (6), (7) with short
time $D=1$ and $\expectationvalue{L_{k}}=1$. For various permeabilities
$\kappa$ the diffusivity $D(t)$ starts at $D(0)=1$ and decays, reaching the
level close to the value of theoretical prediction (12) derived for the
annealed process. The range of $D_{\text{eff}}/D$ shown was suggested by the
experiments of the hop diffusion Murase  _et al._ (2004).

It is worth to note that this process is Fickian (has linear square
displacement), but the transition from short-time diffusion constant $D$ to
the long-time $D_{\text{eff}}=\expectationvalue{\kappa_{B}}D$ causes the mean
square displacement to be a convex function. This behaviour can be seen in
Fig. 7 and bears great resemblance to contemporary experimental results for
porous media Murase  _et al._ (2004); Wang  _et al._ (2019). The observed
shape is easy to be mistaken for subdiffusion (MSD $\propto
t^{\alpha},\alpha<1$); this “intermediate subdiffusive transport” Höfling and
Franosch (2013) suggest paying strong attention towards this possibility in
the applied works as it could be misleading when observed without the entire
six orders of magnitude shown on the time axis.

Figure 7: The MSD of the simulated transiently confined diffusion,
$D=1,\expectationvalue{L_{k}}=1,\kappa=1/50,\Delta t=0.001$ (blue line)
together with the short time prediction of the unconstrained Brownian motion
with $\delta_{X}^{2}(t)\sim 2Dt$ (red dashed line) and the long time
prediction (12) with $\delta_{X}^{2}(t)\sim 2\kappa Dt$ (green dashed line).

Figure 7 also indicates the time range in which we can expect the non-Gaussian
behaviour. For short times the particle did not yet explore its initial domain
and behaves like undisturbed Brownian motion with diffusivity $D$. For very
long times the particle crossed a large number of domains, the environment
becomes homogenised in this scale and the movements are again Brownian, but
with lowered diffusivity $D_{\text{eff}}$. But, in the range between, where
the MSD in Fig. 7 bends, the particle’s movements are strongly affected by the
local environment, which originates the non-Gaussianity.

##  V Bessel and gamma waiting times CTRWs

It seems reasonable that the dependence between waiting times and jumps may
not be a crucial aspect of the dynamics. We will remove it and analyse the
resulting CTRW revealing its non-Gaussian behaviour.

For a particle exploring a medium with no long time correlation sources (such
as strong traps) the dependence becomes irrelevant at long times and even at
short times its is much weaker than in models such as Lévy walks. For the
latter class, the amplitudes of jumps and waiting times are identical and the
joint distribution is degenerate Zaburdaev  _et al._ (2015). In our case it
has a smooth PDF $p_{J,T}(x,\tau)=\exp(-2|x|)\exp(-\tau/|x|)/|x|$ and the
correlation between the jump amplitudes and waiting times has a lower value
around 0.5.

Similarly, we may ignore the uniform distribution within the last domain and
say that the particle just stops in the middle of the last domain. Using the
same line of thought, we can say that the particle started in the middle of
the first domain which therefore has form $[-L_{0}/2,L_{0}/2]$. Conveniently,
it makes the distribution of the initial stationary state and of the
subsequent jumps the same, that is the Laplace distribution
$\mathcal{L}ap(0,1/2)$ derived in (1). It clearly introduces noticeable error
for very short times when there is a significant probability mass in the
initial domain. The uncorrelated CTRW PDF $p_{X}$ can be easily corrected by
the formula

$p_{X}^{\prime}(x;t)=\operatorname*{P}(T_{1}<t)p_{X}^{\text{stat}}(x;t)+\operatorname*{P}(T_{1}\geq
t)p_{X}(x;t),$ (13)

where $p_{X}^{\text{stat}}$ is the stationary PDF in the uncentred domain, see
the footnote Note (1). For the clarity of the presentation we will ignore this
correction further on, but we note it makes the agreement shown in Fig. (10)
better for the short times.

After all these steps, we end up with a CTRW

$X_{t}=\sum_{k=0}^{N_{t}}J_{k},\quad N_{t}=\\#\\{k\colon
T_{1}+\ldots+T_{k}\leq t\\}$ (14)

with jumps $J_{k}\overset{d}{=}\mathcal{L}ap(0,1/2)$ and independent waiting
times $T_{k}=E_{k}L_{k}/2$ . The sum starting from $k=0$ accounts for the
initial distribution. Conditioning and direct integration shows that $T_{k}$
have PDF which can be expressed by the Bessel function,
$p_{T}(\tau)=4K_{0}(\sqrt{8\tau})$. This distribution or the distribution of
$\sqrt{T_{k}}$ appear in various sources as the Bessel distribution or _K-
distribution_ Jakeman and Pusey (1978); Kotz  _et al._ (2001), we will use the
former term. Its peculiar property is that this PDF has logarithmic
singularity at $0^{+}$.

The obtained Bessel waiting times CTRW is simple enough that one can use
Montroll-Weiss formula to obtain expressions for the moments and PDF in
Laplace and Fourier-Laplace spaces Metzler and Klafter (2000). In particular
it shows that the MSD has a logarithmic cusp $\propto t\ln(1/t)$ at small
times $t\to 0^{+}$, see blue line in Fig. 8. However, we will not pursue this
route, instead we will use a method which circumvents the use of Laplace
transform and provides formulas which work globally with respect to $t$.

This is made possible if we approximate the waiting times distribution with a
similar one for which the PDF of the counting process can be more easily
managed. Counting processes with infinitely divisible waiting times are known
to be well-studied; from this class short tailed processes are commonly
modelled using gamma distribution (especially in finance Papapantoleon
(2008)). From those, $\mathcal{G}am(1/2,1)$ is remarkably close to the Bessel
PDF which is reflexed in very small Kolmogorov distance $\approx 0.0395$. To
put this number in a perspective, in standard hypothesis testing setting, we
would need samples with around 2000 observations to notice the difference.
Crucially, this distribution also has the same mean, which is necessary to
preserve the Brownian limit of the original process.

We remark that such a significant similarity seems to stem from deeper
properties of the distributions in question. Gamma distribution has wrong both
$x\to 0^{+}$ asymptotics ($1/\sqrt{x}$ instead of logarithmic) and
$x\to\infty$ asymptotic (exponential instead of $\exp(-\sqrt{x})$). However,
it nearly does not matter, as $\mathcal{G}am(\alpha,1)$ for $\alpha=1/2$
balances these two discrepancies so they cancel themselves out for the vast
majority of the probability mass. The simple value of the parameter
$\alpha=1/2$ seems to be a coincidence, careful numerical study suggests that
the value $\alpha\approx 0.494$ is marginally better. For any practical
purpose the difference is insignificant but it is worth to add that the rest
of the argument below does not depend on this particular value $\alpha=1/2$
but rather only on $\alpha$ being sufficiently close to 1.

We start with a counting process $N_{t}$ with $\mathcal{G}(1/2,1)$ waiting
times. Knowing the exact distribution of the sums $T_{1}+\ldots+T_{n}$ we can
determine the PDF of $N_{t}$ using the relation

$\displaystyle\operatorname*{P}(N_{t}\geq
n)=\operatorname*{P}(T_{1}+\ldots+T_{n}\leq t)$ (15)
$\displaystyle=\operatorname*{P}(\mathcal{G}am(n/2,1)\leq
t)=\widetilde{\Gamma}(n/2,t);$

here $\widetilde{\Gamma}$ is the regularised gamma function,
$\widetilde{\Gamma}(a,t)\mathrel{\mathop{:}}=\int_{t}^{\infty}\differential{s}s^{a-1}\exp(-s)/\Gamma(a)$.
From this we get PDF expressed as a forward difference

$\displaystyle p_{N}(n;t)$
$\displaystyle=\operatorname*{P}(N_{t}=n)=\operatorname*{P}(N_{t}\geq
n)-\operatorname*{P}(N_{t}\geq n+1)$
$\displaystyle=\widetilde{\Gamma}(n/2,t)-\widetilde{\Gamma}((n+1)/2,t).$ (16)

The PDF of the whole CTRW can be linked to $p_{N}$ if we condition it by the
number of jumps performed,

$\displaystyle p_{X}(x;t)=\sum_{n=0}^{\infty}p_{S}(x;n)p_{N}(n),$ (17)

where $S_{n}$ is, as before, a partial sum process,
$S_{n}\mathrel{\mathop{:}}=J_{0}+\ldots+J_{n}$. To obtain its PDF we represent
each jump as a difference of two independent exponential variables
$J\overset{d}{=}\mathcal{J}^{+}-\mathcal{J}^{-}$,
$\mathcal{J}^{\pm}\overset{d}{=}\mathcal{E}xp(2)$. Because of this their sum
can also be split into
$S_{n}=(\mathcal{J}_{0}^{+}+\ldots+\mathcal{J}_{n}^{+})-(\mathcal{J}_{0}^{-}+\ldots+\mathcal{J}_{n}^{-})=\mathcal{S}_{n}^{+}-\mathcal{S}_{n}^{-}$.
Now, we can use the fact that $\mathcal{S}_{n}^{\pm}$ have a simple
distribution, i.e. $\mathcal{G}am(n+1,2)$, and calculate $p_{S}$ as a
convolution

$\displaystyle
p_{S}(x;n)=\int_{-\infty}^{\infty}\\!\\!\differential{y}p_{\mathcal{S}^{+}}(y;n)p_{\mathcal{S}^{-}}(-(|x|-y);n)$
$\displaystyle=\int_{0}^{\infty}\\!\\!\differential{z}p_{\mathcal{S}^{+}}(|x|+z;n)p_{\mathcal{S}^{-}}(z;n)$
(18)
$\displaystyle=4\mathrm{e}^{-2|x|}\frac{4^{n}}{(n!)^{2}}\int_{0}^{\infty}\\!\\!\differential{z}(z+|x|)^{n}z^{n}\mathrm{e}^{-4z}$
$\displaystyle=\frac{1}{n!}\sum_{k=0}^{n}\frac{(2n-k)!}{k!(n-k)!4^{n-k}}|x|^{k}\mathrm{e}^{-2|x|}.$

Finally, (V) and (V) substituted into (17) form the exact series
representation of the PDF of CTRW with gamma waiting times and Laplace jumps.
It can be easy plotted and compared with the data. Alas, because of the
complicated form of (V) and (V) it is unwieldy for analytical investigation.
To proceed, we will introduce another two approximations. The first yields the
Fourier transform of the PDF and moments.

We will exploit the fact that the Poisson counting process has a simple PDF
$t^{n}\exp(-t)/n!$ and survival function $\widetilde{\Gamma}(n,t)$. It has the
same shape of graph as the distribution (15), i.e.
$\widetilde{\Gamma}(n/2,t)$, as we see it is only rescaled by 1/2. So, for the
PDF expressed as a forward difference between this function at $n+1$ and $n$
we can replace the result by the Poisson PDF rescaled by 1/2

$p_{N}(n;t)\approx\frac{1}{2}\frac{t^{n/2}}{\Gamma(n/2+1)}\mathrm{e}^{-t}.$
(19)

As mentioned previously, the same argument works also for more general
$\mathcal{G}(\alpha,\beta)$ waiting times for which the PDF would be rescaled
by $\alpha$ instead of $1/2$. What we did here is basically a perturbation of
the counting process around the Poisson one which states that as long as the
power law $t^{\alpha-1}$ at $t\to 0^{+}$ is not to far away from $\alpha=1$
the obtained distribution is the Poisson one rescaled.

This procedure does not preserve the normalisation of the PDF. Correcting it
leads to

$p_{N}(n;t)\approx\frac{t^{n/2}}{\Gamma(n/2+1)}\frac{1}{E_{1/2}\big{(}\sqrt{t}\big{)}},$
(20)

where $E_{1/2}$ is the Mittag-Leffler function; for this particular
parametrisation it can be also expressed as
$E_{1/2}(x)=\exp(x^{2})\operatorname{erfc}(-x)$.

The error we make in this approximation is a differentiation error; it depends
only on the smoothness of the underlying function. Therefore, it will become
smaller when $t$ increases and the distribution spreads. So for large times
this correction of normalisation is unimportant, but it helps for the small
times, making the approximation globally efficient.

To simplify the formulas, for the remainder of this section let us mark out
the initial condition, i.e. decompose $X_{t}=\mathcal{X}_{t}+X_{0}$. By
conditioning, we relate the distribution of $\mathcal{X}_{t}$ in the Fourier
space to the transform of the jumps’ PDF
$\widehat{p}_{J}(\omega)=1/(1+\omega^{2}/4)$,

$\displaystyle\widehat{p}_{\mathcal{X}}(\omega;t)=\sum_{n=0}^{\infty}\big{(}\widehat{p}_{J}(\omega)\big{)}^{n}p_{N}(n;t)$
(21)
$\displaystyle\approx\sum_{n=0}^{\infty}\frac{z^{n}}{\Gamma(n/2+1)}\frac{1}{E_{1/2}\big{(}\sqrt{t}\big{)}}=\frac{E_{1/2}(z)}{E_{1/2}\big{(}\sqrt{t}\big{)}},$

where for brevity we denoted
$z\mathrel{\mathop{:}}=\sqrt{t}/(1+\omega^{2}/4)$. This formula can be
directly compared against experimental data by using the sample estimate of
$\expectationvalue{\cos(\omega X_{t})}$.

We want to reveal the non-Gaussian nature of the process and a convenient form
to do is log-characteristic function (LCF)
$\zeta^{\omega}_{\mathcal{X}}(t)\mathrel{\mathop{:}}=-2\ln\expectationvalue{\cos(\omega\mathcal{X}_{t})}/\omega^{2}$.
It measures dispersion of the displacements like MSD but gives bigger emphasis
of the spread of probability bulk, so it is expected to be smaller than MSD
for processes with a broader tails and more peaky PDF (for any $\omega$)
Ślęzak  _et al._ (2019). Indeed, for

$\zeta^{\omega}_{\mathcal{X}}(t)\approx
t\frac{1+\omega^{2}/8}{(1+\omega^{2}/4)^{2}}+\ln\frac{\operatorname{erfc}(-z)}{\operatorname{erfc}\big{(}{-\sqrt{t}}\big{)}}$
(22)

the first linear term is dominating for the large times and the second one
$\sim\sqrt{t}\pi^{-1/2}/(1+\omega^{2}/4)$ is dominating for the short times;
both have scaling coefficients smaller than those of the MSD, which we
calculate as the second derivative of (21)

$\displaystyle\delta^{2}_{\mathcal{X}}(t)$
$\displaystyle=\expectationvalue{\mathcal{X}_{t}^{2}}=-\partialderivative{{}^{2}}{\omega^{2}}\widehat{p}_{\mathcal{X}}(\omega;t)\Big{|}_{\omega=0}$
$\displaystyle\approx
t+\pi^{-1/2}\frac{\sqrt{t}}{\operatorname{erfc}\big{(}{-\sqrt{t}}\big{)}}\mathrm{e}^{-t}.$
(23)

Variable $\mathcal{X}_{t}$ corresponds to $X_{0}=0$ initial condition, but to
account for the non-zero one we only need to add
$\expectationvalue{X_{0}^{2}}=1/2$ to the MSD and
$-2\ln\expectationvalue{\cos(\omega
X_{0})}/\omega^{2}=2\ln(1+\omega^{2}/4)/\omega^{2}$ to the LCF (it is as
expected always $\leq 1/2$). For an illustration of the behaviour of these
dispersion measures see Fig. 8. It is also interesting to note that the MSDs
shown there prove that in practice it is very hard to see the difference
between the logarithmic cusp (at $t\to 0^{+}$) of the Bessel waiting times
CTRW and the square-root cusp of the gamma waiting times CTRW whereas the
difference in asympotics might have suggested otherwise.

Figure 8: Comparison of the MSD and LCF with $\omega=1/2,1$ for simulations of
Bessel waiting times CTRW, gamma waiting times CTRW, and their analytical
approximations (22) and (V). LCF being lower than MSD shows that the
distribution has more spread out bulk than Gaussian.

Continuing the analysis of non-Gaussianity, it can be also described using
higher moments and one of particular interest is the excess kurtosis

$\mathcal{K}_{X}(t)\mathrel{\mathop{:}}=\frac{\expectationvalue{X_{t}^{4}}}{\expectationvalue{X_{t}^{2}}^{2}}-3=\frac{\expectationvalue{\mathcal{X}_{t}^{4}}-3\expectationvalue{\mathcal{X}_{t}^{2}}^{2}+\expectationvalue{X_{0}^{4}}-3\expectationvalue{X_{0}^{2}}^{2}}{\big{(}\\!\expectationvalue{\mathcal{X}_{t}^{2}}+\expectationvalue{X_{0}^{2}}\\!\big{)}^{2}}.$
(24)

For any Gaussian variable it is 0; variables with tails broader than Gaussian
are expected to have positive excess kurtosis (to be “leptokurtic”) and in
particular for any Laplace variable it equals 3. We get the fourth moment from
the fourth derivative,

$\displaystyle\expectationvalue{\mathcal{X}_{t}^{4}}$
$\displaystyle=\partialderivative{{}^{4}}{\omega^{4}}\widehat{p}_{\mathcal{X}}(\omega;t)\Big{|}_{\omega=0}$
$\displaystyle\approx
t(t+9/2)+3\pi^{-1/2}\frac{(t+1)\sqrt{t}}{\operatorname{erfc}\big{(}{-\sqrt{t}}\big{)}}e^{-t}.$
(25)

Substituting the calculated averages into (24) leads to a rather bulky, but
purely elementary formula for $\mathcal{K}_{X}(t)$, which is compared against
simulations in Fig. 9. This function starts from $\mathcal{K}_{X}(0)=3$ and
then monotonically decays, at the beginning with rate
$\differential\mathcal{K}_{X}/\differential t=6(1-4/\pi)$ and as time grows
the decrease becomes faster, reaching asymptotic $\sim 9/(2t)$. This behaviour
reflects the initial Laplace regime which transitions to the long-time
Gaussian relaxation.

Figure 9: Excess kurtosis calculated from the simulations of CTRW with Bessel
waiting times, gamma waiting times compared with the approximation (V),(V).
This type of shape is shows relaxation from Laplace to Gaussian distribution.

As a side remark, we note that the kurtosis of $\mathcal{X}_{t}$ diverges at
0, $\mathcal{K}_{\mathcal{X}}(0^{+})=\infty$ and then it converges to 0 in the
same manner as $\mathcal{K}_{X}(t)$. This happens because kurtosis measures
the broadness of the distribution’s tails and also its spikiness at $x=0$. In
this case the latter is the culprit: the initial condition $\mathcal{X}_{0}=0$
causes the distribution of $\mathcal{X}_{t}$ to have Dirac delta at $x=0$ as
$\mathcal{X}_{t}=0$ unless $T_{1}>t$. For a system like ours it is clearly an
unphysical artefact, which shows the importance of a proper choice of the
initial condition.

##  VI Compound Poisson approximation

Figure 10: Comparison of the correlated CTRW and further approximations.
Bessel and gamma waiting times CTRWs are nearly undistinguishable, all others
are also close. These random walks have the same Gaussian limit, as shown for
$t=8$.

The methods which we used in the former section were successful in revealing
the broad tails of the studied diffusion, but it would be beneficial to show
the exponential decay of the PDF directly 333We note that this route will
bring stronger results than using the Fourier transform (21) directly, as it
is unfortunately hard to inverse analytically. One can use the method of
steepest descent to uncover the tail behaviour of $p_{X}$, but the result does
not seem to be very practical, describing only a far away range not easily
available in experiments.. To achieve this we are introducing a second CTRW
approximation with a very simple memory structure which will yield a full
asymptotic series expansion of the PDF.

The insight that we use is a particular interpretation of the gamma waiting
times CTRW. The PDF of the gamma waiting times process agrees with the Poisson
counting process at even $n$. This is because a sum of any pair of
$\mathcal{G}am(1/2,1)$ waiting times has exponential waiting time,
$T_{k}+T_{k+1}\overset{d}{=}\mathcal{G}am(1/2,1)+\mathcal{G}am(1/2,1)\overset{d}{=}\mathcal{E}xp(1)$.
Essentially, if we ignore odd numbers of jumps, the process behaves exactly
like the one with exponential waiting times which makes double jump each time.
Previously, we accounted for the odd numbers of jumps by interpolating the PDF
with the Poisson formula.

The second possibility is to replace the counting process with double jumps by
the one with twice the intensity but only single jumps. This way the average
number of jumps remains the same (we just spread them), however we distort
their variability. For this reason the second approximation is expected to be
less accurate than the first one, but it makes up for it by providing useful
representation of the PDF.

Thus, we have come to consider the CTRW with Laplace jumps and the Poisson
counting process $\mathcal{P}oiss(2t)$. For a comparison of this process to
the other CTRWs considered in this work see Fig. 10. In the decomposition
$X_{t}=\mathcal{X}_{t}+X_{0}$ the term $\mathcal{X}_{t}$ now belongs to the
class of _compound Poisson processes_. These processes are very regular, being
Markovian, infinitely divisible and having independent increments. As one of
the consequences, $\mathcal{X}_{t}$ must have linear MSD, precisely
$\delta^{2}_{\mathcal{X}}(t)=t$. It means that in this approximation we
completely neglect the non-linearity of the MDS present in the more detailed
models. However even in those, the linear range was appearing quickly (even
when the motion was still highly non-Gaussian), see Fig. 8, so the error
caused by this is not huge.

The Fourier space representation of the compound Poisson PDF often has a sleek
form; in our case

$\widehat{p}_{\mathcal{X}}(\omega;t)=\exp(2t\mathopen{}\mathclose{{}\left(\frac{1}{1+\omega^{2}/4}-1}\right)).$
(26)

Again, we may use it to calculate excess kurtosis, which is a simple rational
function

$\mathcal{K}_{X}(t)=\frac{3t+3/4}{(t+1/2)^{2}}.$ (27)

This function decays to 0 even slower than in the case of the gamma waiting
times, $\differential\mathcal{K}_{X}/\differential t=0$ at $t=0^{+}$ and the
asymptotic decay is $\sim 3/t$, but overall the shape of this function is not
much different than before.

Now, to get the PDF of $X_{t}$ in the position space we will use a particular
representation available only for the compound Poisson processes. Instead of
dividing each jump into the difference of two variables like before, we
separate them into two categories corresponding to which are positive and
which are negative for a given trajectory; the same is made to the initial
condition. This can always be done, but only for the Poisson process the two
thinned counting processes obtained,
$N_{t}^{+}\mathrel{\mathop{:}}=\\#\\{k\colon T_{1}+\ldots+T_{k},J_{k}>0\\}$
and $N_{t}^{-}\mathrel{\mathop{:}}=\\#\\{k\colon
T_{1}+\ldots+T_{k},J_{k}<0\\}$ are independent and Poisson with the twice
smaller intensity, $N_{t}^{\pm}\overset{d}{=}\mathcal{P}oiss(t)$ (this is a
particular consequence of the process being Markovian). Thus, we may represent
the total displacements as a difference of two positive, independent
processes, $X_{t}=X_{t}^{+}-X_{t}^{-}$, where

$X_{t}^{\pm}=\sum_{k=0}^{N_{t}^{\pm}}\mathcal{J}_{k}^{\pm},\quad\text{i.i.d.\
}\mathcal{J}_{k}^{\pm}\overset{d}{=}\mathcal{|}J_{k}|\overset{d}{=}\mathcal{E}xp(2).$
(28)

We again use the fact that
$\mathcal{J}_{0}+\ldots+\mathcal{J}_{n}\overset{d}{=}\mathcal{G}am(n+1,2)$;
conditioning by the number of jumps leads to a surprisingly elegant formula
for the PDF of $X_{t}^{\pm}$ Buchak and Sakhno (2017),

$\displaystyle
p_{X^{\pm}}(x;t)=\sum_{n=0}^{\infty}p_{\mathcal{G}am(n+1,2)}(x)p_{N^{\pm}}(n;t)$
(29)
$\displaystyle=2\mathopen{}\mathclose{{}\left(\sum_{n=0}^{\infty}\frac{(2xt)^{n}}{(n!)^{2}}}\right)\mathrm{e}^{-2x}\mathrm{e}^{-t}=2I_{0}\big{(}2\sqrt{2tx}\big{)}\mathrm{e}^{-2x}\mathrm{e}^{-t}.$

From this we immediately determine the part of $p_{X}$ corresponding to the
series of jumps in a one direction, these far tails are given by
$\operatorname*{P}(X_{0}>0)p_{N^{-}}(0;t)p_{X^{+}}(x;t)=I_{0}(2\sqrt{2tx})\mathrm{e}^{-2x}\mathrm{e}^{-2t}$.
Of course the full PDF is significantly broader. As in (V) it is given by the
convolution

$\displaystyle
p_{X}(x;t)=\int_{0}^{\infty}\\!\\!\differential{z}p_{X^{+}}(|x|+z;t)p_{X^{-}}(z;t)$
(30)
$\displaystyle=\mathrm{e}^{-2|x|}\mathrm{e}^{-2t}\int_{0}^{\infty}\\!\\!\differential{z}I_{0}\big{(}2\sqrt{2t(|x|+z)}\big{)}I_{0}\big{(}2\sqrt{2tz}\big{)}4\mathrm{e}^{-4z}.$

Let us denote the integral on the right in the above by
$\mathcal{I}=\mathcal{I}(x;t)$. As Bessel function $I_{0}$ increases
monotonically, but slower than exponentially, $\mathcal{I}(x;t)$ is also
increases slower than exponentially with respect to $t$ and $x$. The
dominating factor of the tails is thus $\exp(-2|x|)$ which is completely
static and is inherited from a single jump PDF, the rest is some time
dependent factor and second order correction with respect to $x$. These are
important as they contain information about the dynamics. The rest of the
section is devoted to their derivation.

The integral $\mathcal{I}$ can be expanded into a series. One method is using
Taylor expansion of $I_{0}$; this leads to a formula similar to the one shown
in the previous section. However, more efficient expansion can be obtained if
we integrate by parts.

It is easier to explain the general method first: under mild regularity
assumptions on $f$, an integration with exponent can be expressed as the
infinite order differentiation operator

$\displaystyle\int_{0}^{\infty}\\!\\!\differential{z}f(z)a\mathrm{e}^{-az}=f(0)+\int_{0}^{\infty}\\!\\!\differential{z}\derivative{z}f(z)\mathrm{e}^{-az}$
$\displaystyle=f(0)+\frac{1}{a}\derivative{z}f(z)\Big{|}_{z=0}+\frac{1}{a}\int_{0}^{\infty}\\!\\!\differential{z}\derivative{{}^{2}}{z^{2}}f(z)\mathrm{e}^{-az}=\ldots$
$\displaystyle=\sum_{n=0}^{\infty}\frac{1}{a^{n}}\derivative{{}^{n}}{z^{n}}f(z)\Big{|}_{z=0}=\frac{1}{1-\frac{1}{a}\derivative{z}}f(z)\Big{|}_{z=0}.$
(31)

In our specific case we need to take
$f(z)=I_{0}\big{(}2\sqrt{2t(|x|+z)}\big{)}\times
I_{0}\big{(}2\sqrt{2tz}\big{)}$ and $a=4$. Function
$I_{0}\big{(}\sqrt{z}\big{)}$ is absolutely monotonic, so the obtained series
consists of positive terms, i.e. we divide the probability mass into Bessel-
like modes.

It can be further simplified if one notes that the derivatives of
$I_{0}\big{(}2\sqrt{2tz}\big{)}$ reduce to factors depending only on $t$. We
expand each derivative using the binomial formula

$\derivative{{}^{m}}{z^{m}}\big{(}f(z)g(z)\big{)}=\sum_{k=1}^{m}\binom{m}{k}\derivative{{}^{k}}{z^{k}}f(z)\derivative{{}^{m-k}}{z^{m-k}}g(z)$
(32)

and then rearrange the terms according to the order of derivative acting on
$I_{0}\big{(}2\sqrt{2t(|x|+z)}\big{)}$. The coefficient before the $n$th one
is

$\displaystyle\sum_{k=n}^{\infty}\binom{k}{n}\frac{1}{4^{k}}\derivative{{}^{k-n}}{z^{k-n}}I_{0}\big{(}2\sqrt{2tz}\big{)}\Big{|}_{z=0}\\!\\!=\sum_{k=n}^{\infty}\binom{k}{n}\frac{1}{4^{k}}\frac{(2t)^{k-n}}{(k-n)!}$
$\displaystyle=\frac{1}{4^{n}}{}_{1}F_{1}(1+n,1;t/2)=\frac{1}{4^{n}}\mathrm{e}^{t/2}\sum_{k=0}^{n}\frac{(n-k+1)_{k}}{(k!)^{2}}\mathopen{}\mathclose{{}\left(\frac{t}{2}}\right)^{k},$
(33)

where ${}_{1}F_{1}$ is Kummer’s function and $(x)_{k}$ is Pochhammer symbol.
Denoting the polynomial on the right by $q_{n}(t)$ the integral of interest
takes the form

$\displaystyle\mathcal{I}$
$\displaystyle=\mathrm{e}^{t/2}\sum_{n=0}^{\infty}\frac{q_{n}(t)}{4^{n}}\derivative{{}^{n}}{x^{n}}I_{0}\big{(}2\sqrt{2t|x|}\big{)}$
(34)
$\displaystyle=\mathrm{e}^{t/2}\sum_{n=0}^{\infty}\frac{q_{n}(t)}{4^{n}}\mathopen{}\mathclose{{}\left(\frac{2t}{|x|}}\right)^{n/2}I_{n}\big{(}2\sqrt{2t|x|}\big{)}.$

Alternatively, instead of using $I_{n}$, one can also repeatedly apply
equalities $\differential I_{0}(z)/\differential z=I_{1}(z)$, $\differential
I_{1}(z)/\differential x=I_{0}(z)-I_{1}(z)/z$ and express the result as a
mixture of $I_{0}$ and $I_{1}$ multiplied by powers of $t$ and $x$. In any
case, the Bessel modes decay with respect to $n$ and have prefactors of type
$t^{\alpha}/x^{\beta}$, which shows that this expansion may converge slowly
for large $t$ and small $x$ but will converge quickly for small $t$ and large
$x$.

Going back to the PDF under consideration, the expansion becomes

$\displaystyle
p_{X}(x;t)=\sum_{n=0}^{\infty}\frac{q_{n}(t)}{4^{n}}\mathopen{}\mathclose{{}\left(\frac{2t}{|x|}}\right)^{n/2}I_{n}\big{(}2\sqrt{2t|x|}\big{)}\mathrm{e}^{-2|x|}\mathrm{e}^{-3t/2}$
$\displaystyle=I_{0}(2\sqrt{2t|x|})\mathrm{e}^{-2|x|}\mathrm{e}^{-3t/2}$ (35)
$\displaystyle+\frac{1}{4}(1+t/2)\sqrt{\frac{2t}{|x|}}I_{1}\big{(}2\sqrt{2t|x|}\big{)}\mathrm{e}^{-2|x|}\mathrm{e}^{-3t/2}+\ldots$

Using the tail asymptotic $I_{n}(z)\sim\mathrm{e}^{z}/\sqrt{2\pi z}$ we get
the leading behaviour for large $|x|$

$p_{X}(x;t)\sim\big{(}4\pi\sqrt{2t|x|}\big{)}^{-1/2}\mathrm{e}^{2\sqrt{2t|x|}}\mathrm{e}^{-2|x|}\mathrm{e}^{-3t/2}$
(36)

confirming the exponential decay of the PDF visible from the simulations, see
Fig. 11 for a comparison.

Figure 11: Comparison of the simulated compound Poisson process, expansion
(VI) and asymptotical formula (36) for $t=2$. For smaller times $t\lesssim 1$
the fit is even better. For larger $t$, as Gaussian regime starts the
exponential behaviour is pushed towards larger $|x|$.

Equation (36) is similar to the large deviation property. We have shown that
$\ln p_{X}(x;t)\sim-tf(|x|/t)$ with the rate function
$f(y)=3/2+2y-2\sqrt{2y}$. This works for large $|x|$ and any $t$, but for the
larger $t$ the convergence becomes progressively slower; this is caused by the
coefficients $q_{n}(t)$ increasing with $t$. Nonetheless, it is not always
necessarily the best formulation of the problem. For small $t$ the rate
function $f(|x|/t)$ diverges, but we can still use Eq. (VI) because it is
exact. For example, if $|x|\to\infty$ and $t|x|\to\infty$ the result is the
same, but if we go with $t\to 0$ fast enough such that $t|x|\to 0$ we get $\ln
p_{X}(x;t)\sim 2t|x|-2|x|-3t/2$.

In any case, the most important factor is the exponential decay
$\exp(-2|x|)\exp(-3t/2)$ which dominates the shape of the PDF in the semi-log
scale. The rest are corrections which mostly move the straight line
$-2|x|-3t/2$ around the coordinate frame. This non-Gaussian Laplace behaviour
is observed in the “interim” time range in which the particle explores the
neighbourhood of several domains and the effective diffusivity has already
been reduced close to its long-time limit; it is the middle part of the range
shown in Figs 6 and 7.

##  VII Discussion

A robust and growing class of observations is exhibiting linear or power-law
mean-square displacement joined by the pronouncedly non-Gaussian probability
density, which most commonly has the exponential (Laplace) shape Wang  _et
al._ (2012); Metzler (2017); Bhattacharya  _et al._ (2013); Hapca  _et al._
(2009); Lampo  _et al._ (2017); Wang  _et al._ (2009) (see also the list of
references in Chechkin  _et al._ (2017)). Explaining this phenomena related to
Brownian, yet non-Gaussian diffusion is crucial for understanding transport
within these media and their biochemical properties, a timely issue especially
in biology and medicine. In spite of this demand, most of the studied
analytical models are very case-specific Kurtuldu  _et al._ (2011); Menzel and
Goldenfeld (2011) or do not provide description which would reflect the real
structure of the considered material, the most prominent examples being the
superstatistical Beck and Cohen (2003); Beck (2006, 2011); Ślęzak  _et al._
(2018) and the diffusing diffusivity Chubynsky and Slater (2014); Lanoiselée
_et al._ (2018); Chechkin  _et al._ (2017); Sposini  _et al._ (2020)
approaches. The line of research presented here follows the widespread
conviction that the physical origin of these observations lays in the
heterogeneity of the studied media and try to emulate it by randomizing the
parameters which appear in the dynamical equations of the diffusion.

In this work we solve a microscopical model of the heterogeneous medium which
is compartmentalised by the thick and thin barriers that strongly impede
diffusion and cause the transient-confinement of the diffusing particles. This
choice of model is motivated by the ubiquity of measurements showing this so-
called hop diffusion observed in biology Fujiwara  _et al._ (2016); Sadegh
_et al._ (2017); Chein  _et al._ (2019); Kusumi  _et al._ (2005); Krapf
(2018); Fujiwara  _et al._ (2002); Murase  _et al._ (2004) and the cage effect
in chemistry Weeks and Weitz (2002); Chaudhuri  _et al._ (2007); Weeks  _et
al._ (2000). Non-Gaussian diffusion was observed together with these phenomena
in multiple works, see the review Höfling and Franosch (2013). In many other
experiments the researchers were restricting their attention to only one of
these two aspects of transport, which leaves an open and intriguing
possibility of this link being much more prominent.

In our model the exponential tails of the probability density essentially stem
from the lack of spatial correlations in the random medium. For media with
more ordered microscopic structure, in particular with the periodically placed
barriers, the observed probability density is expected to be much more
Gaussian Wang  _et al._ (2019) (see Tanner (1978) for theoretical discussion).
However, if the barriers are placed independently, the distances between them
are exponential. Any given temporarily trapped particle then has uniform
distribution over its current domain, but the whole ensemble as a mixture of
those exhibits exponential tails.

The simplicity of this argument makes it quite general, but it is not enough
to specify the dynamics. For this we show how the local behaviour of the
particle is equivalent to a classical heat transfer problem and then derive
the escape times from the domains.

With the known domain sizes and transitions between them it becomes
straightforward to interpret the “hops” of the hop diffusion as the jumps of a
random walk. Monte Carlo simulations then show that this approximation is
indeed very close to the original diffusion and both are the cases of non-
Gaussian diffusion which at long times reaches Brownian limit.

This process is Fickian in both short and long timescales, but they differ by
the diffusivity, which transitions from the initial “unhindered” short-time
$D$ to the lower $D_{\text{eff}}=\expectationvalue{\kappa_{B}}D$; Such decay
of diffusivity was reported in numerous experimental works Murase  _et al._
(2004). The mean Biot number $\expectationvalue{\kappa_{B}}$ is determined by
distances between barriers and their permeability,
$\expectationvalue{\kappa_{B}}=\kappa\expectationvalue{L}$ and our results
extends the analytical proofs for the systems with periodically placed
barriers Tanner (1978); Powles  _et al._ (1992). This link between micro- and
macroscopic properties of the system may be verified experimentally by
comparing the estimated diffusivity to the measured parameters of the barriers
and their distribution in space.

It also causes the mean square displacement to be a convex function and the
diffusion becomes subdiffusive in the transient time range of the diffusivity
decrease. This is a well known phenomenon Höfling and Franosch (2013), but our
model also suggest that this behaviour is linked to the non-Gaussian
probability density with exponential tails which may be observed in the same
time regime. Subdiffusive non-Gaussian observations were reported in Lampo
_et al._ (2017); Weeks and Weitz (2002); Chaudhuri  _et al._ (2007); Weeks
_et al._ (2000). The derived leading behaviour of the probability density
(36), formulas for kurtosis (V) and logarithm of the characteristic function
(22) may be of use in analysing this type of experimental data.

The model which we consider is one-dimensional so one must be careful with
applying it to systems with higher dimensionality. Nevertheless, it provides
strong insights for the crucial aspects of the dynamics applicable to the less
idealised cases: the random walk behaviour should be observed for more general
systems divided into the localised, strongly confining domains. What changes
are the domain sizes and escape times. But, during our analytical calculations
of the non-Gaussianity measures and the probability density we introduced few
subsequent simplifications to the model, removing the dependence of the
environment, correlation and even approximating transition times. Yet, all the
obtained random walks were remarkably similar, showing a surprising level of
universality for our results.

It should be stressed that by this argument we actually establish not one, but
three connections, each remarkable on its own: between systems with barriers
and random walks, between random walks and Brownian, yet non Gaussian
diffusion, and finally, between non Gaussian diffusion and systems with
barriers. This approach opens new research possibilities in each of these
areas, but also offers immediate scientific benefits: we manage to express the
probability density of the displacements, non-Gaussianity measures and the
effective diffusion coefficients using three parameters which describe the
medium at microscopic level: in-domain diffusion coefficient, the average
distance between barriers and their permeability. These relations form a
valuable link between the molecular structure of the systems in question and
the measurements of the non-Gaussian diffusion.

## Appendix: Notation

  * •

i.i.d. independent and identically distributed

  * •

$\overset{d}{=}$ equal in distribution, has a distribution

  * •

$\mathcal{U}nif(a,b)$ uniform distribution,
$p_{\mathcal{U}nif(a,b)}(x)=\boldsymbol{1}_{(a,b)}(x)/(b-a)$

  * •

$\mathcal{G}am(\alpha,\lambda)$ gamma distribution,
$p_{\mathcal{G}am(\alpha,\lambda)}(x)=\frac{\lambda^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}\mathrm{e}^{-\lambda
x}$; it has mean $\alpha/\lambda$

  * •

$\mathcal{E}xp(\lambda)$ exponential distribution,
$\mathcal{E}xp(\lambda)=\mathcal{G}(1,\lambda)$; it has mean $1/\lambda$ thus
$\lambda$ is decay rate.

  * •

$\mathcal{L}ap(\mu,s)$ Laplace distribution,
$p_{\mathcal{L}ap(\mu,s)}(x)=\frac{1}{2s}\mathrm{e}^{-|x-\mu|/s}$; it has
variance $2s^{2}$

  * •

$I_{n}(x),K_{n}(x)$ modified Bessel functions of the second kind

###### Acknowledgements.

This work was supported by the Pazy foundation grant 61139927.

## References

## References

  * Feynman (1964) R. P. Feynman, The Brownian movement, in  _Feynman Lectures of Physics_ (Addison-Wesley, Boston, 1964) Chap. 6, pp. 41–45.
  * Humphries  _et al._ (2010) N. E. Humphries, N. Queiroz, J. R. M. Dyer, N. G. Pade, M. K. Musyl, K. M. Schaefer, D. W. Fuller, J. M. Brunnschweiler, T. K. Doyle, J. D. R. Houghton, G. C. Hays, C. S. Jones, L. R. Noble, V. J. Wearmouth, E. J. Southall, and D. W. Sims, Environmental context explains Lévy and Brownian movement patterns of marine predators, Nature 465, 1066 (2010).
  * Merton (1971) R. C. Merton, Optimum consumption and portfolio rules in a continuous-time model, J. Econ. Theory 3, 373 (1971).
  * Einstein (1905) A. Einstein, Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen, Ann. der Physik 322, 549 (1905).
  * Billingsley (1999) P. Billingsley,  _Convergence of Probability Measures_ , 2nd ed. (Wiley, 1999).
  * Shlesinger  _et al._ (1999) M. F. Shlesinger, J. Klafter, and G. Zumofen, Above, below and beyond Brownian motion, Am. J. Phys. 67, 1253 (1999).
  * Metzler and Klafter (2000) R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339, 1 (2000).
  * Coffey  _et al._ (1996) W. T. Coffey, Y. P. Kalmykov, and J. T. Waldron,  _The Langevin Equation_ (Word Scientific, 1996).
  * Wang  _et al._ (2012) B. Wang, J. Kuo, S. C. Bae, and S. Granick, When Brownian diffusion is not Gaussian, Nat. Mat. 11, 481 (2012).
  * Metzler (2017) R. Metzler, Gaussianity fair: The riddle of anomalous yet non-Gaussian diffusion, Biophys. J. 112, 413 (2017).
  * Bhattacharya  _et al._ (2013) S. Bhattacharya, D. Sharma, S. Saurabh, S. De, A. Sain, A. Nandi, and A. Chowdhury, Plasticization of poly(vinylpyrrolidone) thin films under ambient humidity: Insight from single-molecule tracer diffusion dynamics, J. Phys. Chem. B 117, 7771 (2013).
  * Hapca  _et al._ (2009) S. Hapca, J. W. Crawford, and I. M. Young, Anomalous diffusion of heterogeneous populations characterized by normal diffusion at the individual level, J. Roy. Soc. Interface 6, 111 (2009).
  * Lampo  _et al._ (2017) T. J. Lampo, S. Stylianidou, M. P. Backlund, P. A. Wiggins, and A. J. Spakowitz, Cytoplasmic RNA-protein particles exhibit non-Gaussian subdiffusive behavior, Biophys. J. 112, 532 (2017).
  * Wang  _et al._ (2009) B. Wang, S. M. Anthony, S. C. Bae, and S. Granick, Anomalous yet Brownian, Proc. Natl. Acad. Sci. U.S.A. 106, 15160 (2009).
  * Beck and Cohen (2003) C. Beck and E. G. Cohen, Superstatistics, Physica A 322, 267 (2003).
  * Beck (2006) C. Beck, Superstatistical Brownian motion, Prog. Theor. Phys. Suppl. 162, 29 (2006).
  * Beck (2011) C. Beck, Generalized statistical mechanics for superstatistical systems, Philos. T. Roy. Soc. A 369, 453 (2011).
  * Ślęzak  _et al._ (2018) J. Ślęzak, R. Metzler, and M. Magdziarz, Superstatistical generalised Langevin equation: non-Gaussian viscoelastic anomalous diffusion, New J. Phys. 20, 023026 (2018).
  * Chubynsky and Slater (2014) M. V. Chubynsky and G. W. Slater, Diffusing diffusivity: A model for anomalous, yet Brownian, diffusion, Phys. Rev. Lett. 113, 098302 (2014).
  * Lanoiselée  _et al._ (2018) Y. Lanoiselée, N. Moutal, and D. S. Grebenkov, Diffusion-limited reactions in dynamic heterogeneous media, Nat. Commun. 9, 4398 (2018).
  * Chechkin  _et al._ (2017) A. V. Chechkin, F. Seno, R. Metzler, and I. M. Sokolov, Brownian yet non-Gaussian diffusion: From superstatistics to subordination of diffusing diffusivities, Phys. Rev. X 7, 021002 (2017).
  * Ediger (2000) M. D. Ediger, Spatially heterogeneous dynamics in supercooled liquids, Ann. Rev. Phys. Chem. 51, 99 (2000).
  * He  _et al._ (2016) W. He, H. Song, Y. Su, L. Geng, B. J. Ackerson, H. B. Peng, and P. Tong, Dynamic heterogeneity and non-Gaussian statistics for acetylcholine receptors on live cell membrane, Nat Commun. 7, 178 (2016).
  * Sokolov (2010a) I. M. Sokolov, Ito, Stratonovich, Hänggi and all the rest: the thermodynamics of interpretation ., Chem. Phys. 375 (2010a).
  * (25) E. B. Postnikov, A. Chechkin, and I. M. Sokolov, Brownian yet non-Gaussian diffusion in heterogeneous media: from superstatistics to homogenization, New J. Phys. 22, 063046.
  * Novikov  _et al._ (2011) D. S. Novikov, E. Fieremans, J. H. Jensen, and J. A. Helpern, Random walks with barriers, Nat. Phys. 7, 508 (2011).
  * Song  _et al._ (2000) Y.-Q. Song, S. Ryu, and P. N. Sen, Determining multiple length scales in rocks, Nature 406, 11701 (2000).
  * Mair  _et al._ (1999) R. W. Mair, G. P. Wong, D. Hoffmann, M. D. Hürlimann, S. Patz, L. M. Schwartz, and R. L. Walsworth, Probing porous media with gas diffusion NMR, Phys. Rev. Lett. 83, 3324 (1999).
  * Fujiwara  _et al._ (2016) T. K. Fujiwara, K. Iwasawa, Z. Kalay, T. A. Tsunoyama, Y. Watanabe, Y. M. Umemura, H. Murakoshi, K. G. N. Suzuki, Y. L. Nemoto, N. Morone, and A. Kusumi, Confined diffusion of transmembrane proteins and lipids induced by the same actin meshwork lining the plasma membrane, Mol. Biol. Cell 27, 1101 (2016).
  * Sadegh  _et al._ (2017) S. Sadegh, J. L. Higgins, P. C. Mannion, M. M. Tamkun, and D. Krapf, Plasma membrane is compartmentalized by a self-similar cortical actin meshwork, Phys. Rev. X 7, 011031 (2017).
  * Chein  _et al._ (2019) M. Chein, E. Perlson, and Y. Roichman, Flow arrest in the plasma membrane, Biophys. J. 117, 1 (2019).
  * Kusumi  _et al._ (2005) A. Kusumi, C. Nakada, K. Ritchie, K. Murase, K. Suzuki, H. Murakoshi, R. S. Kasai, J. Kondo, and T. Fujiwara, Paradigm shift of the plasma membrane concept from the two-dimensional continuum fluid to the partitioned fluid: High-speed single-molecule tracking of membrane molecules, Annu. Rev. Biophys. Biomol. Struct. 34, 351 (2005).
  * Syková and Nicholson (2008) E. Syková and C. Nicholson, Diffusion in brain extracellular space, Physiol. Rev. 88, 1277 (2008).
  * Cory and Garroway (1990) D. G. Cory and A. N. Garroway, Measurement of translational displacement probabilities by NMR: an indicator of compartmentation., Magn. Reson. Med. 14, 435 (1990).
  * Yablonskiy  _et al._ (2002) D. A. Yablonskiy, A. L. Sukstanskii, J. C. Leawoods, D. S. Gierada, G. L. Bretthorst, S. S. Lefrak, J. D. Cooper, and M. S. Conradi, Quantitative in vivo assessment of lung microstructure at the alveolar level with hyperpolarized 3He diffusion MRI, Proc. Natl. Acad. Sci. U.S.A. 99, 3111 (2002).
  * Krapf (2018) D. Krapf, Compartmentalization of the plasma membrane, Curr. Opin. Cell Biol. 53, 15 (2018).
  * Fujiwara  _et al._ (2002) T. Fujiwara, K. Ritchie, H. Murakoshi, K. Jacobson, and A. Kusumi, Phospholipids undergo hop diffusion in compartmentalized cell membrane, J. Cell Biol. 157, 1071 (2002).
  * Murase  _et al._ (2004) K. Murase, T. Fujiwara, U. Yasuhiro, K. Suzuki, R. Iino, H. Yamashita, M. Saito, H. Murakoshi, K. Ritchie, and A. Kusumi, Ultrafine membrane compartments for molecular diffusion as revealed by single molecule techniques, Biophys. J. 86, 4075 (2004).
  * Latour  _et al._ (1994) L. L. Latour, K. Svoboda, P. P. Mitra, and S. C. H., Time-dependent diffusion of water in a biological model system, Proc. Natl. Acad. Sci. U.S.A. 91, 1229 (1994).
  * Sen (2003) P. N. Sen, Time-dependent diffusion coefficient as a probe of the permeability of the pore wall, J. Chem. Phys. 119, 9871 (2003).
  * Wang  _et al._ (2019) D. Wang, H. Wu, L. Liu, J. Chen, and D. K. Schwartz, Diffusive escape of a nanoparticle from a porous cavity, Phys. Rev. Lett. 123, 118002 (2019).
  * Weeks and Weitz (2002) E. R. Weeks and D. Weitz, Subdiffusion and the cage effect studied near the colloidal glass transition, Chem. Phys 284, 361 (2002), Strange Kinetics.
  * Chaudhuri  _et al._ (2007) P. Chaudhuri, L. Berthier, and W. Kob, Universal nature of particle displacements close to glass and jamming transitions, Phys. Rev. Lett. 99, 060604 (2007).
  * Weeks  _et al._ (2000) E. R. Weeks, J. C. Crocker, A. C. Levitt, A. Schofield, and D. A. Weitz, Three-dimensional direct imaging of structural relaxation near the colloidal glass transition, Science 287, 627 (2000).
  * Ebeling and Sokolov (2005) W. Ebeling and I. M. Sokolov,  _Statistical Thermodynamics and Stochastic Theory of Nonequilibrium Systems_ (World Scientific, 2005).
  * Kallenberg (2017) O. Kallenberg, Poisson and related processes, in  _Random Measures, Theory and Applications._ (Springer, Cham, 2017) Chap. 3, pp. 70–108.
  * Cox (1962) D. R. Cox,  _Renewal Theory_ (Methuen & Co., 1962).
  * (48) DLMF, NIST Digital Library of Mathematical Functions, http://dlmf.nist.gov/, Release 1.0.21 of 2018-12-15, f. W. J. Olver, A. B. Olde Daalhuis, D. W. Lozier, B. I. Schneider, R. F. Boisvert, C. W. Clark, B. R. Miller and B. V. Saunders, eds.
  * Note (1) In this more realistic scenario the left end of the initial domain lays at $\Theta L$ and the right one at $(\Theta-1)L$; variable $\Theta$ has the uniform distribution $\mathcal{U}nif(0,1)$. Calculation analogical to (1) yields incomplete gamma PDF $p_{X}(x)=\Gamma(0,|x|)/(2|x|)$ which has tails $\sim\exp(-|x|)/(2|x|)$, they are even thicker than $\exp(-2|x|)$. On the other hand this distribution has logarithmic singularity at $x=0$ which make the motion more constrained at short distances.
  * Sokolov (2010b) I. M. Sokolov, Ito, Stratonovich, Hänggi and all the rest: The thermodynamics of interpretation, Chem. Phys. 375, 359 (2010b), Stochastic processes in Physics and Chemistry (in honor of Peter Hänggi).
  * Sánchez-Palencia (1980) E. Sánchez-Palencia,  _Non-Homogeneous Media and Vibration Theory_ , Lecture Notes in Physics, Vol. 127 (Springer, 1980).
  * Mandrekar and Pilipenko (2016) V. Mandrekar and A. Pilipenko, On a Brownian motion with a hard membrane, Stat. Probil. Lett. 113, 62 (2016).
  * Erban and Chapman (2007) R. Erban and S. J. Chapman, Reactive boundary conditions for stochastic simulations of reaction–diffusion processes, Phys. Biol. 4, 16 (2007).
  * Andrews (2009) S. S. Andrews, Accurate particle-based simulation of adsorption, desorption and partial transmission, Phys. Biol. 6, 046015 (2009).
  * Note (2) In the overall valuable work of Powles _et al._ Powles  _et al._ (1992) a system with the regularly placed barriers at $x_{k}\in\\{\ldots,-1,0,1,2,\ldots\\}$ is considered. They claim they obtained the exact PDF for this particular case but no proof is provided and the result may be doubted.
  * Lejay (2006) A. Lejay, On the constructions of the skew Brownian motion, Probab. Surv. 3, 413 (2006).
  * Lejay (2016) A. Lejay, The snapping out Brownian motion, Ann. Appl. Probab. 26, 1727 (2016).
  * Lejay (2018) A. Lejay, A Monte Carlo estimation of the mean residence time in cells surrounded by thin layers, Math. Comput. Simulat. 143, 65 (2018), 10th IMACS Seminar on Monte Carlo Methods.
  * Nagylaki (1976) T. Nagylaki, Clines with variable migration, Genetics 83, 867 (1976).
  * Barton (2008) N. H. Barton, The effect of a barrier to gene flow on patterns of geographic variation, Genet. Res. 90, 139–149 (2008).
  * Borodin (2017) A. N. Borodin, Brownian local time, in  _Stochastic Processes._ (Birkhäuser, 2017) pp. 359–438.
  * Barthelemy  _et al._ (2008) P. Barthelemy, J. Bertolotti, and D. A. Wiersma, A lévy flight for light, Nature 453, 495 (2008).
  * Burioni  _et al._ (2010) R. Burioni, L. Caniparoli, and A. Vezzani, Lévy walks and scaling in quenched disordered media, Phys. Rev. E 81, 060101 (2010).
  * Burioni  _et al._ (2014) R. Burioni, E. Ubaldi, and A. Vezzani, Superdiffusion and transport in two-dimensional systems with lévy-like quenched disorder, Phys. Rev. E 89, 022135 (2014).
  * Grebenkov (2006) D. S. Grebenkov, Partially reflected Brownian motion: A stochastic approach to transport phenomena, in  _Focus on Probability Theory_ (Nova Science Publishers, 2006) pp. 135–169.
  * Bouchaud and Georges (1990) J.-P. Bouchaud and A. Georges, Anomalous diffusion in disordered media: Statistical mechanisms, models and physical applications, Phys. Rep. 195, 127 (1990).
  * Magdziarz  _et al._ (2015) M. Magdziarz, H. Scheffler, P. Straka, and P. Żebrowski, Limit theorems and governing equations for Lévy walks, Stoch. Proc. Appl. 125, 4021 (2015).
  * Tanner (1978) J. E. Tanner, Transient diffusion in a system partitioned by permeable barriers. Application to NMR measurements with a pulsed field gradient, J. Chem. Phys. 69, 1748 (1978).
  * Powles  _et al._ (1992) J. G. Powles, M. J. D. Mallett, G. Rickayzen, and W. A. B. Evans, Exact analytic solutions for diffusion impeded by an infinite array of partially permeable barriers, Proc. Roy. Soc. A 436, 391 (1992).
  * Crick (1970) F. Crick, Diffusion in embryogenesis, Nature 225, 420 (1970).
  * Grebenkov  _et al._ (2014) D. S. Grebenkov, D. Van Nguyen, and J.-R. Li, Exploring diffusion across permeable barriers at high gradients. I. Narrow pulse approximation, J. Magn. Reson. 248, 153 (2014).
  * Höfling and Franosch (2013) F. Höfling and T. Franosch, Anomalous transport in the crowded world of biological cells, Rep. Prog. Phys. 76, 046602 (2013).
  * Zaburdaev  _et al._ (2015) V. Zaburdaev, S. Denisov, and J. Klafter, Lévy walks, Rev. Mod. Phys. 87, 483 (2015).
  * Jakeman and Pusey (1978) E. Jakeman and P. N. Pusey, Significance of $K$ distributions in scattering experiments, Phys. Rev. Lett. 40, 546 (1978).
  * Kotz  _et al._ (2001) S. Kotz, T. J. Kozubowski, and K. Podgórski,  _The Laplace Distribution and Generalizations_ (Springer, 2001).
  * Papapantoleon (2008) A. Papapantoleon, An introduction to Lévy processes with applications in finance, arXiv:0804.0482 (2008), Vienna University of Technology, lecture notes.
  * Ślęzak  _et al._ (2019) J. Ślęzak, R. Metzler, and M. Magdziarz, Codifference can detect ergodicity breaking and non-Gaussianity, New J. Phys. 21, 053008 (2019).
  * Note (3) We note that this route will bring stronger results than using the Fourier transform (21) directly, as it is unfortunately hard to inverse analytically. One can use the method of steepest descent to uncover the tail behaviour of $p_{X}$, but the result does not seem to be very practical, describing only a far away range not easily available in experiments.
  * Buchak and Sakhno (2017) K. Buchak and L. Sakhno, Compositions of Poisson and Gamma processes, Mod. Stoch.: Theory Appl. 4, 161 (2017).
  * Kurtuldu  _et al._ (2011) H. Kurtuldu, J. S. Guasto, K. A. Johnson, and J. P. Gollub, Enhancement of biomixing by swimming algal cells in two-dimensional films, Proc. Natl. Acad. Sci. U.S.A. 108, 10391 (2011).
  * Menzel and Goldenfeld (2011) A. M. Menzel and N. Goldenfeld, Effect of Coulombic friction on spatial displacement statistics, Phys. Rev. E 84, 011122 (2011).
  * Sposini  _et al._ (2020) V. Sposini, D. S. Grebenkov, R. Metzler, G. Oshanin, and F. Seno, Universal spectral features of different classes of random-diffusivity processes, New J. Phys. 22, 063056 (2020).