# Fractional Advection Diffusion Asymmetry Equation, derivation, solution and
application

Wanli Wang1,2 and Eli Barkai 2 1 Department of Applied Mathematics, Zhejiang
University of Technology, Hangzhou 310023, China
2 Department of Physics, Advanced Materials and Nanotechnology Institute, Bar
Ilan University, Ramat-Gan 52900, Israel
wanliiwang<EMAIL_ADDRESS>
###### Abstract

The non-Markovian continuous-time random walk model, featuring fat-tailed
waiting times and narrow distributed displacements with a non-zero mean, is a
well studied model for anomalous diffusion. Using an analytical approach, we
recently demonstrated how a fractional space advection diffusion asymmetry
equation, usually associated with Markovian Lévy flights, describes the
spreading of a packet of particles. Since we use Gaussian statistics for jump
lengths though fat-tailed distribution of waiting times, the appearance of
fractional space derivatives in the kinetic equation demands explanations
provided in this manuscript. As applications we analyse the spreading of
tracers in two dimensions, breakthrough curves investigated in the field of
contamination spreading in hydrology and first passage time statistics. We
present a subordination scheme valid for the case when the mean waiting time
is finite and the variance diverges, which is related to Lévy statistics for
the number of renewals in the process.

  *

##  1 Introduction

Continuous-time random walk [1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11] (CTRW) is a
stochastic jump process for a random walk that jumps instantaneously from one
site to another, following a sojourn period on a site. See also recent works
in [12, 13, 14, 15, 16, 17, 18, 19, 20, 21]. For example, it was used to
describe dispersive transport in time of flight experiment of charge carriers
in disordered systems [22]. The non-Markovian characteristic of CTRW becomes
apparent, especially when the disorder is introduced into the system, as
observed in the sense of ensemble average [10]. Fractional kinetic equations
[23, 4, 24, 25, 26, 27, 28, 29] are a popular tool used to describe various
anomalous phenomena [30, 31]. When the bias of the system is a constant, the
fractional time diffusion equation follows [32, 33]

$\frac{\partial}{\partial
t}\mathcal{P}(x,t)=~\\!{}_{0}\mathcal{D}_{t}^{1-\alpha}\left[\frac{\partial^{2}}{\partial
x^{2}}-\frac{\partial}{\partial x}\right]\mathcal{P}(x,t),$ (1)

describing anomalous dynamics with $0<\alpha<1$. Here we make the assumption
that both the generalized mobility and the generalized diffusion coefficient
are set to one. The time Riemann-Liouville operator
$~\\!{}_{0}\mathcal{D}_{t}^{1-\alpha}$ introduces a convolution integral
involving a power-law kernel that decays slowly over time, which is related to
memory effects of this nonequilibrium system. In addition,
$~\\!{}_{0}\mathcal{D}_{t}^{1-\alpha}$ is directly governed by heavy-tailed
power law waiting times distributions

$\phi(\tau)\propto\tau^{-\alpha-1},~{}~{}~{}~{}\tau\to\infty,$ (2)

in the context of the underlying picture, the CTRW model.

Due to heavy-tailed waiting times, characterized by an infinite mean, Eq. (1)
illustrates sub-diffusion behaviors. However, our focus is on the super-
diffusion scenario. Our recent findings demonstrated that the positional
distribution, $\mathcal{P}(x,t)$, of the mentioned CTRW model follows the
fractional advection diffusion asymmetry equation (FADAE)

$\frac{\partial}{\partial t}\mathcal{P}=D\frac{\partial^{2}}{\partial
x^{2}}\mathcal{P}-V\frac{\partial}{\partial
x}\mathcal{P}+S\\!~{}\frac{\partial^{\alpha}}{\partial(-x)^{\alpha}}\mathcal{P}$
(3)

for $1<\alpha<2$.

Figure 1: A random walker moves on a one-dimensional lattice with fat-tailed
displacements that describe statistics of Lévy flights [see subplot (a)].
While, subplot (b) displays a trajectory of a CTRW walker with thin tailed
displacements, e.g. Gaussian, but the waiting time probability density
functions (PDFs) are fat-tailed. Under certain conditions, the packets of
spreading particles can appear the same, despite the significant differences
of the underlying paths.

Here $D$, $V$, and $S$ are transport coefficients to be discussed, and the
operator $\frac{\partial^{\alpha}}{\partial(-x)^{\alpha}}$ is the right
Riemann-Liouville derivative [34, 35] in regard to space. Thus the fat tailed
distribution of sojourn times Eq. (2), which leads to the fractional time
Fokker-Planck equation (1), when $0<\alpha<1$, can lead as shown here to the
fractional space differential equation (3) if $1<\alpha<2$. There is a
transition in the basic kinetic description of the transport. When the mean
waiting time diverges, we get the mentioned fractional time Fokker-Planck
equation. However, if the mean sojourn time is finite but the variance
diverges, we obtain a fractional space description. Our first goal is to
clarify this point.

Fractional space transport equations [34, 36, 37, 4, 38, 30, 39, 40] are
usually associated with Lévy flights [35, 41, 42, 43, 44, 45, 46, 47, 48]. See
subplot (a) in Fig. 1. The term Lévy flights, also referred to as Lévy motion,
was coined by Benoit Mandelbrot in honour of the French mathematician Paul
Lévy. The individual jump has a length that is distributed with the PDF
decaying at large positions as a fat-tailed jump length distribution with
diverging variance of the size of the jumps. Importantly, the Lévy flight
process is a Markovian process, hence the appearance of the fractional
framework in the context of a model with fat tailed distributions of waiting
times is of interest. In particular, any fitting of data using the fractional
space kinetic equation, cannot be used as strong evidence for Lévy flights, or
for a basic Markovian underlying process.

Subordination methods, based on an integral transformation, present a way of
solving fractional kinetic equations [3, 49, 50, 51, 52, 53, 54, 55, 56, 57].
To be more exact, when the mean waiting time diverges [33], subordination
methods map a classical Fokker-Planck equation onto a fractional diffusion one
with the fractional time derivative. For this case, the term ”inverse Lévy
transform” is sometimes used due to the number of renewals that follows the
inverse one-sided Lévy distribution. This well-known framework works when
$0<\alpha<1$. We will propose a new subordination method to solve a wide range
of problems valid when the mean waiting time is finite, but the variance is
diverging, namely the case $1<\alpha<2$.

The paper is organized as follows. We start by introducing the CTRW model and
give the corresponding concepts in Sec. 2. We derive the corresponding FADAE
and give our explanations in Sec. 3. In Sec. 4, the applications and
extensions of Eq. (3) are considered, ranging from the FADAE in two
dimensions, the positional distribution with the time-dependent bias and
variance, breakthrough curves, and the first passage time obtained from
subordination methods. Three limiting laws of FADAE, describing the large
asymmetric parameter $S$, the large diffusion term $D$, and the general
typical fluctuations, are discussed in Sec. 5. To conclude, we summarize our
findings in Sec. 6.

##  2 CTRW model

Now we define the CTRW model [4, 58, 10] discussed in this manuscript. Let us
consider the process starting at $t=0$ with the initial position $x=0$. A
walker is trapped on the origin for time $\tau_{1}$, then makes a jump and the
displacement is $\chi_{1}$; the walker is further trapped on $\chi_{1}$ for
time $\tau_{2}$, and then jumps to a new position; this process is then
repeated. Thus, the process is characterized by a set of waiting times
$\\{\tau_{1},\tau_{2},\cdots,\tau_{N},B_{t}\\}$ and the displacements
$\\{\chi_{1},\chi_{2},\cdots,\chi_{N}\\}$, where $B_{t}$ is the backward
recurrence time [59, 60] and the time dependent $N$ is the number of renewals
from the observation time zero to $t$. Specifically,
$\sum_{i=1}^{N}\tau_{i}+B_{t}=t$. These random variables, i.e.,
$\\{\tau_{1},\tau_{2},\cdots,\tau_{N}\\}$ and
$\\{\chi_{1},\chi_{2},\cdots,\chi_{N}\\}$ are mutually independent and
identically distributed (IID) random variables with common PDFs $\phi(\tau)$
and $f(\chi)$, respectively. Consider the broad distribution characterized by
a fat tail

$\displaystyle\phi(\tau)=\left\\{\begin{array}[]{ll}0,&\hbox{$\tau<\tau_{0}$;}\\\
\alpha\frac{\tau_{0}^{\alpha}}{\tau^{1+\alpha}},&\hbox{$\tau\geq\tau_{0}$,}\end{array}\right.$
(4)

where $\tau_{0}$ is a time scale and the power-law index $1<\alpha<2$. This
indicates that the mean
$\langle\tau\rangle=\int_{0}^{\infty}\tau\phi(\tau)d\tau=\alpha\tau_{0}/(\alpha-1)$
is finite, but not the variance. The Laplace transform will be used in solving
our problems. It is defined by
$\widehat{g}(s)=\int_{0}^{+\infty}\exp(-s\tau)g(\tau)d\tau$ for a well behaved
function $g(\tau)$. The Laplace transform of $\phi(\tau)$, $\tau\to s$, is [4,
59]

$\widehat{\phi}(s)\sim 1-\langle\tau\rangle s+b_{\alpha}s^{\alpha}$ (5)

with $b_{\alpha}=(\tau_{0})^{\alpha}|\Gamma(1-\alpha)|$. For $s\to 0$, we can
check that $\widehat{\phi}(0)=0$, which indicates that $\phi(\tau)$ in Eq. (4)
is normalized.

For the displacement PDF $f(\chi)$, we assume $\chi$ has a finite mean $a>0$
and the variance $\sigma^{2}$. In this manuscript, we consider a Gaussian
distribution

$f(\chi)=\frac{1}{\sqrt{2\sigma^{2}\pi}}\exp\left[-\frac{(\chi-a)^{2}}{2\sigma^{2}}\right].$
(6)

In Fourier space, $\widetilde{f}(k)=\exp(-\sigma^{2}k^{2}/2-ika)$ with
$\widetilde{f}(k)=\int_{-\infty}^{\infty}\exp(-ikx)f(x)dx$. Thus, the small
$k$ expansion of $\widetilde{f}(k)$ reads

$\widetilde{f}(k)\sim 1-ika-\frac{1}{2}(\sigma^{2}+a^{2})k^{2}.$ (7)

When the sojourn time PDF has a fat tail and the variance of the displacement
is finite, a wide range of anomalous behaviors emerge [3, 61]. As mentioned,
here we will focus on the widely applicable case $1<\alpha<2$ [62, 63, 64, 65,
66, 67, 68, 69].

##  3 Fractional-Space Advection Diffusion Asymmetry Equation

###  3.1 Calculation of the Positional distribution

We now consider the CTRW model and note that the PDF of the position of a
walker at time $t$ is

$P(x,t)_{\rm{CTRW}}=\sum_{N=0}^{\infty}Q_{t}(N)\chi_{N}(x)\to\int_{0}^{\infty}Q_{t}(N)\chi_{N}(x){\rm
d}N,$ (8)

where $Q_{t}(N)$ denotes the probability of the number of events during the
time interval $(0,t)$ and $\chi_{N}(x)$ is the probability that the walker is
on $x$ conditioned that it made $N$ steps. Equation 8 is known as the
subordination of the spatial process $x$ by the temporal process for $N$ [3,
33, 70, 54]. The technique was used in computer simulation in the context of
fractional Fokker-Planck dynamics [71], random diffusivity [56], population
heterogeneity [72] and one-dimensional Brownian search problem [55]. As
mentioned before, time intervals and displacements of walkers under
investigation are IID with common PDFs $\phi(\tau)$ and $f(\chi)$,
respectively. Note the relation between $\chi_{N-1}(x)$ and $\chi_{N}(x)$,
i.e., $\chi_{N}(x)=\int_{-\infty}^{\infty}\chi_{N-1}(y)f(x-y){\rm d}y$, and
the density of the particle’s displacement after the first step is given
exactly by $f(x)$, i.e., $\chi_{1}(x)=f(x)$. Then using the convolution of the
Fourier transform, it follows that as well known

$\chi_{N}(x)=\frac{1}{\sqrt{2\pi\sigma^{2}N}}\exp\left(-\frac{(x-aN)^{2}}{2\sigma^{2}N}\right).$
(9)

Based on the renewal process [59], in Laplace space $\widehat{Q}_{s}(N)$
follows $\widehat{Q}_{s}(N)=\widehat{\phi}(s)^{N}(1-\widehat{\phi}(s))/s$.
Considering the random variable $\epsilon=N-t/\langle\tau\rangle$ and taking
the inverse transform, then typical fluctuations follow [73, 59, 60]

$Q_{t}(N)\sim\frac{1}{(t/\overline{t})^{1/\alpha}}\mathcal{L}_{\alpha}\left(\frac{N-t/\langle\tau\rangle}{(t/\overline{t})^{1/\alpha}}\right)$
(10)

with

$\overline{t}=\langle\tau\rangle^{1+\alpha}/((\tau_{0})^{\alpha}|\Gamma(1-\alpha)|).$
(11)

Eq. (10) is valid in the limit of
$N-t/\langle\tau\rangle\propto(t/\overline{t})^{1/\alpha}$. Rare fluctuations
describing the scaling of $N-t/\langle\tau\rangle\propto t/\langle\tau\rangle$
were investigated in Ref. [60]. Here $\mathcal{L}_{\alpha}(x)$ is the Lévy
distribution [34, 35], defined by

$\mathcal{L}_{\alpha}(x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\exp(ikx)\exp[(-ik)^{\alpha}]{\rm
d}k.$ (12)

Utilizing Eqs. (8), (9) and (10), we have

$P(x,t)_{\rm{CTRW}}\sim\frac{1}{(t/\overline{t})^{1/\alpha}}\int_{0}^{\infty}\mathcal{L}_{\alpha}\left(\frac{N-t/\langle\tau\rangle}{(t/\overline{t})^{1/\alpha}}\right)\frac{\exp\Big{(}-\frac{(x-aN)^{2}}{2\sigma^{2}N}\Big{)}}{\sqrt[]{2\pi\sigma^{2}N}}{\rm
d}N.$ (13)

Note that Eq. (13) can be extended to more general situations, e.g. if the
jump size distribution is fat-tailed. Starting from Eq. (13) and changing
variable $y=(N-t/\langle\tau\rangle)/(t/\overline{t})^{1/\alpha}$, we get

$P(x,t)_{\rm{CTRW}}\sim\int_{-\frac{t}{\langle\tau\rangle}(\overline{t}/t)^{1/\alpha}}^{\infty}\mathcal{L}_{\alpha}(y)\frac{\exp\Big{(}-\frac{(x-a\frac{t}{\langle\tau\rangle}-ay(\frac{t}{\overline{t}})^{1/\alpha})^{2}}{2\sigma^{2}(t/\langle\tau\rangle+y(t/\overline{t})^{1/\alpha})}\Big{)}}{\sqrt[]{2\sigma^{2}\pi(\frac{t}{\langle\tau\rangle}+y(\frac{t}{\overline{t}})^{1/\alpha})}}{\rm
d}y,$ (14)

describing the scaling when $x-at/\langle\tau\rangle\propto
a(t/\bar{t})^{1/\alpha}$. Note that in the long time limit, the lower limit of
the integral in Eq. (14) can be replaced with $-\infty$; see Eq. (21) below.
The reason is as follows the dimensionless parameter

$\frac{t}{\langle\tau\rangle}\left({\overline{t}\over
t}\right)^{1/\alpha}=\left({t\over\tau_{0}}\right)^{1-1/\alpha}\left({\alpha\over\Gamma(2-\alpha)}\right)^{1/\alpha}\to\infty$
(15)

for $1<\alpha<2$ when $t\to\infty$. The scaling form of the position reads

$P_{\rm{CTRW}}(\xi,t)\sim\int_{-\frac{t}{\langle\tau\rangle}(\frac{\overline{t}}{t})^{\frac{1}{\alpha}}}^{\infty}\mathcal{L}_{\alpha}(y)\frac{\exp(-\frac{(\xi-y)^{2}}{2gg(y,t)})}{\sqrt[]{2\pi
gg(y,t)}}{\rm d}y$ (16)

with

$\xi=(x-at/\langle\tau\rangle)/(a(t/\overline{t})^{1/\alpha})$ (17)

and

$gg(y,t)=\frac{\sigma^{2}(t/\langle\tau\rangle+y(t/\overline{t})^{1/\alpha})}{a^{2}(t/\overline{t})^{2/\alpha}}.$
(18)

In the long time limit, Eq. (16) reduces to the well-known result [73, 74]

$P_{\rm{CTRW}}(\xi)\sim\mathcal{L}_{\alpha}(\xi),$ (19)

where we used the fact that $\exp[-z^{2}/(2\rho^{2})]/\sqrt{2\pi\rho^{2}}$
tends to $\delta(z)$ when $\rho\to 0^{+}$ from Eq. (19). More precisely, Eq.
(19) is valid for the fixed $\sigma$, $a$, and the long observation time $t$.
At the same time, we will consider a second dimensionless parameter
$\overline{\sigma}^{2}$,

$gg(y,t)\sim\frac{\sigma^{2}}{a^{2}}\left(\frac{t}{\tau_{0}}\right)^{1-\frac{2}{\alpha}}\left(\frac{\alpha}{\alpha-1}\right)^{\frac{2+\beta}{\beta}}|\Gamma(1-\alpha)|^{-\frac{2}{\alpha}}=\overline{\sigma}^{2},$
(20)

which quantifies the convergence to the Lévy distribution. We will then
discuss three limits. If $\overline{\sigma}^{2}\to 0$, Eq. (16) gives Eq.
(19). However, assume a situation when the bias is weak meaning
$\sigma^{2}/a^{2}$ is very large, this is important in many applications that
we consider, and corresponds to the linear response regime when the external
driving force determined by $a$ is small. Then the last limit will be to
consider a long time but $\bar{\sigma}^{2}$ finite. In this case, Eq. (14) is
valid. Refer to Appendix A for more detailed discussions on
$\overline{\sigma}^{2}$.

Figure 2: PDF of $\xi$ for various $\sigma$ listed in figure. Here the first
moment of waiting times, $\langle\tau\rangle=0.3$, implies that
$t/\langle\tau\rangle\approx 3333.3\gg 1$. The simulation results, represented
by the blue symbols, were obtained by averaging $5\times 10^{5}$ particles
with $\alpha=1.5$ and $t=1000$. The corresponding theoretical prediction,
shown by the solid red lines, was derived from Eq. (21), i.e., the exact
solution of Eq. (3), and it exhibits a perfect match with no fitting. Notice
that what the limit theorem Eq. (19) predicts here is the top dashed black
line, while our results valid for a variety of $a$ and $\sigma$ show the
deviations from the Lévy stable distribution $\mathcal{L}_{\alpha}(\xi)$ when
$\sigma=2/\sqrt{2},4/\sqrt{2},6/\sqrt{2}$, and $16/\sqrt{2}$. See further
discussion in Appendix A.

###  3.2 Derivation of Fractional-Space Advection Diffusion asymmetry Equation

Rewriting Eq. (14) with a lower limit $-\infty$, in the long time limit we get

$\mathcal{P}(x,t)=\int_{-\infty}^{\infty}\mathcal{L}_{\alpha}(y)\frac{\exp\left(-\frac{(x-a\frac{t}{\langle\tau\rangle}-ay(\frac{t}{\overline{t}})^{\frac{1}{\alpha}})^{2}}{2\sigma^{2}t/\langle\tau\rangle}\right)}{\sqrt[]{2\sigma^{2}\pi
t/\langle\tau\rangle}}{\rm d}y;$ (21)

see Fig. 2. Here we use the expression $\mathcal{P}(x,t)$ to denote the
positional distribution $P_{{\rm CTRW}}(x,t)$ in the long time limit. The
statistic of $\mathcal{P}(x,t)$ at a short time $t$ is now discussed. For this
case, Eq. (21) reduces to

$\mathcal{P}(x,t)\simeq\frac{\exp(-\frac{x^{2}}{2\sigma^{2}t/\langle\tau\rangle})}{\sqrt{2\pi\sigma^{2}t/\langle\tau\rangle}},$
(22)

where we ignored the terms $at/\langle\tau\rangle$ and
$ay(t/\overline{t})^{1/\alpha}$ in Eq. (21). It implies that for small $t$ the
Gaussian distribution is found. Besides, using $\delta(x)=\lim_{\rho\to
0^{+}}\exp(-x^{2}/(2\rho^{2}))/\sqrt{2\pi\rho^{2}}$, we can see that the
initial position of the process is $\mathcal{P}(x,t\to 0)=\delta(x)$.
Let us proceed with the discussion of the FADAE. Taking the Fourier transform
of Eq. (21) with respect to $x$ gives

$\displaystyle\begin{array}[]{ll}\mathcal{\widetilde{P}}(k,t)&=\displaystyle\int_{-\infty}^{\infty}\mathcal{L}_{\alpha}(y)\exp\left(-i\left(a\frac{t}{\langle\tau\rangle}+ay\left({t\over\overline{t}}\right)^{1/\alpha}\right)k-k^{2}\frac{\sigma^{2}t}{2\langle\tau\rangle}\right)dy\\\
&=\exp\left(-\frac{t\sigma^{2}}{2\langle\tau\rangle}k^{2}-iak\frac{t}{\langle\tau\rangle}+(-ik)^{\alpha}a^{\alpha}{t\over\overline{t}}\right).\end{array}$
(23)

Let $t$ be zero, Eq. (23) results in $\mathcal{\widetilde{P}}(k,t\to 0)=1$,
illustrating the fact that $\mathcal{P}(x,t\to 0)=\delta(x)$. By application
of the differential operator $\partial/\partial t$, we find FADAE in Fourier
space

$\frac{\partial}{\partial
t}\mathcal{\widetilde{P}}(k,t)=-\frac{\sigma^{2}}{2\langle\tau\rangle}k^{2}\mathcal{\widetilde{P}}(k,t)-ik\frac{a}{\langle\tau\rangle}\mathcal{\widetilde{P}}(k,t)+(-ik)^{\alpha}\frac{a^{\alpha}}{\overline{t}}\mathcal{\widetilde{P}}(k,t).$
(24)

Taking the inverse Fourier transform, in $(x,t)$ space, we get Eq. (3) given
in the introduction [53]. The transport constants of Eq. (3) are as follows

$D=\frac{\sigma^{2}}{2\langle\tau\rangle},~{}~{}~{}~{}V=\frac{a}{\langle\tau\rangle},~{}~{}~{}~{}S=\frac{a^{\alpha}}{\overline{t}}.$
(25)

Here $D$ carrying the dimension $[\sigma^{2}/\langle\tau\rangle]={\rm
cm}^{2}{\rm s}^{-1}$ is the diffusion constant. The drift of the process is
$V$. The last constant $S$ on the right-hand side of Eq. (3) describes the
asymmetric property of the process. The operator
$\partial^{\alpha}/\partial(-x)^{\alpha}$ is the right Riemann-Liouville
derivative [34, 35] with respect to space whose Fourier transform is
$\mathcal{F}[\frac{\rm{d}^{\alpha}}{\rm{d}(-x)^{\alpha}}g(x)]=(-ik)^{\alpha}\widetilde{g}(k)$.
In real space

$\frac{{\rm d}^{\alpha}g(x)}{{\rm
d}(-x)^{\alpha}}=\frac{(-1)^{n}}{\Gamma(n-\alpha)}\frac{{\rm d}^{n}}{{\rm
d}x^{n}}\int_{x}^{\infty}(y-x)^{n-\alpha-1}g(y)\rm{d}y,$ (26)

where $n$ is the smallest integer larger than $\alpha$. We have used $a>0$ in
Eq. 3, but as discussed below, when $a<0$, the left Riemann-Liouville
derivative comes into being. We note by passing that mathematically the left
Riemann-Liouville derivative is defined as [34, 35]

$\frac{{\rm d}^{\alpha}}{{\rm
d}x^{\alpha}}g(x)=\frac{1}{\Gamma(n-\alpha)}\frac{{\rm d}^{n}}{{\rm
d}x^{n}}\int_{-\infty}^{x}(x-y)^{n-\alpha-1}g(y)\rm{d}y$ (27)

with
$\mathcal{F}[\frac{\rm{d}^{\alpha}}{\rm{d}x^{\alpha}}g(x)]=(ik)^{\alpha}\widetilde{g}(k)$.

Here we present a summary of key insights extracted from Eq. (3):

  1. 1.

The first two constants, i.e., $D$ and $V$, are standard relations describing
classical advection-diffusion equation. If $a\to 0$, both the asymmetry
parameter $S$ and the fractional space derivative become negligible, see Eq.
25. The same holds for the drift $V$.

  2. 2.

Note that the solution $\mathcal{P}(x,t)$ of Eq. (3) is normalized for any
observation time $t$ since $\mathcal{\widetilde{P}}(k=0,t)=1$ according to Eq.
(23).

  3. 3.

We can check that the exact solution of Eq. (3) is Eq. (21). To be more
exactly, using Eq. (23), we find that the solution is the convolution of Lévy
distribution and Gaussian distribution, namely

$\mathcal{P}(x,t)=\frac{\mathcal{L}_{\alpha}\left(\frac{x}{(St)^{\alpha}}\right)}{(St)^{\alpha}}\bigotimes\frac{\exp\left(-{(x-Vt)^{2}\over
4Dt}\right)}{\sqrt{4\pi Dt}},$ (28)

where $\bigotimes$ is the convolution symbol of the Fourier transform with
respect to $x$.

  4. 4.

When particles are only allowed to move in one direction, for example,
$f(\chi)=\delta(\chi-a)$. For this case, $\chi_{N}(x)$ behaves as
$\chi_{N}(x)=\delta(x-aN)$. Utilizing Eqs. (8) and (10), we have

$P_{\rm{CTRW}}(x,t)\sim\frac{1}{a(t/\bar{t})^{1/\alpha}}\mathcal{L}_{\alpha}\left(\frac{\frac{x}{a}-\frac{t}{\langle\tau\rangle}}{(t/\bar{t})^{1/\alpha}}\right).$

This is the same as Eq. (19) after changing variables. Here $\chi_{N}(x)$ can
be treated as a delta function when the variance of displacement is finite.
The fractional equation related to this case is given by Eq. (3), where the
diffusion term is absent, indicating a strong bias scenario.

  5. 5.

Under the influence of the bias, the fractional derivative operator
$\frac{{\rm d}^{\alpha}}{{\rm d}(-x)^{\alpha}}$ is with respect to space,
though the distribution of the step length is Gaussian. In that sense, the
fat-tailed displacement is not a basic requirement for the fractional space
derivative operator.

  6. 6.

The order of taking limits of $t$ approaching infinity and $a$ going to zero
cannot be interchanged. As the bias of the system approaches zero, the
solution of Eq. (3) becomes a Gaussian distribution for any observation time
$t$, indicating normal diffusion. However, in the presence of any disturbance,
we observe super-diffusion in the long-time limit.

  7. 7.

According to Eq. (3) or (23), we get the asymptotic behavior of the mean of
the position $\langle x(t)\rangle\sim at/\langle\tau\rangle,$ growing linearly
with time $t$. As expected, Eq. (3) gives an infinite MSD, which is certainly
not possible. In Ref. [68] we find that the MSD is sensitive to rare events,
while Eq. (3) deals with typical fluctuations describing the central part of
the positional distribution. A detailed discussion about the mean squared
displacement (MSD) and fractional moments will be given in Appendix B.

  8. 8.

According to Eqs. (23) and (25), we have

$\mathcal{P}(k,t)=\exp\left[-Dk^{2}t-ikVt+S(-ik)^{\alpha}t\right].$ (29)

Instead of considering $\mathcal{P}(x,t)$, we investigate

$\mathcal{Q}(k,t)=\exp(ikVt)\mathcal{P}(k,t)=\exp\left[-Dk^{2}t+S(-ik)^{\alpha}t\right].$
(30)

Performing the Laplace transform with respect to $t$, we find

$\mathcal{Q}(k,s)=\mathcal{P}\left(k,s-iVk\right)=\frac{1}{s+\frac{\sigma^{2}}{2\langle\tau\rangle}k^{2}-\frac{a^{\alpha}}{\bar{t}}(-ik)^{\alpha}}.$
(31)

The last term of Eq. (31) demonstrates the competition between normal
diffusion and sup-diffusion. It can be seen that when
$\sigma^{2}k^{2}/(2\langle\tau\rangle)\gg|(-ik)^{\alpha}|a^{\alpha}/\bar{t}$,
i.e.,
$|k|\gg(a^{\alpha}2\langle\tau\rangle/(\bar{t}\sigma^{2}))^{1/(2-\alpha)}$, as
mentioned we get the Gaussian distribution. For the opposite limit, the Lévy
stable law is derived.

  9. 9.

The Lévy stable law and Gaussian distribution, which are the limit laws for
solving Eq. (3), are controlled by the dimensionless parameter
$\overline{\sigma}^{2}$. More details on this can be found in Sec. 5 in terms
of CTRW.

  10. 10.

We further investigate the situation when the bias is negative, i.e., $a<0$.
Substituting Eq. (25) into Eq. (29), we have

$\widetilde{\mathcal{P}}(k,t)=\exp\left(-\frac{\sigma^{2}}{2\langle\tau\rangle}k^{2}t-ik\frac{a}{\langle\tau\rangle}t+\frac{a^{\alpha}}{\bar{t}}(-ik)^{\alpha}t\right).$

Arranging the above equation, we find

$\widetilde{\mathcal{P}}(k,t)=\exp\left(-\frac{\sigma^{2}k^{2}t}{2\langle\tau\rangle}-ik\frac{at}{\langle\tau\rangle}+\frac{(-a)^{\alpha}}{\bar{t}}(ik)^{\alpha}t\right).$
(32)

Taking derivative with respect to $t$ and then performing the inverse Fourier
transform, Eq. (32) leads to

$\frac{\partial}{\partial
t}\mathcal{P}=\frac{\sigma^{2}}{2\langle\tau\rangle}\frac{\partial^{2}}{\partial
x^{2}}\mathcal{P}-\frac{a}{\langle\tau\rangle}\frac{\partial}{\partial
x}\mathcal{P}+\frac{(-a)^{\alpha}}{\bar{t}}\frac{\partial^{\alpha}}{\partial
x^{\alpha}}\mathcal{P}.$ (33)

In the case of a positive bias, we have the right Riemann-Liouville derivative
in Eq. (3). While, if the bias is negative, the fractional operator in Eq. (3)
is replaced by the left Riemann-Liouville derivative [see Eq. (27)],
controlling the right fat tail of the positional distribution.

  11. 11.

When waiting times have an infinite mean, i.e., Eq. (4) with $0<\alpha<1$, and
the displacement is still drawn from Gaussian distribution Eq. (6), the
corresponding fractional advection diffusion equation [32, 75, 4, 33, 30] is

$\frac{\partial}{\partial
t}\mathcal{P}(x,t)=~\\!{}_{0}\mathcal{D}_{t}^{1-\alpha}\Big{[}D\frac{\partial^{2}}{\partial
x^{2}}-V\frac{\partial}{\partial x}\Big{]}\mathcal{P}(x,t),$ (34)

where $D=\sigma^{2}/b_{\alpha}$ and $V=a/b_{\alpha}$ with
$b_{\alpha}=\tau_{0}^{\alpha}\Gamma(1-\alpha)$. In Eq. (34), the fractional
operator $\\!{}_{0}\mathcal{D}_{t}^{1-\alpha}$ is in reference to time rather
than the space. The fractional time operator
${}_{0}\mathcal{D}_{t}^{1-\alpha}$ called the Riemann-Liouville time
derivative is as follows [37, 76]

${}_{0}\mathcal{D}_{t}^{1-\alpha}g(t)=\frac{1}{\Gamma(\alpha)}\frac{\partial}{\partial
t}\int_{0}^{t}\frac{g(\tau)}{(t-\tau)^{1-\alpha}}d\tau$ (35)

with $0<\alpha<1$. Further, when $0<\alpha<1$, the positional distribution has
a right fat tail. While, in the context of $1<\alpha<2$, the positional
distribution has a left one. It can be seen that the forms of fractional
equations with different $\alpha$ show great differences.

  12. 12.

Subordination discussed here is vastly different from the case for
$0<\alpha<1$. When $1<\alpha<2$, the number of renewals follows the two-sided
Lévy stable law with an infinite variance, instead of the one-sided Lévy
distribution (also called the Mittag-Leffler distribution).

  13. 13.

Solutions we discussed are valid for free boundary conditions. One can use the
subordination scheme Eq. (13) to consider other cases. The key idea is that
the number of steps in the long time limit is distributed according to the
Lévy law given in Eq. (10). Then a normal process is considered, for example
diffusion in a finite medium, with the operational time $N$. The Lévy central
limit theorem is then used to transform the operational time $N$ (also the
number of steps) to the laboratory time $t$.

  14. 14.

Benson et al. proposed an equation called fractional advection dispersion
equation to describe contaminant source [44, 46]

$\frac{\partial}{\partial t}\mathcal{P}=-V\frac{\partial}{\partial
x}\mathcal{P}+K\left[p\frac{\partial^{\alpha}}{\partial
x^{\alpha}}\mathcal{P}+q\frac{\partial^{\alpha}}{\partial(-x)^{\alpha}}\mathcal{P}\right],$
(36)

which is obtained from Lévy flights; see Fig. 1 (a). If $D=0$ and $p=0$, Eq.
(3) is consistent with Eq. (36). In Ref. [77], the authors use $V=0.8$m/h,
$D=0$, $K=2.8{\rm m}^{1.51}/{\rm h}$, $q=1$, and $p=0$ to fit experimental
contaminant data in terms of the breakthrough curves. According to this, one
may interpret the data as coming from Lévy flights. On the other hand, here we
want to stress that these data may stem from CTRW model with a strong bias in
the long time limit, i.e., these experimental data are consistent with CTRW
with a fat-tailed waiting time distribution. In that sense, it indicates that
both Lévy flights and CTRW can be used as a tool to describe the contaminant
in disorder systems. This holds as long as we have information on the
spreading packets and no physical insight into the underlying trajectories.

##  4 Applications

As mentioned before, fractional diffusion equations are widely used as a tool
to describe the dispersion of contaminants and other phenomena observed in
non-equilibrium systems. For that, we explore several applications of FADAE
(3). These applications encompass calculating FADAE in two dimensions,
determining the positional distribution with the time-dependent bias and
variance, analyzing breakthrough curves [66, 78], estimating the first passage
time, and comparing them to CTRW dynamics.

###  4.1 FADAE in two dimensions

Here we consider the two dimensional FSDAE which is not only of theoretical
significance, but also potential application value [79, 80, 81, 66]. For
simplification, the bias of the system is only in $x$ direction and the
equation reads

$\frac{\partial}{\partial
t}\mathcal{P}(x,y,t)=D_{x}\frac{\partial^{2}}{\partial
x^{2}}\mathcal{P}(x,y,t)-V_{x}\frac{\partial}{\partial
x}\mathcal{P}(x,y,t)+S_{\vskip 3.0pt plus 1.0pt minus
1.0ptx}\frac{\partial^{\alpha}}{\partial(-x)^{\alpha}}\mathcal{P}(x,y,t)+D_{y}\frac{\partial^{2}}{\partial
y^{2}}\mathcal{P}(x,y,t).$ (37)

Here $D_{x}$, $V_{x}$, $S_{x}$, and $D_{y}$ are constants. Taking double
Fourier transforms, $x\to k_{x}$ and $y\to k_{y}$, we get

$\widetilde{\mathcal{P}}(k_{x},k_{y},t)=\exp\left(-D_{x}(k_{x})^{2}t-iV_{x}k_{x}t+(-ik_{x})^{\alpha}tS_{x}-D_{y}(k_{y})^{2}t\right).$
(38)

Using convolution properties of Fourier transform, Eq. (38) yields

$\mathcal{P}(x,y,t)=\int_{-\infty}^{\infty}\mathcal{L}_{\alpha}(z)\frac{\exp\left(-\frac{(x-V_{x}t-(tS_{x})^{1/\alpha}z)^{2}}{4D_{x}t}\right)}{\sqrt{4\pi
D_{x}t}}\frac{\exp\left(-\frac{y^{2}}{4D_{y}t}\right)}{\sqrt{4\pi D_{y}t}}dz,$
(39)

where $\mathcal{L}_{\alpha}(z)$ is the non-symmetric Lévy stable distribution
Eq. (12). Below we obtain four transport coefficients given in Eq. (37). In
the language of the CTRW model, the displacement follows

$f(x,y)=\frac{1}{\sqrt{2(\sigma_{x})^{2}\pi}}\exp\left(-\frac{(x-a_{x})^{2}}{2(\sigma_{x})^{2}}\right)\times\frac{1}{\sqrt{2(\sigma_{y})^{2}\pi}}\exp\left(-\frac{y^{2}}{2(\sigma_{y})^{2}}\right),$
(40)

where $a_{x}$, $\sigma_{x}\neq 0$, and $\sigma_{y}\neq 0$ are constants. In
double Fourier spaces, $x\to k_{x}$ and $y\to k_{y}$, we get

$\widetilde{f}(k_{x},k_{y})=\exp\left(-ik_{x}a_{x}-\frac{1}{2}(\sigma_{x})^{2}(k_{x})^{2}-\frac{1}{2}(\sigma_{y})^{2}(k_{y})^{2}\right).$
(41)

For the waiting time, we continue to utilize the fat-tailed power law
distribution Eq. (4). Using Eq. (25) and the same arguments as before, we find

$D_{x}=\frac{\sigma_{x}^{2}}{2\langle\tau\rangle},~{}~{}~{}V_{x}=\frac{a_{x}}{\langle\tau\rangle},~{}~{}~{}S_{x}=\frac{(a_{x})^{\alpha}}{\overline{t}},~{}~{}~{}D_{y}=\frac{\sigma_{y}^{2}}{2\langle\tau\rangle}.$
(42)

It can be seen that when the external force is only in the $x$ direction, the
fat tail of the packet of spreading particle $\mathcal{P}(x,y,t)$ is with
respect to $x$. Note that the exact solution of Eq. (37) is Eq. (39), i.e.,
$P_{\rm{CTRW}}(x,y,t)$ in the long time limit. What we want to stress is that
when the bias is also in the $y$ direction, the fractional space operator with
respect to $y$ should be added. See the solution of Eq. (37) in Fig. 3 and the
marginal distribution in Fig. 4.

Figure 3: The dynamic of $\mathcal{P}(x,y,t)/{\rm max}\mathcal{P}(x,y,t)$
obtained from Eq. (37) for the fractional advection-diffusion process in two
dimensions. However, in the inset, the results are projected. The parameters
of CTRW are $\tau_{0}=0.1$, $\alpha=1.5$, $a_{x}=2$, $t=1000$ and
$\sigma_{x}=\sigma_{y}=5$.
Figure 4: Plot of marginal distribution Eq. (37) compared with CTRW
simulations. Here we denote $\mathcal{P}_{2}(x,t)$ and $\mathcal{P}_{2}(y,t)$
as the marginal distribution of $\mathcal{P}(x,y,t)$ with respect to $x$ and
$y$, respectively. It can be seen that $\mathcal{P}_{2}(x,t)$ is the same as
the solution of Eq. (3) and $\mathcal{P}_{2}(y,t)$ is a Gaussian distribution.
The parameters are the same as in Fig. 4.

###  4.2 Propagator with the time-dependent bias and variance

Figure 5: Illustrations of the time-dependent bias, e.g., $a(t)$, when the
process begins at time $t=0$ and ends at time $t$. In our simulations, the
mean displacements of four states, i.e., bias, are determined by
$\\{a,4a\mu/(\mu+1),a/(2/3+\mu),a\\}$.

Motivated by [66], we consider the time-dependent bias and variance determined
by four states. In some sense, assumptions are more closely aligned with the
real-world spread of pollutants. This idea yields some interesting results,
for example complex structure of breakthrough curves; see below. We suppose
that the rapid injection of particles is done immediately after starting
observing the process. In other words, the initial positions of all particles
are $\mathcal{P}(x,0)=\delta(x)$. Here we simulate the particles consisting of
four states: (i) after the injection of the particles at time $t=0$, they
undergo CTRW processes with a constant bias determined by $a$ in time interval
$(0,t_{1})$ (here $a$ is the mean of displacements), (ii) during the time
interval $t_{1}<t<t_{2}$ we increase the force, namely replace the bias $a$
with $4a\mu/(\mu+1)$ for $\mu\geq 1/3$, (iii) we further decrease the bias to
$a/(2/3+\mu)$ from time $t_{2}$ to $t_{3}$, and (iv) finally we return to the
state (i) from time $t_{3}$; see Fig. 5 for the statistics of the bias versus
time. Here $\mu$ is a constant that controls the strength of the force. Let
$\sigma=5$ for states (i), (ii), and (iv), i.e.,
$\sigma_{1}=\sigma_{2}=\sigma_{4}=5$. While, for state (iii), we use
$\sigma_{3}=0.1$. In particular, when $\mu=1/3$, the bias of all states
mentioned above is the same but the variance is still time-dependent. Note
that the elapsed time of each state should be a bit long, otherwise one can
find that the difference between them is not very large. In Fourier space, the
initial condition satisfies $\widetilde{\mathcal{P}}(k,0)=1$. In view of the
special expression of $\widetilde{\mathcal{P}}(k,t)$,
$\widetilde{\mathcal{P}}(k,t)$ for different states is as follows

$\widetilde{\mathcal{P}}(k,t)=\left\\{\begin{array}[]{ll}\exp(-c_{11}k^{2}-ic_{12}k+c_{13}(-ik)^{\alpha}),&\hbox{$0<t\leq
t_{1}$;}\\\
\exp(-c_{21}k^{2}-ic_{22}k+c_{22}(-ik)^{\alpha}),&\hbox{$t_{1}<t\leq
t_{2}$;}\\\
\exp(-c_{31}k^{2}-ic_{32}k+c_{33}(-ik)^{\alpha}),&\hbox{$t_{2}<t\leq
t_{3}$;}\\\
\exp(-c_{41}k^{2}-ic_{42}k+c_{43}(-ik)^{\alpha}),&\hbox{$t_{3}<t$,}\end{array}\right.$
(43)

with

$c_{m1}=\left\\{\begin{array}[]{ll}t\frac{\sigma_{1}^{2}}{2\langle\tau\rangle},&\hbox{$m=1$;}\\\
t\frac{\sigma_{2}^{2}}{2\langle\tau\rangle}+\frac{t_{1}(\sigma_{1}^{2}-\sigma_{2}^{2})}{2\langle\tau\rangle},&\hbox{$m=2$;}\\\
t\frac{\sigma_{3}^{2}}{2\langle\tau\rangle}+\frac{t_{1}(\sigma_{1}^{2}-\sigma_{2}^{2})}{2\langle\tau\rangle}+\frac{t_{2}(\sigma_{2}^{2}-\sigma_{3}^{2})}{2\langle\tau\rangle},&\hbox{$m=3$;}\\\
t\frac{\sigma_{4}^{2}}{2\langle\tau\rangle}+\frac{t_{1}(\sigma_{1}^{2}-\sigma_{2}^{2})}{2\langle\tau\rangle}+\frac{t_{2}(\sigma_{2}^{2}-\sigma_{3}^{2})}{2\langle\tau\rangle}+\frac{t_{3}(\sigma_{3}^{2}-\sigma_{4}^{2})}{2\langle\tau\rangle},&\hbox{$m=4$,}\end{array}\right.$
(44)

$c_{m2}=\left\\{\begin{array}[]{ll}a_{1}\frac{t}{\langle\tau\rangle},&\hbox{$m=1$;}\\\
\left[a_{1}t_{1}+a_{2}(t-t_{1})\right]\frac{1}{\langle\tau\rangle},&\hbox{$m=2$;}\\\
\left[a_{1}t_{1}+a_{2}(t_{2}-t_{1})+a_{3}(t-t_{2})\right]\frac{1}{\langle\tau\rangle},&\hbox{$m=3$;}\\\
\left[a_{1}t_{1}+a_{2}(t_{2}-t_{1})+a_{3}(t_{3}-t_{2})+a_{4}(t-t_{3})\right]\frac{1}{\langle\tau\rangle},&\hbox{$m=4$,}\end{array}\right.$
(45)

and

$c_{m3}=\left\\{\begin{array}[]{ll}a_{1}^{\alpha}t\frac{1}{\overline{t}},&\hbox{$m=1$;}\\\
\left[a_{1}^{\alpha}t_{1}+a_{2}^{\alpha}(t-t_{1})\right]\frac{1}{\overline{t}},&\hbox{$m=2$;}\\\
\left[a_{1}^{\alpha}t+a_{2}^{\alpha}(t_{2}-t_{1})+a_{3}^{\alpha}(t-t_{2})\right]\frac{1}{\overline{t}},&\hbox{$m=3$;}\\\
\left[a_{1}^{\alpha}t_{1}+a_{2}^{\alpha}(t_{2}-t_{1})+a_{3}^{\alpha}(t_{3}-t_{2})+a_{4}^{\alpha}(t-t_{3})\right]\frac{1}{\overline{t}},&\hbox{$m=4$.}\end{array}\right.$
(46)

Here $m=1,2,3,4$ is related to the number of states. Recall that
$\bar{t}=\langle\tau\rangle^{1+\alpha}/[(\tau_{0})^{\alpha}|\Gamma(1-\alpha)|]$.
In calculations, the main idea is that the final position of each state is
treated as the initial position of the next stage. The inverse Fourier
transform of Eq. (43) yields

$\mathcal{P}(x,t)=\int_{-\infty}^{\infty}\frac{1}{\sqrt{4\pi
c_{m1}}}\exp\left[-\frac{(x-y-
c_{m2})^{2}}{4c_{m1}}\right]\frac{1}{(c_{m3})^{1/\alpha}}\mathcal{L}_{\alpha}\left[\frac{y}{(c_{m3})^{1/\alpha}}\right]dy,$
(47)

describing the positional distribution of the mentioned four states. Note that
when $0<t<t_{1}$, the solution Eq. (47) is the same as Eq. (21) or that of Eq.
(3). Below, we use Eq. (47) to predict breakthrough curves [82].

###  4.3 Breakthrough curves

Contamination spreading, as one of the most crucial problems ranging from
environment to agriculture, has attracted a lot of attention [80, 81, 66, 83,
84]. In real systems, how to quantity the contaminant transport is quite a
challenge for researchers since the processes are so complex. In [66], Nissan,
Dror and Berkowitz considered the spreading particle with changing conditions
to illustrate the positive and negative effects of the ambient environment. It
is no wonder that the changing force fields are quite common and vital in the
real world. Theoretical predictions about breakthrough curves, which can be
directly measured in experiments, were made using the CTRW particle tracking
approach in [66]. While, here we choose Eq. (47) as a tool to describe
breakthrough curves. The breakthrough curves are measured in the sense of the
distribution, indicating the probability of particles being at $x_{b}$ at a
specific time $t$. The differences and similarities between the constant and
the time-dependent force are quite interesting. When particles experience
small disturbances or the total time of the first three states is short, the
response of time-dependent force, calculated by the density of the position,
is weak except for a ‘shift’. However, when we increase the total time of the
first three states, the mentioned two cases exhibit great differences; see
Fig. 6. Note that in our discussions the detection point is $x_{b}=1600$.
Choosing the suitable detection point, i.e., $x_{b}$, is a critical factor in
observing the distinctive structure of the particle packet. For a short time
$t$, such as $t=10$, it can be seen that $\mathcal{P}(x_{b},t)\to 0$ since all
the particles are on the way to the site $x_{b}$ pushed by the force and need
more time to visit the position $x_{b}$. As particles enter the second state,
they are moving faster than the state (i) for $\mu=6$, resulting in a rapid
increase in the probability of the position. During the time interval
$(t_{2},t_{3})$, the disturbance and the variance are weak, leading to a more
or less flat breakthrough curve. At the end of this state, there are still
numerous particles being located at the left side of $x_{b}$. In the last
state, $\mathcal{P}(x_{b},t)$ slowly tends to zero.

Figure 6: Theoretical predictions of breakthrough curves together with the
simulations plotted by the symbols. The solid lines are the theories based on
Eq. (47). Time intervals of the four states are $100,100,100,700$,
respectively. Packets of particles are asymmetric in the long time limit,
which are the same for both constant and time-dependent cases. Semi-log scales
are represented to show asymmetric properties and the heavy tail; see the
inset. Breakthrough curves are observed at the position $x_{b}=1600$. Here the
parameters are $\mu=1/3$ (red line), $\mu=6$ (black line),
$\tau_{0}=0.1,\alpha=1.5$, $\sigma_{1}=\sigma_{2}=\sigma_{4}=5$, and
$\sigma_{3}=0.1$.

###  4.4 First passage time

So far we considered problems with free boundary conditions. The boundary
conditions for fractional space diffusion equations are generally non-trivial.
This is because fractional operators are non-local in space. Our approach to
the problem is to use the subordination method, namely rely on Eqs. (8,10).
The idea is as follows: We showed already that fluctuations of the number of
renewals $N$ is given by Lévy statistics. We may think about the whole problem
as a normal process, where $N$ is an operational time. Then this normal
process is transformed from $N$ to the laboratory time $t$. By normal process,
we mean, for example, Brownian diffusion with advection in a finite domain,
with reflecting or absorbing boundary conditions.

Here we consider as an example the case of the first passage time problem [85,
86, 87, 88, 89]. The first passage time serves as a significant tool for
describing the time required by a particle to reach a specific point or an
absorbing barrier. For standard diffusion this means that we use an absorbing
boundary condition (see below). We will later show how to use subordination to
obtain the solution for the anomalous process.

Before doing so we will briefly delve into the method of images. It was shown
previously, that for Lévy flights the image method fails [90, 91, 92, 93].
However, in our case, we do not have any Lévy flight as an underlying process.
Instead, we use fat-tailed distributions for waiting times. We will first show
that the method of images generally fails in our case, however when the bias
is weak it is a reasonable approximation. We will then turn to the more
general solution using subordination.

####  4.4.1 Recap for the image method and it’s failure

Let us determine the first passage probability based on Eq. (3) for a
diffusing particle that starts at $x_{0}=0$. As mentioned, we use the image
method and show that while it is generally not a good approximation, it works
well when the bias is small. Here we assume that the absorbing boundary
condition is $x_{ab}$ and use $C(x,t)$ to denote the concentration. Thus, when
$x\geq x_{ab}$, we have $C(x,t)=0$. The most appealing approach, dealing with
the absorbing boundary condition, is the method of images, which stems from
electrostatics. As mentioned we will soon show, following Eq. (3), that this
method does not work when the bias is not terribly small. In real space, the
image method leads to [33, 85],

$C(x,t)=\mathcal{P}(x,t)-W\mathcal{P}(x-2x_{ab},t),~{}-\infty<x\leq x_{ab},$
(48)

where the unknown parameter $W$ is defined below. Here $\mathcal{P}(x,t)$
denotes the positional distribution without absorbing condition, for example,
below we use the solution of Eq. (3) to check the validity of the image
method. It can be seen that the concentration is the difference between
$\mathcal{P}(x,t)$ and $\mathcal{P}(x-2x_{ab},t)$ with a weight $W$ determined
by $x_{ab}$. Recall that $C(x,t)$ vanish on $x_{ab}$, i.e., $C(x,t)=0$. Then
$W$ obeys

$\mathcal{P}(x_{ab},t)-W\mathcal{P}(-x_{ab},t)=0.$

Thus, we get with the image method

$W=\frac{\mathcal{P}(x_{ab},t)}{\mathcal{P}(-x_{ab},t)}.$ (49)

Consider the well-studied case of a particle undergoing normal diffusion (Eq.
(3) with $S=0$) with an absorbing condition at $x=x_{ab}$. Thus,
$\mathcal{P}(x,t)$ in the absence of the absorbing boundary condition follows
Gaussian distribution

$\mathcal{P}(x,t)=\frac{1}{\sqrt[]{4\pi
Dt}}\exp\Big{(}-\frac{(x-Vt)^{2}}{4Dt}\Big{)},$ (50)

which leads to

$W=\exp\left(x_{ab}\frac{V}{D}\right)$ (51)

according to Eq. (49). Note that in the general case $W$ is determined by
$x_{ab}$, $D$, $V$ and $S$.

Based on Eqs. (48) and (49), mathematically, the concentration is as follows

$C(x,t)=\mathcal{P}(x,t)-\frac{\mathcal{P}(x_{ab},t)}{\mathcal{P}(-x_{ab},t)}\mathcal{P}(x-2x_{ab},t)$
(52)

with $-\infty<x\leq x_{ab}$. Equation (52) is demonstrated in Fig. 14 in
Appendix C for different absorbing conditions. We see here that for weak
fields the image theory works since the process is almost Gaussian.

####  4.4.2 Subordinating the first passage problem

To solve this problem, we use a subordination technique that involves
utilizing the method of images on the $\chi_{N}(x)$. Based on Eq. (8), the
discrete form of $C(x,t)$ follows

$C(x,t)\sim\sum_{N=0}^{\infty}Q_{t}(N)\chi^{*}_{N}(x),$ (53)

describing the subordinated process $x$ as a function of time $N$. Here
$Q_{t}(N)$ is the PDF of the number of the renewals given in Eq. (10) and
$\chi^{*}_{N}(x)$ is the solution of the ordinary Fokker-Planck equation

$\frac{\partial}{\partial N}P(x,N)=a\frac{\partial}{\partial
x}P(x,N)-\frac{\delta^{2}}{2}\frac{\partial^{2}}{\partial x^{2}}P(x,N)$ (54)

with the initial conditions $P(x,N=0)=\delta(x)$ and different boundary
conditions, be specific, absorbing boundary condition at $x=x_{ab}$. The same
approach can be used for other types of boundary conditions. Setting the
absorbing condition $x=x_{ab}$ on Eq. (54) and using the method of images for
the solution of Eq. (54) yield

$\chi^{*}_{N}(x)=\frac{\exp\Big{(}-\frac{(x-aN)^{2}}{2\sigma^{2}N}\Big{)}}{\sqrt[]{2\pi\sigma^{2}N}}-\exp\left(\frac{2ax_{ab}}{\delta^{2}}\right)\frac{\exp\Big{(}-\frac{(x-2x_{ab}-aN)^{2}}{2\sigma^{2}N}\Big{)}}{\sqrt[]{2\pi\sigma^{2}N}}.$
(55)

In the long time limit, the continuous form of Eq. (53) follows

$C(x,t)=\left(\frac{\overline{t}}{t}\right)^{\frac{1}{\alpha}}\int_{0}^{\infty}\mathcal{L}_{\alpha}\left(\frac{N-t/\langle\tau\rangle}{(t/\overline{t})^{1/\alpha}}\right)\chi^{*}_{N}(x){\rm
d}N.$ (56)

Equation (56) is verified in Fig. 7 showing a perfect match. We now compare
solution Eq. (56) and that of Eq. (3) with free boundary conditions. For that,
the solution of Eq. (3) without absorbing condition was plotted by dashed
black lines. As expected, when $x$ is much smaller than $x_{ab}$, $C(x,t)$ is
consistent with Eq. (3), i.e., the random particles are not affected by the
absorbing condition, or most of particles have not yet reached the position
near $x_{ab}$. In addition, when the bias is strong, Eq. (3) agrees with Eq.
(56) for $x<x_{ab}$ at least to the naked eye. Roughly speaking, this is
because when the bias is strong, no particles are coming back. Thus, the
absorbing condition under study loses its role. While, if $x$ approaches to
$x_{ab}$ and the bias is weak, Eq. (3) and Eq. (56) illustrate two different
behaviors. See (d) in Fig. 7.

We further consider the first passage time using the survival probability. The
survival probability $\mathcal{S}(t)$ describing the probability that the
particles do not arrive on the position $x_{ab}$ until the time $t$ reads

$\mathcal{S}(t)=\int_{-\infty}^{x_{ab}}C(x,t)dx.$ (57)

Let $t_{f}$ be the time to visit the position $x_{ab}$ for the first time.
Utilizing Eq. (57), the PDF of $t_{f}$ reads

$\varphi(t_{f})=-\frac{d\mathcal{S}(t_{f})}{dt_{f}}.$ (58)

From Eqs. 53 and 58, we have

$\displaystyle\begin{array}[]{ll}\varphi(t_{f})&\sim-\sum_{N=1}^{\infty}\frac{d}{dt_{f}}Q_{t_{f}}(N)\int_{-\infty}^{x_{ab}}\chi_{N}^{*}(x)dx\\\
&=\frac{1}{2}\sum_{N=1}^{\infty}\frac{d}{dt_{f}}Q_{t_{f}}(N)\Big{(}-{\rm
erfc}\left(\frac{x_{ab}-aN}{\sqrt{2}\delta\sqrt{N}}\right)+e^{\frac{2ax_{ab}}{\delta^{2}}}{\rm
erfc}\left(\frac{aN+x_{ab}}{\sqrt{2}\delta\sqrt{N}}\right)\Big{)}\end{array}$
(59)

with

$\frac{d}{dt_{f}}Q_{t_{f}}(N)=-\mathcal{L}_{\alpha}\left(\frac{N-\frac{t_{f}}{\langle\tau\rangle}}{(\frac{t_{f}}{\bar{t}})^{1/\alpha}}\right)\left(\frac{\bar{t}}{t_{f}}\right)^{1/\alpha}\frac{1}{\alpha
t_{f}}+\left(\frac{t_{f}}{\bar{t}}\right)^{1/\alpha}\frac{d}{dt_{f}}\mathcal{L}_{\alpha}\left(\frac{N-\frac{t_{f}}{\langle\tau\rangle}}{(\frac{t_{f}}{\bar{t}})^{1/\alpha}}\right).$

Here ${\rm erfc}(z)$ denotes the complementary error function, i.e., ${\rm
erfc}(z)=1-{\rm erf}(z)$ with error function ${\rm
erf}(z)=\frac{2}{\sqrt{\pi}}\int_{0}^{z}\exp(-\tau^{2})d\tau$. Eq. 59 is
confirmed in Fig. (8) displaying a perfect match.

Figure 7: Plot of concentration $C(x,t)$ using Eq. (56) for different biases
and locations of absorbing boundary conditions. Here we choose $t=1000$,
$\tau_{0}=0.1$, $\alpha=3/2$, and $10^{7}$ realizations. We use
$a=2,x_{ab}=7300$; $a=1,x_{ab}=3700$; $a=1/2,x_{ab}=1850$; and
$a=3/10,x_{ab}=1100$ for subplot (a) to (d), respectively. The Solution of Eq.
(3) without absorbing condition is plotted by the dashed lines.
Figure 8: Behaviors of the first passage time PDF generated by trajectories of
particles for various $a$. The parameters are the same as in Fig. 7. The
different lines describe the theoretical result obtained from Eq. (59).

##  5 Three limiting laws of the positional distribution

Using Eq. (3), we derive three limiting laws that characterize the CTRW model
under different conditions: when $S$ tends towards infinity, when $D$
approaches infinity, and when both $S$ and $D$ are constants. More exactly,
these laws include the Lévy stable distribution showing the statistics when
$x-Vt\propto(St)^{1/\alpha}$, Gaussian distribution in the context of
$x-Vt\propto\sqrt{Dt}$, and a general expression that encompasses all values
of $D,V$ and $S$.

Taking Laplace transform on Eq. (29) with respect to $t$, we have

$\mathcal{P}(k,s)=\frac{1}{s+Dk^{2}+ikV-S(-ik)^{\alpha}}.$ (60)

We are interested in the statistics of $x$ in the long time limit, which means
that both $s$ and $k$ go to zero in Eq. (60). Based on the above equation, the
mentioned three laws will be discussed.

###  5.1 Lévy stable distribution when $S\to\infty$

Let $S\to\infty$, Eq. (60) yields

$\mathcal{P}(k,s)\sim\frac{1}{s+iVk-S(-ik)^{\alpha}},$ (61)

where we dropped the term $Dk^{2}$ due $Dk^{2}\ll iVk,~{}S(-ik)^{\alpha}$. The
inverse Laplace-Fourier transforms of Eq. (61) lead to

$\mathcal{P}(x,t)\sim\frac{1}{(St)^{1/\alpha}}\mathcal{L}_{\alpha}\left(\frac{x-Vt}{(St)^{1/\alpha}}\right).$
(62)

Integrating Eq. (62) from negative infinity to positive infinity implies that
$\mathcal{P}(x,t)$ is a normalized density. See also the equivalent expression
given in Eq. (19). This is illustrated in Fig. 9, where the PDF of
$\xi=(x-Vt)/(St)^{1/\alpha}$ with $\alpha=1.5$ is plotted. Obviously, the PDF
of $\xi$ is asymmetric with respect to $\xi$, whose two tails show two
different dynamic behaviors, namely the right-hand side of the tail tends to
zero rapidly but the other one decays slowly like a power law; see the inset
in Fig. 9. More precisely, based on Eq. (12), we get
$\mathcal{L}(\xi)\sim(-\xi)^{-1-\alpha}/\Gamma(-\alpha)$, being the same as
the tail of the symmetric Lévy stable distribution, for $\xi\to-\infty$. Thus,
when $q>\alpha$, the integral
$\int_{-\infty}^{\infty}|\xi|^{q}\mathcal{L}_{\alpha}(\xi){\rm d}\xi$
diverges. This means that Eq. (62) does not give any information to the MSD,
which will be discussed in Appendix B.

Figure 9:  PDF $\mathcal{L}_{\alpha}(\xi)$ versus
$\xi=(x-Vt)/(a(St)^{1/\alpha})$ for a biased CTRW model with the waiting time
distribution Eq. (4) and Gaussian distribution of the step length Eq. (6). The
red line is the theoretical prediction Eq. (62), showing the limit of typical
fluctuations when $x-Vt$ is of the order of $t^{1/\alpha}$. The corresponding
simulation result, represented by symbols ‘$\circ$’, is obtained by averaging
$5\times 10^{6}$ trajectories. Here the parameters are $t=1000$, $a=2$,
$\sigma=0.07$, and $\tau_{0}=0.1$. We will show how the limiting law Eq. (62)
fails for a finite time $t$ or a large $D$ and give a general law, together
with a slow convergence; see Fig. 2 and below discussions.

###  5.2 Gaussian distribution when $D\to\infty$

Note that Eq. (62) is independent of $D$, and when $D$ is large, we expect
this approximation to fail for a large but finite $t$. Here we focus on the
case when $s+ikV\propto k^{2}$. As previously discussed, it is related to the
linear response regime. In this limit, according to Eq. (60)

$\mathcal{P}(k,u)\sim\frac{1}{s+Dk^{2}+ikV}.$ (63)

By inversion, using the shifting property of the inverse Fourier transform,
yields the limiting distribution of $x$. The scaling form of
$\varsigma=(x-Vt)/\sqrt[]{2Dt}$ gives the Gaussian distribution with mean zero
and variance one

$\mathcal{P}(\varsigma)\sim\frac{1}{\sqrt{2\pi}}\exp\left(-\frac{\varsigma^{2}}{2}\right);$
(64)

see Fig. 10. As expected, the simulations are consistent with the theoretical
result Eq. (64) for a large $D$. Recall that in Figs. 9 and 10, we used the
same observation time $t$, i.e., $t=1000$, but propagators are totally
different (one is the Lévy stable law and the other is Gaussian distribution).
It indicates that these laws are determined by the relationship between $D$
and $S$. The current challenge is to determine a universal form that can
represent both of the mentioned scalings.

Figure 10:  The dynamic of $\mathcal{P}(\varsigma)$ with scaling
$\varsigma=(x-Vt)/\sqrt[]{2Dt}$ for different observation times $t$. Here we
choose $\alpha=1.5$, $a=1$, $\sigma=11$, $\tau_{0}=0.1$, and
$\langle\tau\rangle=0.3$. These parameters suggest that the process can be
approximated as nearly Gaussian, with $D$ being approximately $201.7$ and $S$
being approximately $2.3$. The full red line depicts the Gaussian distribution
Eq. (64), which represents the behaviors when $x-Vt$ is of the order of
$\sqrt{t}$. On the other hand, the symbols correspond to simulation results
obtained by averaging $10^{6}$ trajectories.

###  5.3 Typical fluctuations

We have discussed that under certain conditions the positional distribution
follows Lévy stable distribution Eq. (62) or Gaussian distribution Eq. (64),
determined by the relation between $D$ and $S$. Therefore, it would be
beneficial to acquire a universal law that applies to the aforementioned
scenarios. Based on Eq. (60), we have

$\mathcal{P}(x,t)=\frac{1}{\sqrt{4\pi
D}S^{1/\alpha}t^{1/\alpha+1/2}}\int_{-\infty}^{\infty}\mathcal{L}_{\alpha}\left(\frac{y}{(St)^{1/\alpha}}\right)\exp\left(-\frac{(x-Vt-y)^{2}}{4Dt}\right)dy.$
(65)

We denote Eq. (65), i.e., the solution of Eq. (3), as typical fluctuations
describing the central of the positional distribution. As mentioned before, in
the long time limit, Eq. (65) leads to the Lévy stable distribution Eq. (62).
Figure 11 demonstrates the validity of Eq. (65) for various values of $t$. For
comparison, we also plot the solution Eq. (16) in Fig. 12. Eq. (16) agrees
with Eq. (65) for large $x$, it rapidly approaches zero as $x$ approaches
zero, whereas Eq. (65) exhibits a cutoff at the tail. Therefore, the primary
difference between the two equations lies in small values of $x$.

Figure 11:  Graph of densities $\mathcal{P}(\xi,t)$ with
$\xi=(x-at/\langle\tau\rangle)/(a(t/\overline{t})^{1/\alpha})$ for different
times $t$. The red solid lines represent the analytical scaling result
obtained from Eq. (65) and the corresponding blue symbols are simulations of
$5\times 10^{6}$ trajectories. The dashed black line is the limiting law
$\mathcal{L}_{\alpha}(\xi)$ predicted by Eq. (62). Just as the figure shows,
with the increase of the observation time $t$, the PDF $\mathcal{P}(\xi,t)$
tends to $\mathcal{L}_{\alpha}(\xi)$ slowly. The parameters are $\alpha=1.5$,
$\tau_{0}=0.1$, $a=1$, $\langle\tau\rangle=0.3$ and $\sigma=1/\sqrt{2}$. See
also discussions in Appendix A.
Figure 12: Comparison between Eq. (16) and Eq. (65) for different forces. The
symbols are the plot of Eq. (16) and the solid lines describe the solution of
Eq. (3) or Eq. (65). Here we choose $t=1000$, $t_{0}=0.1$, $\alpha=1.5$, and
$\delta=1$.

To summarize, typical fluctuations Eq. (65), giving the information when
$x-Vt\propto t^{1/\alpha}$, are valid for numerous $D$, $V$, $S$ and $t$. One
may wonder “why Eq. (65) is more useful than Eq. (19)?”. The reason is as
follows: Mathematically, Eq. (19) is valid under the condition that
$gg(y,t,\sigma,a,\alpha)\to 0$. This implies that the observation time $t$
should be especially long, for example when $a$ is small and
$\sigma\to\infty$. While, Eq. (65) avoids this, this means that the condition
reduces to $\langle N\rangle\propto t/\langle\tau\rangle$.

##  6 Conclusion and discussion

In the context of the CTRW model, if the displacement follows a narrow
distribution with a non-zero mean and the waiting time has an infinite mean,
the fractional operator of the diffusion equation is linked to the Riemann-
Liouville time derivative, as discussed in Refs. [4, 33]. However, nearly two
decades later, the fractional-space advection-diffusion equation was developed
for situations in which waiting times possess a finite mean and an infinite
variance, together with the mentioned displacements [53]. As an extension of
work [53], here we demonstrated further evidence of Eq. (3), gave more
applications, and showed limiting laws of the kinetic equation, through our
in-depth explanations. At the same time, theoretical predictions are checked
by simulations of the CTRW model.

An issue here is that for Lévy flights, with a fat-tailed distribution of jump
lengths, the method of images was shown to lead to wrong results, due to the
non-local behavior [90, 91, 92]. One may wonder, whether or not the method of
images will work for our process which does not include any fat-tailed jump
length distribution. We show that generally, it does not work, unless for the
case when the bias is small, then the image method is a valid approximation.
For that, a subordination method was used. This method holds in general beyond
the first passage problem, and it is also very different if compared to the
subordination method used for $0<\alpha<1$. In practice we use the image
method on the conditional probability $\chi_{N}(x)$, instead of the positional
distribution $\mathcal{P}(x,t)$. In particular, for a free particle, we assume
Eq. (54) has a free boundary condition. In that sense, $\chi^{*}_{N}(x)$ Eq.
(55) reduces to $\chi_{N}(x)$ given in Eq. (9). The mentioned strategy can
also be extended to a much more general case, for example, the diffusion of
particles on a finite domain. For applications, we also discussed the two-
dimensional diffusion equation and breakthrough through curves with time-
dependent bias and variance.

We demonstrated that the solution of Eq. (3) is obtained by convolving the
Gaussian distribution with the Lévy stable distribution, see Eq. (65). Here we
denote the solution of Eq. (3) as typical fluctuations of CTRW describing a
wide range of $D$, $V$, and $S$. In the language of the biased CTRW, the Lévy
stable distribution stems from the PDF of the number of renewals, and the
Gaussian distribution comes from the conditional distribution $\chi_{N}(x)$.
When the bias is strong, i.e., $S\to\infty$, the diffusion term in Eq. (3) can
be ignored and the solution is an asymmetric Lévy stable distribution Eq.
(62). For the CTRW model, this limit is achieved when the ratio
$\sigma^{2}/a^{2}$ is finite and $t\to\infty$ or $\sigma^{2}/a^{2}\to 0$ with
a finite $t$. On the contrary, when $S\to 0$, the spreading packet follows
Gaussian distribution. It indicates that when $\sigma^{2}/a^{2}\gg 1$, the
well-known Eq. (62) is not useful unless $t$ is extremely long. Furthermore,
we provide a parameter, denoted as $\overline{\sigma}^{2}$, in Eq. (20) to
characterize the convergence.

There are still unanswered questions that require addressing. For example, in
the simulation of particle trajectories, generally, we first generate waiting
times and displacements. Then, upon the completion of the total observation
time, denoted as $t$, we obtain positional statistics. However, when the
observation time is long and the number of realizations is large, such as
$10^{7}$, running codes may require several days to complete. The problem now
is whether we can generate particle statistics based on Eq. (21). More
precisely, we first generate a variable drawn from the Lévy stable
distribution $\mathcal{L}(y)$, and then use it to generate the position at
time $t$ according to a Gaussian distribution. This method seems fine if we
are only interested in the positional distribution, but it fails as expected
when we focus on the MSD. Solving this problem, i.e. formulating Langevin
paths, is a matter that will be considered in the future.

## Acknowledgments

W.W. is supported by the National Natural Science Foundation of China under
Grant No. 12105243 and the Zhejiang Province Natural Science Foundation
LQ$22$A$050002$. E.B. acknowledges the Israel Science Foundations grant
$1614/21$.

##  Appendix A Additional discussions on Fig. 2

From simulations of typical fluctuations in Fig. 2, we can see that we have
$\overline{\sigma}^{2}\simeq 0.1115\sigma^{2}/a^{2}$ defined by Eq. (20) with
$t=1000$, $\alpha=1.5$, and $\tau_{0}=0.1$. In the particular case
$\sigma=0.1/\sqrt{2}$, from Eq. (20) we get a small
$\overline{\sigma}^{2}\simeq 5.57\times 10^{-4}$, this is the reason why Eq.
(16) tends to the limit theorem Eq. (19) or (62) as shown by the black dashed
line in Fig. 2. However, when $\sigma=16/\sqrt{2}$, $\sigma^{2}/a^{2}=128$ and
then $\overline{\sigma}^{2}=14.27$, which is certainly not a small number if
compared with $5.57\times 10^{-4}$. Thus even though the average of renewals
$\langle N(t)\rangle\sim t/\langle\tau\rangle=1000/0.3\approx 3333$, we cannot
say that the limit of $t\to\infty$ is reached, and indeed, just as the bottom
line in Fig. 2 shows, we see for this case nearly Gaussian behavior. For
asymmetry breaking properties, one choice is to consider
$\overline{\sigma}^{2}$. This means that when $\overline{\sigma}^{2}\to 0$,
the corresponding density is asymmetric. While, for a large
$\overline{\sigma}^{2}$, the opposite situation emerges.

##  Appendix B Fractional moments

As mentioned before, Eq. (3) fails to estimate MSD of the CTRW model.
Recently, we give a way to compute the MSD using infinite densities [68]. Here
we show that infinite density frameworks can also used to calculate the
fractional moments in an asymptotic sense. For the sake of completeness, first
let us concentrate on the variance of the walk’s displacement, called MSD,
which is a measure of the deviation of the position of the particle with
respect to the position over time. Suppose that $x=0$ is the initial position
of particles and $N$ is the number of renewals between $0$ and $t$. Thus, the
first and the second order moment of the position are given by

$\langle x\rangle=\left\langle\sum_{j=1}^{N}\chi_{i}\right\rangle$ (66)

and

$\langle
x^{2}\rangle=\left\langle\left(\sum_{j=1}^{N}\chi_{i}\right)^{2}\right\rangle,$
(67)

respectively. Here $\chi_{i}$ denotes the displacement of the particle in the
$i$-th jump and $\langle\cdot\rangle$ describes the average over
$\chi_{1},\chi_{2},\dots,\chi_{N}$. Therefore, the variance is

$V\\!ar(x)=\langle x^{2}\rangle-\langle x\rangle^{2}=\langle
N\rangle(\langle\triangle x^{2}\rangle-\langle\triangle x\rangle^{2})+(\langle
N^{2}\rangle-\langle N\rangle^{2})\langle\triangle x\rangle^{2},$ (68)

where we assumed that $\chi_{i}$ are IID random variables, and
$\langle\triangle x\rangle$ and $\langle\triangle x^{2}\rangle$ correspond to
the first and the second order moments of $\chi_{i}$, respectively. Let us
take the step length distribution as Eq. (6), thus we have $\langle\triangle
x\rangle=a$ and $\langle\triangle x^{2}\rangle=a^{2}+\sigma^{2}$. To obtain
$V\\!ar(x)$, we need to calculate the first and the second moments of $N$.
Using the previous results given in Ref. [59], we have

$\langle\widehat{N}(s)\rangle=\frac{\widehat{\phi}(s)}{s(1-\widehat{\phi}(s))}$
(69)

and

$\langle\widehat{N}^{2}(s)\rangle=\frac{\widehat{\phi}(s)(1+\widehat{\phi}(s))}{s(1-\widehat{\phi}(s))^{2}}.$
(70)

Utilizing Eqs. (69) and (70), Eq. (68) gives

$V\\!ar(x)\sim\frac{2a^{2}\tau_{0}^{\alpha}}{(2-\alpha)(3-\alpha)\langle\tau\rangle^{3}}t^{3-\alpha}+\sigma^{2}\frac{t}{\langle\tau\rangle},$
(71)

which is plotted in Fig. 13. The term $t^{3-\alpha}$ becomes the leading term
of $V\\!ar(x)$ as we increase of the observation time $t$; See also
discussions in Ref. [68]. When the observation time is not very large, the
linear term $\sigma^{2}t/\langle\tau\rangle$ wins since the spreading packet
nearly follows Gaussian distribution; see the dash-dotted line in Fig. 13. For
a fixed observation time $t$, we can find an interesting phenomenon depicting
the competition between anomalous diffusion and normal diffusion, which is
determined by $a$, $\sigma$, and $t$. According to Eq. (71), the transition
point is

$t^{*}\sim\left((2-\alpha)(3-\alpha)\frac{\sigma^{2}\langle\tau\rangle^{2}}{2a^{2}\tau_{0}^{\alpha}}\right)^{\frac{1}{2-\alpha}};$
(72)

see the magenta vertical lines in Fig. 13. When $|a|\to 0$, $t^{*}$ goes to
infinity. It indicates that when the bias is weak, the process needs a long
observation time $t$ to exhibit anomalous behaviors. On the other hand, when
$a$ approaches infinity, the value of $t^{*}$ tends to zero, causing the
transition from normal diffusion to anomalous diffusion rapidly.

Figure 13:  MSD for observation time $t=500$ [left top panel], $t=2000$ [right
top panel], $t=4000$ [right bottom panel], $t=10^{5}$ [left bottom panel],
where the parameters are $\alpha=1.5$, $\sigma=3/\sqrt{2}$, $a=1$, and
$\tau_{0}=0.1$. The solid lines are the theoretical result Eq. (71) and the
dashed curves are the dominating term of Eq. (71), i.e.,
$(2a^{2}\tau_{0}^{\alpha})/((2-\alpha)(3-\alpha)\langle\tau\rangle^{3})t^{3-\alpha}$.
The symbols are simulations obtained by averaging $10^{4}$ realizations.

When $q=4$, the corresponding moment is

$\langle(x-\langle x\rangle)^{4}\rangle=\langle x^{4}\rangle+6\langle
x\rangle^{2}\langle x^{2}\rangle-4\langle x\rangle\langle
x^{3}\rangle-3\langle x\rangle^{4}$ (73)

Using Eq. (8), the integer order moments are easy to get. Then substituting
them into Eq. (73) gives

$\langle(x-\langle x\rangle)^{4}\rangle\sim a^{4}\langle(N-\langle
N\rangle)^{4}\rangle+6a^{2}\sigma^{2}(\langle N^{3}\rangle+\langle
N\rangle^{3}-2\langle N\rangle\langle N^{2}\rangle)$ (74)

Now the calculation of $\langle(x-\langle x\rangle)^{4}\rangle$ is transformed
into the moments of the number of renewals. After some simple arithmetics, we
have

$\displaystyle{\begin{array}[]{ll}\langle(x-\langle
x\rangle)^{4}\rangle&\sim\frac{4a^{4}b_{\alpha}}{|\Gamma(1-\alpha)|\langle\tau\rangle^{5}(5-\alpha)(4-\alpha)}t^{5-\alpha}+\frac{6a^{2}\sigma^{2}b_{\alpha}}{\langle\tau\rangle^{4}}\Big{(}\frac{18}{\Gamma(5-\alpha)}-\frac{8}{\Gamma(4-\alpha)}+\frac{1}{\Gamma(3-\alpha)}\Big{)}t^{4-\alpha}\\\
&~{}~{}~{}+\frac{12a^{4}b_{\alpha}}{\langle\tau\rangle^{4}}\Big{(}\frac{16}{\Gamma(5-\alpha)}-\frac{4}{\Gamma(4-\alpha)}+\frac{1}{\Gamma(3-\alpha)}\Big{)}t^{4-\alpha}.\end{array}}$
(75)

Our next step is to explore $q$ order absolute moments of
$\epsilon=x-at/\langle\tau\rangle$, defined by

$\langle|\epsilon|^{q}\rangle=\int_{-\infty}^{\infty}|\epsilon|^{q}P(\epsilon,t){\rm
d}\epsilon$ (76)

with $q>0$. Let us start from the calculation of low order moments, i.e.,
$q<\alpha$. The low order moments can be obtained by the limit of typical
fluctuations Eq. (19), which reads

$\langle|\epsilon|^{q}\rangle\sim
a^{q}(t/\overline{t})^{q/\alpha}\int_{-\infty}^{\infty}y^{q}\mathcal{L}_{\alpha}(y){\rm
d}y.$ (77)

When $q<\alpha$, the integral
$\int_{-\infty}^{\infty}y^{q}\mathcal{L}_{\alpha}(y){\rm d}y$ is a finite
number which can be evaluated by the asymptotic behaviors of Lévy distribution
$\mathcal{L}_{\alpha}(y)$. We can check that for $0<q<\alpha$,
$\langle|\epsilon|^{q}\rangle\sim t^{q/\alpha}$ is always sublinear in $t$.
Here we would like to consider two special cases: for $q\to 0$, the zeroth
moment of the variable $\epsilon$ is, evidently, one; the other case is the
Gaussian limit, the linear time dependence of the MSD is recovered, namely
$\lim_{\alpha,q\to 2}\langle|\epsilon|^{q}\rangle\sim t$. For $q>\alpha$, Eq.
(77) diverges. It indicates that the normalized density Eq. (19) can not give
a valid prediction for high order moments; see Eq. (68).

Now we consider high-order moments. As mentioned before, for this case,
typical fluctuations does not work. In [68], we demonstrate the existence of
an additional limiting law when the second time scale is taken into account.
This law, denoted as Eq. (78), characterizes the scaling behavior when
$x-at/\langle\tau\rangle$ is of the order of $t/\langle\tau\rangle$. Namely,
the density of $\eta=(x-at/\langle\tau\rangle)/(at/\langle\tau\rangle)$ is

$P(\eta,t)\sim\frac{b_{\alpha}t^{1-\alpha}}{\Gamma(-\alpha)\langle\tau\rangle}(-\eta)^{-\alpha-1}\Big{(}1-\frac{1-\alpha}{\alpha}\eta\Big{)}$
(78)

with $-1<\eta<0$. See numerical simulations and discussions in Ref. [68].

As mentioned, this situation for $q=2$ is especially important in physical and
biological experiments [94, 4]. After rescaling Eq. (78) by
$\epsilon=x-at/\langle\tau\rangle$, it can be seen that $|\epsilon|^{q}$ is
integrable with respect to this non-normalized state. So this can be used to
get high order moments, i.e.,

$\langle|\epsilon|^{q}\rangle\sim\frac{a^{q}q\tau_{0}^{\alpha}}{\langle\tau\rangle^{q+1}(q-\alpha)(q-\alpha+1)}t^{q+1-\alpha}.$
(79)

This demonstrates that, in the long time regime, the asymptotic distribution
of $\epsilon$ is broad with a slowly decaying tail; see Eq. (78). Note that
when $t$ is not sufficient large, the correction term is necessary to be added
since the packet of particles is Gaussian; see left top panel in Fig. 13.

To summarize, the limiting law Eq. (19) does describe the low order moments.
The rare fluctuations, predicted by Eq. (78), give the information to the
events when $\epsilon$ is of the order of $t$. In the long time limit, the
results of $q$ order moments can be generalized to

$\langle|\epsilon|\rangle\sim t^{\varrho(q)},$ (80)

where $\varrho(q)=q/\alpha$ if $q<\alpha$, and $\varrho(q)=q+1-\alpha$ if
$q>\alpha$.

##  Appendix C Plot of concentration with a weak bias

We plot Eq. (52) obtained from the image method for a weak bias; see Fig. 14.
For this case, the asymmetric term $S$ loses its role and the process is
nearly normal. As mentioned in the main text, if the bias is strong, Eq. (52)
breaks down.

Figure 14: Plot of concentration Eq. (52) for different absorbing boundary
conditions, i.e., $x_{ab}=1200,1000,800,600$ for (a), (b), (c), (d)
respectively. The parameters for CTRW simulations are $10^{7}$ realizations,
$t=800$, $\alpha=1.5$, $\sigma=1$, $a=0.3$, and $\tau_{0}=0.1$. Red solid
lines demonstrate the theoretical predication Eq. (52) with $\mathcal{P}(x,t)$
being the solution of Eq. (3) and dashes lines are the solution of Eq. (3),
i.e., without absorbing conditions, for comparison.

## References

  * [1] Shlesinger M F 1974 Asymptotic solutions of continuous-time random walks.  _J. Stat. Phys._ 10 421–434
  * [2] Montroll E W, Weiss G H 1965 Random walks on lattices. II.  _J. Math. Phys._ 6 167–181
  * [3] Bouchaud J P, Georges A 1990 Anomalous diffusion in disordered media: statistical mechanisms, models and physical applications.  _Phys. Rep._ 195 127–293
  * [4] Metzler R, Klafter J 2000 The random walk’s guide to anomalous diffusion: a fractional dynamics approach.  _Phys. Rep._ 339 1–77
  * [5] Chaudhuri P, Berthier L, Kob W 2007 Universal nature of particle displacements close to glass and jamming transitions.  _Phys. Rev. Lett._ 99 060604
  * [6] Weigel A V, Simon B, Tamkun M M, Krapf D 2011 Ergodic and nonergodic processes coexist in the plasma membrane as observed by single-molecule tracking.  _PNAS_ 108 6438—6443
  * [7] Meroz Y, Sokolov I M 2015 A toolbox for determining subdiffusive mechanisms.  _Physics Reports_ 573 1–29
  * [8] Höfling F, Franosch T 2013 Anomalous transport in the crowded world of biological cells.  _Rep. Prog. Phys._ 76 046602
  * [9] Gradenigo G, Bertin E, Biroli G 2016 Field-induced superdiffusion and dynamical heterogeneity.  _Phys. Rev. E_ 93 060105
  * [10] Kutner R, Masoliver J 2017 The continuous time random walk, still trendy: fifty-year history, state of art and outlook.  _Eur. Phys. J. B_ 90 50
  * [11] Comolli A, Dentz M 2018 Impact of diffusive motion on anomalous dispersion in structured disordered media: From correlated Lévy flights to continuous time random walks.  _Phys. Rev. E_ 97 052146
  * [12] Barkai E, Burov S 2020 Packets of diffusing particles exhibit universal exponential tails.  _Phys. Rev. Lett._ 124 060603
  * [13] Michelitsch T M, Polito F, Riascos A P 2020 Biased continuous-time random walks with Mittag-Leffler jumps.  _Fractal and Fractional_ 4 51
  * [14] Vot F L, Abad E, Metzler R, Yuste S B 2020 Continuous time random walk in a velocity field: role of domain growth, Galilei-invariant advection-diffusion, and kinetics of particle mixing.  _New J. Phys._ 22 073048
  * [15] Wang W, Barkai E, Burov S 2020 Large deviations for continuous time random walks.  _Entropy_ 22 697
  * [16] Pacheco-Pozo A, Sokolov I M 2021 Large deviations in continuous-time random walks.  _Phys. Rev. E_ 103 042116
  * [17] Liu J, Zhu P, Bao J D, Chen X 2022 Strong anomalous diffusive behaviors of the two-state random walk process.  _Phys. Rev. E_ 105 014122
  * [18] Michelitsch T M, Polito F, Riascos A P 2022 Asymmetric random walks with bias generated by discrete-time counting processes.  _Commun. Nonlinear Sci. Numer. Simul._ 109 106121
  * [19] Burov S, Wang W, Barkai E Exponential tails and asymmetry relations for the spread of biased random walks.  ArXiv:2209.03410
  * [20] Vitali S, Paradisi P, Pagnini G 2022 Anomalous diffusion originated by two markovian hopping-trap mechanisms.  _J. Phys. A: Math. Theor._ 55 224012
  * [21] Dentz M, Kirchner J W, Zehe E, Berkowitz B 2023 The role of anomalous transport in long-term, stream water chemistry variability.  _Geophys. Res. Lett._ 50 e2023GL104207
  * [22] Scher H, Montroll E W 1975 Anomalous transit-time dispersion in amorphous solids.  _Phys. Rev. B_ 12 2455–2477
  * [23] David S A, Valentim C A, Debbouche A 2022 Fractional modeling applied to the dynamics of the action potential in cardiac tissue.  _Fractal Fract._ 6
  * [24] Deng W 2009 Finite element method for the space and time fractional Fokker-Planck equation.  _SIAM J. Numer. Anal._ 47 204–226
  * [25] Sylvain T T A, Patrice E A, Marie E E J, Pierre O A, Hubert B B G 2022 A unified three-dimensional extended fractional analytical solution for air pollutants dispersion.  _Fractals_ 30 2250065
  * [26] Suleiman K, Song Q, Zhang X, Liu S, Zheng L 2022 Anomalous diffusion in a circular comb with external velocity field.  _Chaos, Solitons & Fractals_ 155 111742
  * [27] Yan Z 2022 New integrable multi-Lévy-index and mixed fractional nonlinear soliton hierarchies.  _Chaos, Solitons & Fractals_ 164 112758
  * [28] Evangelista L R, Lenzi E K 2023 _Fractional Anomalous Diffusion_ , 189–236 (Cham: Springer International Publishing)
  * [29] Tang Y, Gharari F, Arias-Calluari K, Alonso-Marroquin F, Najafi M N 2023 Variable order porous media equations: Application on modeling the S&P500 and bitcoin price return.  ArXiv:2309.04206
  * [30] Uchaikin V, Sibatov R 2013 _Fractional Kinetics in Solids_ (Singapore: World Scientific)
  * [31] Eraso-Hernandez L, Riascos A, Michelitsch T, Wang-Michelitsch J 2021 Random walks on networks with preferential cumulative damage: Generation of bias and aging.  _J. Stat. Mech: Theory Exp._ 2021 063401
  * [32] Metzler R, Barkai E, Klafter J 1999 Anomalous diffusion and relaxation close to thermal equilibrium: a fractional Fokker-Planck equation approach.  _Phys. Rev. Lett._ 82 3563–3567
  * [33] Barkai E 2001 Fractional Fokker-Planck equation, solution, and application.  _Phys. Rev. E_ 63 046118
  * [34] Oldham K B, Spanier J 1974 _The Fractional Calculus_ (Academic Press, New York)
  * [35] Benson D A, Wheatcraft S W, Meerschaert M M 2000 The fractional-order governing equation of Lévy motion.  _Water Resour. Res._ 36 1413–1423
  * [36] Samko S G, Kilbas A A, Marichev O I 1993 _Fractional integrals and derivatives_ (Gordon and Breach Science Publishers, Yverdon)
  * [37] Podlubny I 1999 _Fractional Differential Equations_ (Academic Press, Inc., San Diego)
  * [38] Meerschaert M M, Sikorskii A 2012 _Stochastic Models for Fractional Calculus_ , vol. 43 of _De Gruyter Studies in Mathematics_ (Walter de Gruyter & Co., Berlin)
  * [39] Estrada-Rodriguez G, Gimperlein H, Painter K J, Stocek J 2019 Space-time fractional diffusion in cell movement models with delay.  _Math. Models Methods Appl. Sci._ 29 65–88
  * [40] Zheng X, Wang H 2020 An error estimate of a numerical approximation to a hidden-memory variable-order space-time fractional diffusion equation.  _SIAM J. Numer. Anal._ 58 2492–2514
  * [41] Hufnagel L, Brockmann D, Geisel T 2004 Forecast and control of epidemics in a globalized world.  _PNAS_ 101 15124–15129
  * [42] Lomholt M A, Ambjörnsson T, Metzler R 2005 Optimal target search on a fast-folding polymer chain with volume exchange.  _Phys. Rev. Lett._ 95 260603
  * [43] Brockmann D, Hufnagel L, Geisel T 2006 The scaling laws of human travel.  _Nature_ 439 462–465
  * [44] Benson D A, Schumer R, Meerschaert M M, Wheatcraft S W 2001 Fractional dispersion, Lévy motion, and the MADE tracer tests.  _Transp. Porous Media_ 42 211–240
  * [45] Zaburdaev V, Denisov S, Klafter J 2015 Lévy walks.  _Rev. Mod. Phys._ 87 483–530
  * [46] Zhang Y, Meerschaert M M, Neupauer R M 2016 Backward fractional advection dispersion model for contaminant source prediction.  _Water Resour. Res._ 52 2462–2473
  * [47] Guo Y, Wang L 2021 Heat current flows across an interface in two-dimensional lattices.  _Phys. Rev. E_ 103 052141
  * [48] Sin C S 2023 Cauchy problem for fractional advection-diffusion-asymmetry equations.  _Results in Mathematics_ 78 111
  * [49] Saichev A I, Zaslavsky G M 1997 Fractional kinetic equations: solutions and applications.  _Chaos_ 7 753–764
  * [50] Gorenflo R, Mainardi F, Vivoli A 2007 Continuous-time random walk and parametric subordination in fractional diffusion.  _Chaos Solitons Fractals_ 34 87–103
  * [51] Weron K, Jurlewicz A, Magdziarz M, Weron A, Trzmiel J 2010 Overshooting and undershooting subordination scenario for fractional two-power-law relaxation responses.  _Phys. Rev. E_ 81 041123
  * [52] Dybiec B, Gudowska-Nowak E 2010 Subordinated diffusion and continuous time random walk asymptotics.  _Chaos: An Interdisciplinary Journal of Nonlinear Science_ 20 043129
  * [53] Wang W, Barkai E 2020 Fractional advection-diffusion-asymmetry equation.  _Phys. Rev. Lett._ 125 240606
  * [54] Chechkin A, Sokolov I M 2021 Relation between generalized diffusion equations and subordination schemes.  _Phys. Rev. E_ 103 032133
  * [55] Zhou T, Trajanovski P, Xu P, Deng W, Sandev T, Kocarev L 2022 Generalized diffusion and random search processes.  _J. Stat. Mech.: Theory Exp._ 2022 093201
  * [56] Wang X, Chen Y 2022 Ergodic property of random diffusivity system with trapping events.  _Phys. Rev. E_ 105 014106
  * [57] Defaveri L, dos Santos M A F, Kessler D A, Barkai E, Anteneodo C 2023 Non-normalizable quasiequilibrium states under fractional dynamics.  _Phys. Rev. E_ 108 024133
  * [58] Berkowitz B, Klafter J, Metzler R, Scher H 2002 Physical pictures of transport in heterogeneous media: Advection-dispersion, random-walk, and fractional derivative formulations.  _Water Resour. Res._ 38 9
  * [59] Godrèche C, Luck J M 2001 Statistics of the occupation time of renewal processes.  _J. Stat. Phys._ 104 489–524
  * [60] Wang W, Schulz J H P, Deng W H, Barkai E 2018 Renewal theory with fat-tailed distributed sojourn times: Typical versus rare.  _Phys. Rev. E_ 98 042139
  * [61] Metzler R, Jeon J H, Cherstvy A G, Barkai E 2014 Anomalous diffusion models and their properties: Non-stationarity, non-ergodicity, and ageing at the centenary of single particle tracking.  _Phys. Chem. Chem. Phys._ 16 24128–24164
  * [62] Margolin G, Berkowitz B 2002 Spatial behavior of anomalous transport.  _Phys. Rev. E_ 65 031101
  * [63] Levy M, Berkowitz B 2003 Measurement and analysis of non-Fickian dispersion in heterogeneous porous media.  _J. Contam. Hydrol._ 64 203–226
  * [64] Schroer C F E, Heuer A 2013 Anomalous diffusion of driven particles in supercooled liquids.  _Phys. Rev. Lett._ 110 067801
  * [65] Luo L, Tang L H 2015 Sample-dependent first-passage-time distribution in a disordered medium.  _Phys. Rev. E_ 92 042137
  * [66] Nissan A, Dror I, Berkowitz B 2017 Time-dependent velocity-field controls on anomalous chemical transport in porous media.  _Water Resour. Res._ 53 3760–3769
  * [67] Vezzani A, Barkai E, Burioni R 2019 Single-big-jump principle in physical modeling.  _Phys. Rev. E_ 100 012108
  * [68] Wang W, Vezzani A, Burioni R, Barkai E 2019 Transport in disordered systems: The single big jump approach.  _Phys. Rev. Res._ 1 033172
  * [69] Zhang C, Hu Y, Liu J 2022 Correlated continuous-time random walk with stochastic resetting.  _J. Stat. Mech: Theory Exp._ 2022 093205
  * [70] Fogedby H C 1994 Langevin equations for continuous time Lévy flights.  _Phys. Rev. E_ 50 1657–1660
  * [71] Magdziarz M, Weron A, Weron K 2007 Fractional Fokker-Planck dynamics: Stochastic representation and computer simulation.  _Phys. Rev. E_ 75 016708
  * [72] Fedotov S, Han D 2023 Population heterogeneity in the fractional master equation, ensemble self-reinforcement, and strong memory effects.  _Phys. Rev. E_ 107 034115
  * [73] Kotulski M 1995 Asymptotic distributions of continuous-time random walks: A probabilistic approach.  _J. Stat. Phys._ 81 777–792
  * [74] Burioni R, Gradenigo G, Sarracino A, Vezzani A, Vulpiani A 2013 Rare events and scaling properties in field-induced anomalous dynamics.  _J. Stat. Mech. Theory Exp._ 2013 P09022
  * [75] Metzler R, Barkai E, Klafter J 1999 Deriving fractional Fokker-Planck equations from a generalised master equation.  _EPL_ 46 431–436
  * [76] Deng W, Hou R, Wang W, Xu P 2020 _Modeling Anomalous Diffusion: From Statistics to Mathematics_ (Singapore: World Scientific)
  * [77] Zhang X, Crawford J W, Deeks L K, Stutter M I, Bengough A G, Young I M 2005 A mass balance based numerical method for the fractional advection-dispersion equation: Theory and application.  _Water Resour. Res._ 41
  * [78] Metzler R, Rajyaguru A, Berkowitz B 2022 Modelling anomalous diffusion in semi-infinite disordered systems and porous media.  _New J. Phys._ 24 123004
  * [79] Berkowitz B, Kosakowski G, Margolin G, Scher H 2001 Application of continuous time random walk theory to tracer test measurements in fractured and heterogeneous porous media.  _Ground Water_ 39 593–604
  * [80] Berkowitz B, Cortis A, Dentz M, Scher H 2006 Modeling non-Fickian transport in geological formations as a continuous time random walk.  _Rev. Geophys._ 44
  * [81] Berkowitz B, Scher H 2009 Exploring the nature of non-Fickian transport in laboratory experiments.  _Adv. Water Resour._ 32 750–755
  * [82] Cortis A, Berkowitz B 2004 Anomalous transport in “Classical” soil and sand columns.  _Soil Sci. Soc. Am. J._ 68 1539–1548
  * [83] de Moraes R, Ozelim L, Cavalcante A 2022 Generalized skewed model for spatial-fractional advective–dispersive phenomena.  _Sustainability_ 14
  * [84] Doerries T J, Chechkin A V, Schumer R, Metzler R 2022 Rate equations, spatial moments, and concentration profiles for mobile-immobile models with power-law and mixed waiting time distributions.  _Phys. Rev. E_ 105 014105
  * [85] Redner S 2001 _A Guide to First-Passage Processes_ (Cambridge University Press, Cambridge)
  * [86] Metzler R, Chechkin A V, Klafter J 2012 _Lévy Statistics and Anomalous Transport: Lévy Flights and Subdiffusion_ , 1724–1745 (New York: Springer)
  * [87] Wardak A 2020 First passage leapovers of Lévy flights and the proper formulation of absorbing boundary conditions.  _J. Phys. A: Math. Theor._ 53 375001
  * [88] Zan W, Xu Y, Metzler R, Kurths J 2021 First-passage problem for stochastic differential equations with combined parametric gaussian and Lévy white noises via path integral method.  _J. Comput. Phys._ 435 110264
  * [89] Höll M, Nissan A, Berkowitz B, Barkai E ArXiv:2208.10262
  * [90] Chechkin A V, Metzler R, Gonchar V Y, Klafter J, Tanatarov L V 2003 First passage and arrival time densities for Lévy flights and the failure of the method of images.  _J. Phys. A: Math. Gen._ 36 L537
  * [91] Sokolov I M, Metzler R 2004 Non-uniqueness of the first passage time density of Lévy random processes.  _J. Phys. A: Math. Gen._ 37 L609–L615
  * [92] Palyulin V V, Blackburn G, Lomholt M A, Watkins N W, Metzler R, Klages R, Chechkin A V 2019 First passage and first hitting times of Lévy flights and Lévy walks.  _New J. Phys._ 21 103028
  * [93] Padash A, Chechkin A V, Dybiec B, Pavlyukevich I, Shokri B, Metzler R 2019 First-passage properties of asymmetric Lévy flights.  _J. Phys. A: Math. Theor._ 52 454004
  * [94] Barkai E, Garini Y, Metzler R 2012 Strange kinetics of single molecules in living cells.  _Phys. Today_ 65 29–35