# On relation between generalized diffusion equations and subordination
schemes

A. Chechkin<EMAIL_ADDRESS>Institute of Physics and Astronomy,
Potsdam University, Karl-Liebknecht-Strasse 24/25, 14476 Potsdam-Golm, Germany
Akhiezer Institute for Theoretical Physics, Akademicheskaya Str. 1, 61108
Kharkow, Ukraine  I.M. Sokolov<EMAIL_ADDRESS>Institut für
Physik and IRIS Adlershof, Humboldt Universität zu Berlin, Newtonstraße 15,
12489 Berlin, Germany

###### Abstract

Generalized (non-Markovian) diffusion equations with different memory kernels
and subordination schemes based on random time change in the Brownian
diffusion process are popular mathematical tools for description of a variety
of non-Fickian diffusion processes in physics, biology and earth sciences.
Some of such processes (notably, the fluid limits of continuous time random
walks) allow for either kind of description, but other ones do not. In the
present work we discuss the conditions under which a generalized diffusion
equation does correspond to a subordination scheme, and the conditions under
which a subordination scheme does possess the corresponding generalized
diffusion equation. Moreover, we discuss examples of random processes for
which only one, or both kinds of description are applicable.

##  I Introduction

In his seminal paper of 1961 Robert Zwanzig has introduced a generalized non-
Markovian Fokker-Planck equation Zwanzig  with a memory kernel, a GFPE in what
follows. The work was much cited, due to the fact that the approach to the
derivation of this equation based on the projection operator formalism has
found its application in the variety of problems in non-equilibrium
statistical physics Grabert . The equation itself was however hardly used,
except for obtaining Markovian approximations. The situation changed when the
generalized Fokker-Planck equations with power-law memory kernels gained
popularity in different fields. This kind of GFPEs, called fractional Fokker-
Planck equations (FFPEs) describe the continuous (long time-space) limit of
continuous time random walks (CTRW) with power-law waiting times Hilfer ;
Compte ; MeerScheff2019 , which gives a physical foundation and explains broad
applicability of FFPE. Such equations proved useful for description of
anomalous transport processes in different media. Applications range from
charge transport in amorphous semiconductors to underground water pollution
and motion of subcellular units in biology, see, e.g., the reviews
MetzlerKlafter ; MeKla-2 ; PhysToday  and references therein, as well as the
Chapters in collective monographs AnoTrans ; FracDyn . Further important
generalizations involve kernels consisting of mixtures of power laws, which
correspond to distributed-order fractional derivatives CheGorSok2002 ;
ChechkinFCAA ; CheKlaSok2003 ; APPB ; Naber2004 ; SoKla2005 ; SokChe2005 ;
UmaGor2005 ; MeerScheff2005 ; MeerScheff2006 ; Langlands ; Hanyga ; MaiPag2007
; MaiPaGo2007 ; CheGorSok2008 ; Kochubei2008 ; Meer2011 ; CheSoKla2012 , or
truncated (tempered) power laws SokCheKlaTruncated ; Stanislavsky1 ; MeerGRL ;
Baeumer . Other kernels in use include combinations of power laws with Mittag-
Leffler functions and with generalized Mittag-Leffler functions (Prabhakar
derivatives) TriChe2015 ; StanWer2016 ; Trifce1 ; Trifce2 ; StanWer2018 ;
Trifce3 ; StanWer2019 .

It is well-known that the Markovian, ”normal” Fokker-Planck equation, can be
obtained from the Langevin equation, the stochastic differential equation for
a Brownian motion under the action of an external force Chandra . Similarly,
the FFPE follows from two Langevin equations, giving a parametric
representation of the time dependence of the coordinate. These equations
describe the evolution of the coordinate and of the physical time in an
internal time-like variable (operational time) Fogedby ; Baule1 ; Baule2 ;
Kleinhans ; Hofmann ; Trifce4 . Such an approach is closely related to the
concept of subordination, i.e. random time change in a random process: A
process $X(\tau(t))$ is said to be subordinated to the process $X(\tau)$ under
operational time $\tau(t)$ being a random process with non-negative increments
Feller . The FFPEs discussed above thus describe the Brownian motion, possibly
in a force field, under a random time change, i.e. a process subordinated to a
Brownian motion with or without drift Meerschaert1 ; Stanislavsky2 ; Gorenflo1
; Gorenflo2 ; Gorenflo . Subordination schemes not only deliver a method of
analytical solution of FFPE SaiZasl ; Eli ; Meerschaert4 ; Meerschaert2  or
GFPE SokSub  by its integral transformation to the usual, Markovian,
counterpart, but also give the possibility of stochastic simulations of the
processes governed by GFPEs Kleinhans ; SaiUt1 ; MaiPa2003 ; Piryatinska ;
MaiPa2006 ; GoMai2007 ; Marcin1 ; Marcin2 ; Gajda ; Meerschaert3 ;
Stanislavsky3 ; Annunziato ; Marcin3 ; Marcin4 ; Stanislavsky4 ; Stanislavsky5
.

The subordination approach was also used in BNG ; Vittoria1 ; Vittoria2  in
describing, within the diffusive diffusivity model, a recently discovered but
widely spread phenomenon of Brownian yet non-Gaussian diffusion, a kind of
diffusion process in which the mean squared displacement (MSD) grows linearly
in time, like in normal, Fickian diffusion, but the probability density of the
particles’ displacements shows a (double-sided) exponential rather than
Gaussian distribution, at least at short or intermediate times Wang1 ; Wang2 ;
Chub ; Sebastian ; Cherail ; Grebenkov1 ; Grebenkov2 ; Korean ; Sandalo ;
RalfEPJB .

This broad use of generalized (not necessarily fractional) Fokker-Planck
equations on one hand, and of random processes subordinated to the Brownian
motion, on the other hand, urges us to put a question on the relation of these
two kinds of description. In other words, the following questions arise: (i)
given a GFPE (with a certain memory kernel), can one find the corresponding
subordinator or show that none exists, and (ii) given a subordination scheme,
can one find the corresponding GFPE (if any), or show that none exists. These
two questions are addressed in our paper.

Our main statements are summarized as follows. Not all valid GFPEs correspond
to subordination schemes. Not all subordination schemes possess a
corresponding GFPE. In our paper we give criteria to check, whether a
subordination scheme possesses a GFPE, and whether a particular GFPE
corresponds to a subordination scheme or not. We moreover discuss examples of
particular stochastic processes of interest by themselves, having one or
another description, or both.

##  II Generalized Fokker-Planck equations and subordination schemes

Let us first present the objects of our investigation: The GFPEs of a specific
form, and the two kinds of subordination schemes as they are discussed in the
literature cited above.

###  II.1 Generalized Fokker-Planck equations

In present paper we discuss equations of the form

$\frac{\partial}{\partial
t}P(\mathbf{x},t)=\hat{\Phi}\mathcal{L}P(\mathbf{x},t),$ (1)

with the linear integrodifferential operator $\hat{\Phi}$ acting on time
variable, and ${\cal L}$ is a time-independent linear operator acting on a
function of spatial variable(s) $\mathbf{x}$. We note that Eq.(1) is the most
popular, but not the most general form of such equations; in Ref. Zwanzig  a
more general form was derived. This includes the possible additional
coordinate or time dependence of $\hat{\Phi}$ and ${\cal L}$, respectively. In
the case of time dependence of ${\cal L}$ or of position-dependence of
$\hat{\Phi}$ the operators may not commute. Such situations were discussed
e.g. in Refs. SokKla  and Inhomogeneous , and are not a topic of present
investigation.

In the time domain the corresponding equations are always representable as a
GFPE

$\frac{\partial}{\partial
t}P(\mathbf{x},t)=\int_{0}^{t}\Phi(t-t^{\prime})\mathcal{L}P(\mathbf{x},t^{\prime})dt^{\prime}$
(2)

where the memory kernel $\Phi(t-t^{\prime})$ can be a generalized function
(contain delta-functions or derivatives thereof). Sometimes the corresponding
equations come in a form

$\frac{\partial}{\partial t}P(\mathbf{x},t)=\frac{\partial}{\partial
t}\int_{0}^{t}M_{R}(t-t^{\prime}){\cal L}P(\mathbf{x},t^{\prime})dt^{\prime},$
(3)

or

$\int_{0}^{t}M_{L}(t-t^{\prime})\frac{\partial}{\partial
t^{\prime}}P(\mathbf{x},t^{\prime})dt^{\prime}={\cal L}P(\mathbf{x},t),$ (4)

however, Eqs. (3) and (4) can be reduced to Eq.(2). Indeed, in the Laplace
domain Eq.(2) reads

$u\tilde{P}(\mathbf{x},u)-P(\mathbf{x},0)=\tilde{\Phi}(u)\mathcal{L}\tilde{P}(\mathbf{x},t)$

with

$\tilde{\Phi}(u)=\left\\{\begin{array}[]{ll}u\tilde{M}_{R}(u)&\mbox{for the
case Eq.(\ref{NMFPE1})}\\\ 1/\tilde{M}_{L}(u)&\mbox{for the case
Eq.(\ref{NMFPE2}),}\end{array}\right.$ (5)

where $\tilde{M}_{...}(u)=\int_{0}^{\infty}M_{...}(t)e^{-ut}dt$. We note that
for special cases of integral kernels corresponding to fractional or
distributed-order derivatives Eqs.(3) and (4) were called “modified” and
“normal” forms of generalized Fokker-Planck equation, respectively APPB .
Essentially, in the most cases the equation can be expressed in the either
form, with the left or the right memory kernel, but sometimes one of the forms
is preferable APPB ; SoKla2005 ; CheSoKla2012 ; Trifce4 . Finally, one can
also consider schemes with integral kernels on both sides, for which case

$\tilde{\Phi}(u)=u\frac{\tilde{M}_{R}(u)}{\tilde{M}_{L}(u)}.$

From now on we will use one-dimensional notation for $x$. Generalization to
higher spatial dimensions is straightforward. In all our examples we will
concentrate mostly on the case of free diffusion

$\mathcal{L}=D\frac{\partial^{2}}{\partial x^{2}}.$ (6)

The coefficient $D$ has a dimension of the normal diffusion coefficient,
$[D]=[\mathrm{L}^{2}/\mathrm{T}]$, the operator $\hat{\Phi}$ is therefore
dimensionless, and its integral kernel has a dimension of the inverse time. We
note however that concentration on free diffusion is not a restriction for the
generality of our approach since it is only about temporal operators and
temporal parts of the subordination procedures.

###  II.2 Kernel of GFPE uniquely defines MSD in free diffusion

Here we show that the form of the memory kernel uniquely defines the mean
squared displacement (MSD) in free diffusion, and vice versa (under mild
restrictions). Let us consider the MSD in free diffusion (i.e. without
external force and in absence of external boundaries). We multiply both parts
of Eq.(1) by $x^{2}$ and integrate over $x$ to get

$\frac{d}{dt}\int_{-\infty}^{\infty}x^{2}P(x,t)dx=D\hat{\Phi}\int_{-\infty}^{\infty}x^{2}\frac{\partial^{2}}{\partial
x^{2}}P(x,t)dx.$ (7)

Integrating the right hand side of Eq.(7) by parts twice and assuming the PDF
$P(x,t)$ to vanish at infinity together with its first derivative we get for
the r.h.s.

$\int_{-\infty}^{\infty}x^{2}\frac{\partial^{2}}{\partial x^{2}}P(x,t)dx=2,$

so that the evolution of the MSD is governed by

$\frac{d}{dt}\langle x^{2}(t)\rangle=2D\hat{\Phi}1,$ (8)

with operator $\hat{\Phi}$ acting on a numeric constant. Passing to the
Laplace representation we obtain

$u\langle x^{2}(u)\rangle=2D\tilde{\Phi}(u)\frac{1}{u}$

where we assumed to start from an initial condition concentrated at the
origin, $\langle x^{2}(t=0)\rangle=0$. This uniquely defines $\tilde{\Phi}(u)$
via the MSD:

$\tilde{\Phi}(u)=\frac{1}{2D}u^{2}\langle x^{2}(u)\rangle.$ (9)

Let us consider our first example. If the MSD in free motion grows linearly in
time, $\langle x^{2}(t)\rangle=2Dt$, we have $\langle
x^{2}(u)\rangle=2D/u^{2}$ and therefore $\tilde{\Phi}(u)\equiv 1$ (a unit
operator, an integral operator with a $\delta$-functional kernel). This means
that the only GFPE leading to the linear growth of the MSD is a usual,
Fickian, diffusion equation, for which the PDF is Gaussian at all times.
Therefore, the GFPE is not a valid instrument to describe the BnG diffusion,
contrary to what is claimed in Korean . On the other hand, the subordination
schemes of Chub ; BNG ; Grebenkov1  do describe the phenomenon.

Now we can turn to the main topic of our work: which processes can and which
can not be described by the GFPEs. To this end we first discuss a specific
integral representation of the solution to a GFPE and its relation to
subordination schemes.

###  II.3 Subordination schemes

Let $X(\tau)$ be a random process parametrized by a “time-like” variable
$\tau$. Let $\tau$ by itself be a random process, parametrized by the physical
time, or clock time, $t$. The random process $X(\tau)$ is called the parent
process, the random variable $\tau$ is called the operational time, and the
random process $\tau(t)$ the directing process, or subordinator. The process
$t(\tau)$ is called the leading process of the subordination scheme, see
Gorenflo  for the consistent explanation of the terminology used. The process
$X(\tau(t))$ is said to be subordinated to $X(\tau)$ under the operational
time $\tau$. The properties of $X(\tau)$ and $\tau(t)$ (or alternatively,
$t(\tau)$) fully define the properties of the composite process $X(\tau(t))$.

The fact that the variable $\tau$ is time-like means that the directing
process $\tau(t)$ preserves the causality: from $t_{2}>t_{1}$ it must follow
that $\tau(t_{2})\geq\tau(t_{1})$: the directing process of a subordination
scheme is increasing at least in a weak sense, that is possesses non-negative
increments. It is moreover assumed that $\tau(0)=0$: the count of operational
time starts together with switching the physical clock.

The composite random function $X(\tau(t))$ can be defined in two ways. The
first way corresponds to an explicit definition by defining the (stochastic)
equations governing $X(\tau)$ and $\tau(t)$. The second way defines the
function parametrically, so that the equations for $X(\tau)$ and $t(\tau)$ are
given.

A process subordinated to a Brownian motion with drift is a process whose
parent process is defined by a stochastic differential equation

$\frac{d}{d\tau}x(\tau)=F(x(\tau))+\sqrt{2D}\xi(\tau)$ (10)

with white Gaussian noise $\xi(\tau)$, with $\langle\xi(\tau)\rangle=0$,
$\langle\xi(\tau_{1})\xi(\tau_{1})\rangle=\delta(\tau_{1}-\tau_{2})$, whose
strength is given by a diffusion coefficient $D$, and with deterministic drift
$F(x(\tau))$ which will be considered absent in our examples concentrating on
free diffusion.

As example of the explicit, or direct, subordination scheme we name the
minimal diffusing diffusivity model,

$\displaystyle\frac{dx(\tau)}{d\tau}$ $\displaystyle=$
$\displaystyle\sqrt{2}\xi(\tau),$ $\displaystyle\frac{d\tau}{dt}$
$\displaystyle=$ $\displaystyle D(t),$

where the random diffusion coefficient $D(t)$ is a squared Ornstein-Uhlenbeck
process BNG .

In the parametric subordination scheme the dependence of $t(\tau)$ is given,
again in a form of an SDE, or via an additional transform of its solution. The
classical Fogedby scheme Fogedby  corresponds to a stochastic differential
equation

$\frac{dt}{d\tau}=\lambda(\tau)$ (11)

with $\lambda(t)$ being a one-sided Lévy noise. The clock time given by the
solution of this equation is
$t(\tau)=\int_{0}^{\tau}\lambda(\tau^{\prime})d\tau^{\prime}$, , and the PDF
$q(t,\tau)$ of the process $t(\tau)$ is given by a one-sided $\alpha$-stable
Lévy law, such that its Laplace transform in $t$ variable reads
$\tilde{q}(u,\tau)=\exp(-u^{\alpha}\tau),0<\alpha<1$. This scheme corresponds
to a diffusive limit of CTRW with a power law waiting time distribution.

The correlated CTRW model of Ref. Hofmann  is another example of the
parametric scheme where $t(\tau)$ is obtained by an additional integration of
$\lambda(\tau)$ above:

$\frac{dt}{d\tau}=\int_{0}^{\tau}\Psi(\tau-\tau^{\prime})\lambda(\tau^{\prime})d\tau^{\prime}.$
(12)

The previous case is restored if the kernel $\Psi(\tau)$ is a
$\delta$-function.

The attractiveness of subordination schemes lays in the fact that if the
solution for the PDFs of the parent process at a given operational time, and
of the directing process at a given physical time are known, the PDF $P(x,t)$
of the subordinated process can be obtained simply by applying the Bayes
formula. Let $f(x,\tau)$ be the PDF of $x=X(\tau)$ for a given value of the
operational time $\tau$, and $p(\tau,t)$ the PDF of the operational time for
the given physical time $t$. Than the PDF of $x=X(t)=X(\tau(t))$ is given by

$P(x,t)=\int_{0}^{\infty}f(x,\tau)p(\tau,t)d\tau,$ (13)

which in this context is called the integral formula of subordination Feller .
The PDF $p(\tau,t)$ of the operational time at a given clock time is delivered
immediately by explicit schemes, and can be obtained for parametric
subordination schemes by using an additional transformation Baule1 ; Gorenflo
, see below. Note that the PDF $f(x,\tau)$ for a process subordinated to a
Brownian motion always satisfies a usual, Markovian Fokker-Planck equation

$\frac{\partial}{\partial\tau}f(x,\tau)={\cal L}f(x,\tau),$ (14)

with $\mathcal{L}f(x,\tau)=-\frac{\partial}{\partial
x}F(x)f(x,\tau)+D\frac{\partial^{2}}{\partial x^{2}}f(x,\tau)$ by virtue of
Eq.(10) of the subordination schemes.

##  III The sufficient condition for GFPE to have a subordination scheme

In Ref. SokSub  it was shown that the formal solution of the GFPE (2) can be
obtained in a form of an integral decomposition

$P(x,t)=\int_{0}^{\infty}f(x,\tau)T(\tau,t)d\tau,$ (15)

where $f(x,\tau)$ is a solution of a Markovian FPE with the same Fokker-Planck
operator ${\cal L}$, Eq.(14), and for the same initial and boundary
conditions. Here the function $T(\tau,t)$ is normalized in its first variable
and connected with the memory kernel of the GFPE, as it is discussed in this
Section below. The corresponding form of solution was obtained in SaiZasl ;
Eli  for the fractional diffusion and Fokker-Planck equations, and was applied
to the fractional Kramers equation in BarSil ; its more general discussion
followed in Sok1 . Equation (15) is akin to the integral formula of
subordination, Eq.(13). However, the PDF $P(x,t)$ in Eq. (15) may or may not
correspond to a PDF of a random process subordinated to the Brownian motion
with drift (as described by the ordinary FPE), since $T(\tau,t)$ may or may
not be a conditional probability density of $\tau$ at time $t$, e.g.
$T(\tau,t)$ may get negative. In Ref. SokSub  the kernels corresponding to
subordination schemes with non-negative $T(\tau,t)$ were called ”safe”, while
the kernels not corresponding to any subordination scheme, for which
$T(\tau,t)$ oscillate, were called “dangerous”. For safe kernels the non-
negativity of solutions of GFPE corresponding to non-negative initial
conditions is guaranteed by virtue of Eq.(13). The sufficient condition for
Eq. (15) to correspond to a subordination scheme will be considered later.

If one assumes that the solution of GFPE (2) can be obtained in the form of
integral decomposition (15) and then insert such form in Eq. (2), one gets the
Laplace transform of the function $\tilde{T}(\tau,t)$ in its second variable,
$\tilde{T}(\tau,u)=\int_{0}^{\infty}T(\tau,t)e^{-ut}dt$, as SokSub

$\tilde{T}(\tau,u)=\frac{1}{\tilde{\Phi}(u)}\exp\left[-\tau\frac{u}{\tilde{\Phi}(u)}\right].$
(16)

This however, does not answer the questions what are the conditions under
which such a solution holds and whether it is unique. Below we discuss these
issues in some detail by presenting an alternative derivation, i.e. explicitly
constructing the solution.

Let us start from our Eq.(2) and integrate its both parts over time getting

$\displaystyle\int_{0}^{t}\frac{\partial}{\partial
t^{\prime\prime}}P(x,t^{\prime\prime})dt^{\prime\prime}=P(x,t)-P(x,0)$
$\displaystyle\qquad=\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}\Phi(t^{\prime\prime}-t^{\prime}){\cal
L}P(x,t^{\prime})dt^{\prime}.$

Now we exchange the sequence of integrations in $t^{\prime}$ and
$t^{\prime\prime}$ on the r.h.s.,

$\displaystyle\int_{0}^{t}dt^{\prime\prime}\int_{0}^{t^{\prime\prime}}K(t^{\prime\prime}-t^{\prime}){\cal
L}P(x,t^{\prime})dt^{\prime}$
$\displaystyle\qquad=\int_{0}^{t}dt^{\prime}{\cal
L}P(x,t^{\prime})\int_{t^{\prime}}^{t}\Phi(t^{\prime\prime}-t^{\prime})dt^{\prime\prime},$

getting the integral form

$P(x,t)-P(x,0)=\int_{0}^{t}K(t-t^{\prime}){\cal L}P(x,t^{\prime})dt^{\prime},$
(17)

with the integral kernel
$K(t)=\int_{0}^{t}\Phi(t^{\prime\prime})dt^{\prime\prime}$ whose Laplace
transform is equal to $\tilde{\Phi}(u)/u$. Using the condition $\tau(t=0)=0$
we substitute the assumed solution form, Eq.(15), into Eq. (17):

$\displaystyle\int_{0}^{\infty}f(x,\tau)T(\tau,t)d\tau-P(x,0)=$
$\displaystyle\qquad\int_{0}^{t}dt^{\prime}K(t-t^{\prime})\int_{0}^{\infty}d\tau\mathcal{L}f(x,\tau)T(\tau,t^{\prime}).$

Now we use the assumption that $f(x,\tau)$ is the solution of a Markovian
Fokker-Planck equation, and make the substitution
$\mathcal{L}f(x,\tau)=\frac{\partial}{\partial\tau}f(x,\tau)$. We get:

$\displaystyle\int_{0}^{\infty}f(x,\tau)T(\tau,t)d\tau-P(x,0)=$
$\displaystyle\qquad\int_{0}^{t}dt^{\prime}K(t-t^{\prime})\int_{0}^{\infty}d\tau\left(\frac{\partial}{\partial\tau}f(x,\tau)\right)T(\tau,t^{\prime}).$

Performing partial integration in the inner integral on the r.h.s., and
interchanging the sequence of integrations in $t^{\prime}$ and in $\tau$ in
the integral which appears in the r.h.s. we arrive at the final expression

$\displaystyle\int_{0}^{\infty}f(x,\tau)T(\tau,t)d\tau-P(x,0)=$
$\displaystyle\qquad\int_{0}^{t}K(t-t^{\prime})\left[f(x,\infty)T(\infty,t^{\prime})-f(x,0)T(0,t^{\prime})\right]dt^{\prime}$
$\displaystyle\qquad-\int_{0}^{\infty}d\tau
f(x,\tau)\frac{\partial}{\partial\tau}\int_{0}^{t}K(t-t^{\prime})T(\tau,t^{\prime})dt^{\prime}.$

Now we request that the l.h.s. and the r.h.s. are equal at any time $t$ for
all admissible functions $f(x,\tau)$ satisfying the Fokker-Planck equation
$\frac{\partial}{\partial\tau}f(x,\tau)=\mathcal{L}f(x,\tau)$ irrespective of
the particular form of the linear operator $\mathcal{L}$ and of the boundary
and initial conditions. This gives us three conditions:

$\displaystyle\int_{0}^{\infty}f(x|\tau)T(\tau,t)d\tau=$
$\displaystyle\qquad-\int_{0}^{\infty}d\tau
f(x|\tau)\frac{\partial}{\partial\tau}\int_{0}^{t}K(t-t^{\prime})T(\tau,t^{\prime}),$
$\displaystyle-P(x,0)=-f(x,0)\int_{0}^{t}K(t-t^{\prime})T(0,t^{\prime}),$
$\displaystyle 0=\int_{0}^{t}K(t-t^{\prime})f(x,\infty)T(\infty,t^{\prime}),$

which can be rewritten as conditions on $T(\tau,t)$ only:

$\displaystyle
T(\tau,t)=-\frac{\partial}{\partial\tau}\int_{0}^{t}K(t-t^{\prime})T(\tau,t^{\prime})dt^{\prime},$
(18) $\displaystyle\int_{0}^{t}K(t-t^{\prime})T(0,t^{\prime})dt^{\prime}=1,$
(19)
$\displaystyle\int_{0}^{t}K(t-t^{\prime})T(\infty,t^{\prime})dt^{\prime}=0.$
(20)

In the Laplace domain Eq.(18) turns to a simple linear ODE

$-\frac{\tilde{\Phi}(u)}{u}\frac{\partial}{\partial\tau}\tilde{T}(\tau,u)=\tilde{T}(\tau,u)$

whose general solution is

$\tilde{T}(\tau,u)=C\cdot\exp\left(-\tau\frac{u}{\tilde{\Phi}(u)}\right)$

with the integration constant $C$. This integration constant is set by the
second equation, Eq.(19), which in the Laplace domain reads

$\frac{\tilde{\Phi}(u)}{u}\tilde{T}(0,u)=\frac{1}{u},$

so that $C=1/\Phi(u)$, and therefore

$\tilde{T}(\tau,u)=\frac{1}{\tilde{\Phi}(u)}\exp\left(-\tau\frac{u}{\tilde{\Phi}(u)}\right),$

which is our Eq.(16). The function $T(\tau,t)$ is normalized in its first
variable, which follows by the direct integration of Eq.(16):

$\int_{0}^{\infty}\tilde{T}(\tau,u)d\tau=\frac{1}{u},$ (21)

so that its inverse Laplace transform to the time domain is unity:

$\int_{0}^{\infty}T(\tau,t)d\tau=1$ (22)

for any $t>0$.

The third condition, Eq.(20), is fulfilled automatically provided
$T(\infty,t^{\prime})=0$ which implies $\tilde{T}(\infty,u)=0$. This is e.g.
always the case for non-negative kernels $\Phi(t)$ (as encountered in all our
examples) whose Laplace transform $\tilde{\Phi}(u)$ is positive for all $u$.
For non-positive kernels the property has to be checked explicitly.

Let us stress again that the solution in form of Eq.(16), which, as we have
seen, is applicable for a wide range of memory kernels $\Phi$, may or may not
correspond to some subordination scheme. We note, however, the fact that the
kernel does not correspond to a subordination scheme does not devaluate the
corresponding GFPE by itself, and does not mean that this leads to negative
probability densities. As an example of a “dangerous” kernel let us consider a
simple exponential kernel, $\Phi(t)=re^{-rt}$ (where the prefactor $r$ of the
exponential is added to keep the correct dimension of $\Phi$). For example, a
generalized diffusion equation with an exponential kernel,

$\frac{\partial}{\partial
t}p(x,t)=\int_{0}^{t}re^{-r(t-t^{\prime})}D\frac{\partial^{2}}{\partial
x^{2}}p(x,t^{\prime})dt^{\prime},$ (23)

is essentially the Cattaneo equation for the diffusion with finite propagation
speed, i.e. a kind of a telegrapher’s equation, as can be seen by taking a
derivative of its both sides w.r.t. $t$:

$\frac{\partial^{2}}{\partial t^{2}}p(x,t)=r\frac{\partial}{\partial
t}p(x,t)+rD\frac{\partial^{2}}{\partial x^{2}}p(x,t).$ (24)

The solutions to Eq.(24) for non-negative initial conditions are known to be
non-negative on the whole real line, but changing the operator ${\cal L}$ from
a diffusion to a more general one (e.g. to diffusion in presence of the
constant force) may lead to oscillating solutions SokSub .

The reason for this, within our line of argumentation, is that the function
$\tilde{T}(\tau,u)$ for equation (23),

$\tilde{T}(\tau,u)=(1+ur^{-1})\exp[-\tau u(1+r^{-1}u)]$ (25)

is not a Laplace transform of a non-negative function. We remind that the
function $\tilde{\phi}(u),0\leq u\leq\infty$, is a Laplace transform of a non-
negative function $\phi(t)$ defined on the non-negative half-axis, if and only
if $\tilde{\phi}(u)$ is completely monotone, i.e. its derivatives satisfy

$(-1)^{n}\phi^{(n)}(u)\geq 0$

for $n=0,1,2,...$ Feller . On the other hand, it is easy to see that the
second derivative of $\tilde{T}(\tau,u)$ changes its sign. Moreover, using the
mean value theorem, it is not hard to show that the Laplace transform of any
non-negative function integrable to unity cannot decay for $u\to\infty$ faster
than exponentially (see Appendix A), which is not true for Eq.(25).

Summarizing the result of this Section, we see that not all GFPEs can
correspond to subordination schemes. The sufficient condition to have such
scheme is the following: The kernel $\Phi(t)$ in the GFPE (2) is such that the
function $\tilde{T}(\tau,u)$ given by Eq.(16) is completely monotone as a
function of $u$. This always corresponds to subordination scheme, for which
$T(\tau,t)$ can be interpreted as a probability density function of the
operational time $\tau$ for given physical time $t$, $T(\tau,t)\equiv
p(\tau,t)$.

##  IV What subordination schemes do have a GFPE?

Let us first perform some simple manipulations while assuming that the
function $T(\tau,t)$ in the integral decomposition formula (15) does
correspond to a subordination scheme, and thus has a meaning of the PDF of
operational time $\tau$ for a given physical time $t$, $T(\tau,t)\equiv
p(\tau,t)$. Then, in the Laplace domain the function $\tilde{p}(\tau,u)$,
Eq.(16), can be represented as

$\tilde{p}(\tau,u)=-\frac{d}{d\tau}u^{-1}\exp\left[-\tau\frac{u}{\tilde{\Phi}(u)}\right],$
(26)

so that in the time domain we have

$p(\tau,t)=-\frac{d}{d\tau}\int_{0}^{t}q(t^{\prime},\tau)dt^{\prime},$

where the function $q(t,\tau)$ is given by the inverse Laplace transform of

$\tilde{q}(u,\tau)=\exp\left[-\tau\frac{u}{\tilde{\Phi}(u)}\right].$ (27)

Thus

$\int_{\tau_{0}}^{\infty}p(\tau,t_{0})d\tau=\int_{0}^{t_{0}}q(t,\tau_{0})dt.$
(28)

Now we proceed to show that, as discussed already in Baule1 ; Gorenflo , the
function $q(t,\tau)$ has a clear physical meaning: this is namely the PDF of
clock times corresponding to the given operational time $\tau$. We note here
that in spite of the fact that Eq.(28) was obtained for a specific form of the
PDF $p(\tau,t)$ given by Eq.(26), it is more general and gives the relation
between the PDFs of a (weakly) increasing process and its inverse.

Indeed, let us consider a set of monotonically non-decaying functions, either
continuous (diffusive limit) or of càdlàg type (genuine continuous time random
walks) on a $(t,\tau)$ plane, see Fig. 1. The integral on the l.h.s. of
Eq.(28) counts all functions (with their probability weights) which cross the
horizontal segment $t\in[0,t_{0}),\tau=\tau_{0}$, the integral on the r.h.s.
counts all functions crossing the semi-infinite vertical segment
$t=t_{0},\tau\in[\tau_{0},\infty)$. The set of these functions is the same:
any monotonically non-decaying function crossing the horizontal segment or
passing from its one side to another side on a jump has to cross the vertical
one. No non-decaying function which never crossed the horizontal segment can
cross the vertical one. Therefore such a monotonicity implies that

$\textrm{Prob}(\tau>\tau_{0}|t_{0})=\textrm{Prob}(t<t_{0}|\tau_{0}),$

where the probabilities are defined on the set of the corresponding
trajectories. The physical meaning of the functions
$\textrm{Prob}(\tau>\tau(t))$ and $\textrm{Prob}(t<t(\tau))$ is that they
represent the survival probability and the cumulative distribution function
for the operational time and for the clock time, respectively. In the
continuous case the PDF of a clock time given operational time is then given
by

$\displaystyle q(t,\tau)$ $\displaystyle=$
$\displaystyle\frac{d}{dt}\textrm{Prob}(t<t(\tau))=\frac{d}{dt}\textrm{Prob}(\tau>\tau(t))$
(29) $\displaystyle=$
$\displaystyle\frac{d}{dt}\int_{\tau}^{\infty}p(\tau^{\prime},t)d\tau^{\prime}.$

This statement allows for immediate transition between direct (random variable
change $\tau(t)$) and parametric (inverse) subordination schemes. This also
gives a necessary condition for a process obeying GFPE to be a subordinated
one.

Figure 1: A schematic picture explaining the nature of Eq.(28). Black lines:
the case of continuous time traces. Here all monotonically non-decaying traces
crossing the horizontal segment $[0,\tau_{0})$ (shown in blue online) also
cross the vertical one $[\tau_{0},\infty)$ (shown in red online), and
therefore contribute equally to the r.h.s. and to the l.h.s. of Eq.(28). A
gray dashed line shows exemplary a càdlàg piecewise constant function with
jumps, the genuine operational time of a CTRW. Any càdlàg trace passing at a
jump from below to above the horizontal segment has to cross the vertical
segment during the waiting time.

To obtain such a necessary condition, let us fix two subsequent non-
intersecting operational time intervals $\tau_{1}$ and $\tau_{2}$
corresponding to the physical time intervals $t_{1}$ and $t_{2}$. Than
$t(\tau_{1}+\tau_{2})=t(\tau_{1})+t(\tau_{2})=t_{1}+t_{2}$. Than it follows
from Eq.(27) that

$\tilde{q}(u,\tau_{1}+\tau_{2})=\exp\left[-\tau_{1}\frac{u}{\tilde{\Phi}(u)}\right]\cdot\exp\left[-\tau_{2}\frac{u}{\tilde{\Phi}(u)}\right],$

or, denoting the Laplace characteristic functions of $t_{1}$ and $t_{2}$ by
$\tilde{\theta}_{t_{1}}(u)$ and $\tilde{\theta}_{t_{2}}(u)$, and the
characteristic function of their sum by $\tilde{\theta}_{t_{1}+t_{2}}(u)$ we
get

$\tilde{\theta}_{t_{1}+t_{2}}(u)=\tilde{\theta}_{t_{1}}(u)\cdot\tilde{\theta}_{t_{2}}(u),$
(30)

that is the random variables $t_{1}$ and $t_{2}$ are sub-independent Hamedani
. This property puts a necessary condition for the possibility to describe the
subordination scheme by a generalized FPE. The condition is always fulfilled
when $t_{1}$ and $t_{2}$ are independent (e.g. for a parametric Fogedby scheme
Fogedby ).

Therefore, the following statement can be made: Only subordination schemes in
which the increments of physical time $t(\tau)$ are sub-independent can be
described by GFPEs.

Combining Eqs.(27) and (30) we find the final criteria for the existence of
GFPE for a given subordination scheme.

  * •

A direct subordination scheme does posses a GFPE only if the Laplace transform
of the PDF $p(\tau,t)$ in its $t$ variable has the form

$\tilde{p}(\tau,u)=\frac{f(u)}{u}\exp(-\tau f(u)).$ (31)

The kernel of the ensuing GFPE is then $\tilde{\Phi}(u)=u/f(u)$. Here it might
be easier to check that the double Laplace transform has a form

$\tilde{\tilde{p}}(s,u)=\int_{0}^{\infty}d\tau\int_{0}^{\infty}dte^{-s\tau}e^{-ut}p(\tau,t)=\frac{f(u)}{u[s+f(u)]},$
(32)

i.e. the function $F(s,u)=1/\tilde{\tilde{p}}(s,u)$ is a linear function in
$s$: $F(s,u)=a(u)s+1$. Using this criterion one can show that the BnG model of
Ref. BNG  does not posses a corresponding GFPE, see Appendix 2.

  * •

A parametric subordination scheme does posses a GFPE only if the
characteristic function (Laplace transform) of the PDF $q(t,\tau)$ in its $t$
variable has a form

$\tilde{q}(u,\tau)=\exp(-\tau f(u)),$ (33)

i.e. the $\tau$-dependence of this function must be simple exponential. The
$t$-variables corresponding to non-intersecting $\tau$ intervals are thus sub-
independent. The kernel of the ensuing GFPE is then $\tilde{\Phi}(u)=u/f(u)$.

Using the last criterion one can easily show that there is no GFPE
corresponding to correlated CTRW of Ref. Hofmann . The distribution of $t$ as
a function of $\tau$ in this model has a Laplace characteristic function

$\tilde{q}(u,\tau)=\exp\left[-u^{\alpha}\phi(\tau)\right]$ (34)

with

$\phi(\tau)=\int_{0}^{\tau}d\tau^{\prime}\left[\int_{\tau^{\prime}}^{\tau}d\tau^{\prime\prime}\Psi(\tau^{\prime\prime}-\tau^{\prime})\right]^{\alpha}$

where $\Psi(\tau)$ denotes the memory function for waiting times along the
trajectory expressed as a function of the number of steps, cf. Eq.(26) of Ref.
Hofmann . To correspond to any GFPE the argument of the exponential in Eq.(34)
has to be linear in $\tau$, i.e. $\phi(\tau)=a\tau$, where $a$=const. This
means that the square bracket in the expression for $\phi$ must be a constant
(equal to $a$) and therefore $\Psi(\tau)=a^{1/\alpha}\delta(\tau)$, which
corresponds to a standard, non-correlated CTRW.

##  V Conclusions

Growing awareness of the complexity of non-Markovian diffusion processes in
physics, earth sciences and biology gave rise to a spark of interest to
mathematical tools capable to describe such non-standard diffusion processes
beyond the Fick’s law. Generalized diffusion equations on one hand, and
subordination schemes, on the other hand, are the two classes of such
instruments, which were successfully used for investigation of a broad variety
of anomalous diffusion processes. For several situations, notably, for
decoupled continuous time random walks, both are applicable, and stand in the
same relation to each other as the Fokker-Planck equation and the Langevin
equation do for the case of normal, Fickian diffusion. In the present work we
address the question, whether this is always the case. The answer to this
question is negative: some processes described by the generalized diffusion
equations do not possess an underlying subordination scheme, i.e. cannot be
described by a random time change in the normal diffusion process. On the
other hand, many subordination schemes do not possess the corresponding
generalized diffusion equation. The example for the first situation is the
Cattaneo equation, which can be represented as a generalized diffusion with
exponential memory kernel, for which no subordination scheme exists. The
example of the second situation is the minimal model of Brownian yet non-
Gaussian diffusion and correlated CTRW, the subordination schemes for which we
show that no corresponding generalized diffusion equation can be put down. We
discuss the conditions under which one or the other description is applicable,
i.e. what are the properties of the memory kernel of the diffusion equation
sufficient for its relation with subordination, and what are the properties of
the random time change in the subordination scheme necessary for existence of
the corresponding generalized diffusion equation.

##  VI Acknowledgements

AC acknowledges support from the NAWA Project PPN/ULM/2020/1/00236.

##  Appendix A Laplace transform of a non-negative function

Let us show that the Laplace transform $\widetilde{\phi}(u)$ of any non-
negative function $\phi(t)$ integrable to a constant, i.e. with $\phi(t)\geq
0$ and with $0<\int_{0}^{\infty}\phi(t)dt<\infty$, cannot decay faster than
exponentially with the Laplace variable $u$:

$\widetilde{\phi}(u)=\int_{0}^{\infty}\phi(t)e^{-ut}dt\geq Ae^{-Bu},$

with a positive constant $A$ and a non-negative constant $B$.

To see this we use the following chain of relations:

$\displaystyle\int_{0}^{\infty}\phi(t)e^{-ut}dt$ $\displaystyle\geq$
$\displaystyle\int_{0}^{C}\phi(t)e^{-ut}dt$ $\displaystyle=$ $\displaystyle
e^{-Bu}\int_{0}^{C}\phi(t)dt=Ae^{-Bu}.$

Here $C$ is some cut-off value which is chosen such that
$\int_{0}^{C}\phi(t)dt>0$, which is always possible due to our assumptions
about the integrability of $\phi(t)$ to a positive constant. The inequality
follows from the mean value theorem for integrals: by assumptions $\phi(t)$ is
non-negative and integralble, and $e^{-ut}$ is, evidently, continuous. The
value of $B$ then follows the inequality $0\leq B\leq C$, i.e. is non-
negative.

For our function $T(\tau,t)$, Eq.(25), the Laplace transform $T(\tau,u)$ is
defined for $u=0$ so that $\int_{0}^{\infty}T(\tau,t)dt=1>0$ for any $\tau$.
Eq.(25) essentially corresponds to a Laplace transform of a function strongly
oscillating at small $t$.

##  Appendix B The minimal model of BNG diffusion

In the BnG model of Ref. BNG  the PDF $p(\tau,t)$ (denoted there as
$T(\tau,t))$ is defined via its Laplace transform in $\tau$ variable,

$\displaystyle\tilde{p}(s,t)=$
$\displaystyle\frac{e^{t/2}}{\sqrt{\frac{1}{2}\left(\sqrt{1+2s}+\frac{1}{\sqrt{1+2s}}\right)\mathrm{sinh}(t\sqrt{1+2s})+\mathrm{cosh}(t\sqrt{1+2s})}}.$

Let us take the double Laplace transform

$\tilde{\tilde{p}}(s,u)=\int_{0}^{\infty}\tilde{p}(s,t)e^{-ut}dt$ (35)

and check whether it has the form of Eq.(32).

We first denote $\alpha=\sqrt{1+2s}>1$ and rewrite the function
$\tilde{p}(s,t)$ as

$\displaystyle\tilde{p}(s,t)$ $\displaystyle=$
$\displaystyle\frac{e^{t/2}}{\sqrt{\frac{1}{2}\left(\alpha+\frac{1}{\alpha}\right)\mathrm{sinh}(\alpha
t)+\mathrm{cosh}(\alpha t)}}$ $\displaystyle=$
$\displaystyle\frac{2\sqrt{\alpha}e^{-\frac{\alpha-1}{2}t}}{(\alpha+1)^{2}}\left[1+\left(\frac{\alpha-1}{\alpha+1}\right)^{2}e^{-2\alpha
t}\right]^{-\frac{1}{2}}.$

Denoting $A=2\sqrt{\alpha}/(\alpha+1)^{2}$ and
$\zeta=\left[(\alpha-1)/(\alpha+1)\right]^{2}$ we get

$\tilde{p}(s,t)=Ae^{-\frac{\alpha-1}{2}t}\left(1+\zeta e^{-2\alpha
t}\right)^{-\frac{1}{2}}.$

Substituting this expression into Eq. (35) and denoting
$\beta=(\alpha-1+2u)/2$, we obtain

$\tilde{\tilde{p}}(s,u)=A\int_{0}^{\infty}\left(1+\zeta e^{-2\alpha
t}\right)^{-\frac{1}{2}}e^{-\beta t}dt.$

Now we change the variable of integration to $x=e^{-2\alpha t}$ to arrive to
the expression

$\tilde{\tilde{p}}(s,u)=\frac{A}{2\alpha}\int_{0}^{1}\left(1+\zeta
x\right)^{-\frac{1}{2}}x^{\frac{\beta}{2\alpha}-1}dx.$

This can be compared with the integral representation of the hypergeometric
function:

$\;{}_{2}F_{1}(a,b;c;z)=$
$\displaystyle\qquad\frac{\Gamma(c)}{\Gamma(b)\Gamma(c-b)}\int_{0}^{1}x^{b-1}(1-x)^{c-b-1}(1-xz)^{-a}dx,$

Eq.(15.3.1) of Ref. AbraSteg , from which we get $a=\frac{1}{2}$,
$b=\frac{\beta}{2\alpha}$, $c=b+1$, and $z=-\zeta$ so that,

$\tilde{\tilde{p}}(s,u)=\frac{A}{\beta}\;_{2}F_{1}\left(\frac{1}{2},\frac{\beta}{2a},1+\frac{\beta}{2a},-\zeta\right).$

Substituting the values of parameters we obtain:

$\displaystyle\tilde{\tilde{p}}(s,u)=\frac{4(1+2s)^{\frac{1}{4}}}{(1+\sqrt{1+2s})^{2}(\sqrt{1+2s}+2u-1)}\times$
$\qquad{}_{2}F_{1}\left[\frac{1}{2},\frac{\sqrt{1+2s}-1+2u}{4\sqrt{1+2s}},\frac{5\sqrt{1+2s}-1+2u}{4\sqrt{1+2s}},\right.$
$\displaystyle\qquad\qquad\left.-\left(\frac{\sqrt{1+2s}-1}{\sqrt{1+2s}+1}\right)^{2}\right].$

The function $F(s,u)=1/\tilde{\tilde{p}}(s,u)$ is not a linear function of $s$
for fixed $u$. This can be clearly seen when plotting the $s$-derivative of
this function for fixed $u$ with the help of Mathematica, see Fig. 2.

Figure 2: The $s$-derivative of the function $F(s,u)$ for $u=0.5,1$ and $2$,
shown by solid, dashed and dotted lines, respectively.

## References

  * (1) R. Zwanzig, Memory Effects in Irreversible Thermodynamics, Phys. Rev. 124, 983 (1961).
  * (2) H. Grabert, Projection operator techniques in nonequilibrium statistical mechanics, Springer, Berlin, 1982.
  * (3) R. Hilfer and L. Anton, Fractional master equations and fractal time random walks, Phys. Rev. E 51, R848 (1995).
  * (4) A. Compte, Stochastic foundations of fractional dynamics, Phys. Rev. E 53 (4), 4191 (1996).
  * (5) M. Meerschaert and H.P. Scheffler, Continuous time random walks and space-time fractional differential equations. In: A. Kochubei, Yu. Luchko (eds.) Handbook of Fractional Calculus with Applications. Volume 2. Fractional Differential Equations (De Gruyter, Berlin, 2019), pp. 385-406.
  * (6) R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Reports, 339, 1-78 (2000).
  * (7) I.M. Sokolov, J. Klafter and A. Blumen, Fractional kinetics, Physics Today 55, 48-54 (2002).
  * (8) R. Metzler and J. Klafter, The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics, J. Phys. A: Math. Gen. 37, R161–R208 (2004).
  * (9) R. Klages, G. Radons and I.M. Sokolov (eds.), Anomalous Transport: Foundations and Applications (Wiley VCH - Verlag, Weinheim, 2004).
  * (10) J. Klafter, S.C. Lim and R. Metzler (eds.), Fractional dynamics: recent advances (World Scientific, Singapore, 2012).
  * (11) A.V. Chechkin, R. Gorenflo, and I.M. Sokolov, Retarding subdiffusion and accelerating superdiffusion governed by distributed-order fractional diffusion equations, Phys. Rev. E 66, 046129 (2002).
  * (12) A.V. Chechkin, R. Gorenflo, I. M. Sokolov, and V. Y. Gonchar, Distributed order time fractional diffusion equation, Fractional Calculus and Applied Analysis 6, 259-280 (2003).
  * (13) A. V. Chechkin, J. Klafter and I. M. Sokolov, Fractional Fokker-Planck equation for ultraslow kinetics, Europhys. Lett. 63, 326-332 (2003).
  * (14) I. M. Sokolov, A. V. Chechkin, and J. Klafter, Distributed-order fractional kinetics, Acta Physica Polonica B 35, 1323-134 (2004).
  * (15) M. Naber, Distributed order fractional sub-diffusion, Fractals 12, 23-32 (2004).
  * (16) I. M. Sokolov and J. Klafter, From diffusion to anomalous diffusion: A century after Einstein’s Brownian motion, Chaos 15, 026103 (2005).
  * (17) I. M. Sokolov and A. V. Chechkin, Anomalous diffusion and generalized diffusion equations, Fluct. Noise Lett. 5, L275-L282 (2005).
  * (18) S. Umarov and R. Gorenflo, Cauchy and Nonlocal Multi-Point Problems for Distributed Order Pseudo-Differential Equations, Part One, J. Anal. Appl. 24, 449-466 (2005).
  * (19) M. M. Meerschaert and H.-P. Scheffler, Limit theorems for continuous time random walks with slowly varying waiting times, Statistics Prob. Lett. 71, 15-22 (2005).
  * (20) M. M. Meerschaert and H.-P. Scheffler, Stochastic model for ultraslow diffusion, Stoch. Proc. Appl. 116, 1215-1235 (2006).
  * (21) T. A. M. Langlands, Solution of a modified fractional diffusion equation, Physica A 367, 136-144 (2006).
  * (22) A. Hanyga, Anomalous diffusion without scale invariance, J. Phys. A: Math. Theor. 40, 5551 (2007).
  * (23) F. Mainardi and G. Pagnini, The role of the Fox–Wright functions in fractional sub-diffusion of distributed order, J. Comput. Appl. Math. 207, 245-257 (2007).
  * (24) F. Mainardi, G. Pagnini and R. Gorenflo, Some aspects of fractional diffusion equations of single and distributed order, Appl. Math. Comput. 187, 295 (2007).
  * (25) A.V. Chechkin, V.Yu. Gonchar, R. Gorenflo, N. Korabel and I.M. Sokolov, Generalized fractional diffusion equations for accelerating subdiffusion and truncated Lévy flights, Phys. Rev. E 78, 021111 (2008).
  * (26) A. N. Kochubei, Distributed order calculus and equations of ultraslow diffusion, J. Math. Anal. Appl. 340, 252-281 (2008).
  * (27) M. M Meerschaert, E. Nane, P. Vellaisamy, Distributed-order fractional diffusions on bounded domains. Journal of Mathematical Analysis and Applications 379, 216-228 (2011).
  * (28) A. V. Chechkin, I. M. Sokolov and J. Klafter, Natural and Modified Forms of Distributed-Order Fractional Diffusion Equations. In FracDyn , Chapter 5, pp. 107-127.
  * (29) I. M. Sokolov, A. V. Chechkin, and J. Klafter, Fractional diffusion equation for a power-law-truncated Lévy process, Physica A 336, 245-251 (2004).
  * (30) A. Stanislavsky, K. Weron, and A. Weron, Diffusion and relaxation controlled by tempered $\alpha$-stable processes, Phys. Rev. E 78, 061106 (2008)
  * (31) M. M. Meerschaert, Y. Zhang, B. Baeumer, Tempered anomalous diffusion in heterogeneous systems. Geophys. Res. Lett. 35, L17403 (2008).
  * (32) B. Baeumer, M. M. Meerschaert, Tempered stable Lévy motion and transient super-diffusion. J. Comp. Appl. Math. 233, 2438-2448 (2010).
  * (33) T. Sandev, A. Chechkin, H. Kantz, R. Metzler, Diffusion and Fokker-Planck equations with Generalized Memory Kernels (Survey paper). Frac. Calc. Appl. Anal. 18 (4), 1006-1038 (2015).
  * (34) A. Stanislavsky, K. Weron, Atypical case of the dielectric relaxation responses and its fractional kinetic equation. Frac. Calc. Appl. Anal. 19, 212 (2016).
  * (35) T. Sandev, I.M. Sokolov, R. Metzler, A. Chechkin, Beyond monofractional kinetics, Chaos, Solitons and Fractals, 102, 210-217 (2017).
  * (36) T. Sandev, R. Metzler, A. Chechkin, From continuous time random walks to the generalized diffusion equation. Frac. Calc. Appl. Anal. (Review paper) 21 (1), 10-28 (2018).
  * (37) A. Stanislavsky and A. Weron, Transient anomalous diffusion with Prabhakar-type memory. J. Chem. Phys. 149, 044107 (2018).
  * (38) T. Sandev, Z. Tomovski, J. L. A. Dubbeldam, and A. Chechkin, Generalized diffusion-wave equation with memory kernel, J. Phys. A: Math. Theor. 52, 015201 (2019)
  * (39) A. Stanislavsky, A. Weron, Control of the transient subdiffusion exponent at short and long times, Phys. Rev. Research 1, 023006 (2019).
  * (40) S. Chandrasekhar, Stochastic Problems in Physics and Astronomy, Rev. Mod. Phys. 15, 1 (1943).
  * (41) H. Fogedby, Langevin equations for continuous time Lévy flights, Phys. Rev. E 50, 1657 (1994).
  * (42) A. Baule and R. Friedrich, Joint probability distributions for a class of non-Markovian processes, Phys. Rev. E 71, 026101 (2005)
  * (43) A. Baule and R. Friedrich, A fractional diffusion equation for two-point probability distributions of a continuous-time random walk. Eur. Phys. Lett. 77, 10002 (2007).
  * (44) D. Kleinhans and R. Friedrich, Continuous-time random walks: Simulation of continuous trajectories. Phys. Rev. E 76, 061102 (2007).
  * (45) A. V. Chechkin, M. Hofmann, and I. M. Sokolov, Continuous-time random walk with correlated waiting times, Phys. Rev. E 80, 031112 (2009).
  * (46) T. Sandev, A. V. Chechkin, N. Korabel, H. Kantz, I. M. Sokolov, and R. Metzler, Distributed-order diffusion equations and multifractality: models and solutions, Phys. Rev. E 92, 042117 (2015)
  * (47) W. Feller, An Introduction to Probability Theory and Its Applications, vol. 1 and 2, John Wiley & Sons, 1968
  * (48) M. M. Meerschaert, D. A. Benson, H.-P. Scheffler, and B. Baeumer, Stochastic solution of space-time fractional diffusion equations, Phys. Rev. E 65, 041103 (2002).
  * (49) A. A. Stanislavsky, Subordinated Brownian Motion and its Fractional Fokker–Planck Equation, Phys. Scr. 67, 265 (2003).
  * (50) R. Gorenflo, F. Mainardi, and A. Vivoli, Continuous-time random walk and parametric subordination in fractional diffusion, Chaos, Solitons and Fractals 34, 87 - 103 (2007).
  * (51) R. Gorenflo and F. Mainardi, Subordination pathways to fractional diffusion, Eur. Phys. J.: Special Topics 193, 119 - 132 (2011).
  * (52) R. Gorenflo and F. Mainardi, Parametric subordination in fractional diffusion processes, in: S.C. Lim, J. Klafter and R. Metzler, eds., Fractional Dynamics, Recent Advances, World Scientific, Singapore, 2012 (Ch.10, pp. 227-261)
  * (53) A.I. Saichev and G.M. Zaslavsky, Fractional kinetic equations: Solutions and applications, Chaos, 7, 753 - 764 (1997).
  * (54) E. Barkai, Fractional Fokker-Planck equation, solution, and application. Phys. Rev. E 63, 046118 (2001).
  * (55) B. Baeumer and M. M. Meerschaert, Stochastic solutions for fractional Cauchy Problems, Frac. Calc. Appl. Anal, 4, 481-500 (2001).
  * (56) B. Baeumer, D. A. Benson, M. M. Meerschaert, and S. W. Wheatcraft, Subordinated advection-dispersion equation for contaminant transport, Water. Res. Research 37, 1543 - 1550 (2001).
  * (57) I.M. Sokolov, Solutions of a class of non-Markovian Fokker-Planck equations, Phys. Rev. E 66, 041101 (2002).
  * (58) F. Mainardi, G. Pagnini, and R. Gorenflo, Mellin transform and subordination laws in fractional diffusion processes, Fract. Calc. Appl. Anal., 6, 441 (2003).
  * (59) A.I. Saichev and S.G. Utkin, Random Walks with Intermediate Anomalous-Diffusion Asymptotics, J. Exp. Theor. Phys., 99, 443 - 448 (2004)
  * (60) A. Piryatinska, A.I. Saichev, and W.A. Woyczynski, Models of anomalous diffusion: the subdiffusive case, Physica A 349, 375 - 420 (2005).
  * (61) F. Mainardi, G. Pagnini, and R. Gorenflo, Mellin convolution for subordinated stable processes, J. Math. Sci. 132 637 (2006).
  * (62) R. Gorenflo, F. Mainardi, A. Vivoli, Continuous-time random walk and parametric subordination in fractional diffusion, Chaos, Solitons & Fractals, 34 87 (2007).
  * (63) M. Magdziarz, Langevin Picture of Subdiffusion with Infinitely Divisible Waiting Times, J. Stat. Phys. 135, 763–772 (2009).
  * (64) M. Magdziarz, Stochastic representation of subdiffusion processes with time-dependent drift, Stoch. Proc. Appl. 119, 3238–3252 (2009).
  * (65) J. Gajda and M. Magdziarz, Fractional Fokker-Planck equation with tempered $\alpha$-stable waiting times: Langevin picture and computer simulation, Phys. Rev. E 82, 011117 (2010).
  * (66) M. M. Meerschaert, P. Straka, Inverse Stable Subordinators, Math. Model. Nat. Phenom. 8, 1-16 (2013).
  * (67) A. Stanislavsky, K. Weron, and A. Weron, Anomalous diffusion with transient subordinators: A link to compound relaxation laws, J. Chem. Phys. 140, 054113 (2014).
  * (68) M. Annunziato, A. Borzi, M. Magdziarz, and A. Weron, A fractional Fokker–Planck control framework for subdiffusion processes, Optim. Control. Appl. Meth. 37, 290-304 (2015).
  * (69) M.Magdziarz and R. L. Schilling, Asymptotic properties of Brownian motion delayed by inverse subordinators, Proc. Amer. Math. Soc. 143, 4485–4501 (2015).
  * (70) M. Magdziarz and T. Zorawik, Stochastic representation of fractional subdiffusion equation. The case of infinitely divisible waiting times, Lévy noise and space-time-dependent coefficients, Proc. Amer. Math. Soc. 144, 1767-1778 (2016).
  * (71) A. Stanislavsky and A. Weron, Control of the transient subdiffusion exponent at short and long times, Phys. Rev. Research 1, 023006 (2019).
  * (72) A. Stanislavsky and A. Weron, Accelerating and retarding anomalous diffusion: A Bernstein function approach, Phys. Rev. E 101, 052119 (2020).
  * (73) 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).
  * (74) V. Sposini, A.V. Chechkin, F. Seno, G. Pagnini, R. Metzler, Random diffusivity from stochastic equations: comparison of two models for Brownian yet non-Gaussian diffusion, New J. Phys. 20, 043044 (2018).
  * (75) V. Sposini, A. Chechkin and R. Metzler, First passage statistics for diffusing diffusivity, J. Phys. A: Math. Theor. 52, 04LT01 (2019).
  * (76) B. Wang, S.M. Antony, S.C. Bae and S. Granick, Anomalous yet Brownian, PNAS 106, 15160 (2009).
  * (77) B. Wang, J. Kuo, S.C. Bae and S. Granick, When Brownian diffusion is not Gaussian, Nature Materials 11, 481 (2012).
  * (78) M.V. Chubynsky and G.W. Slater, Diffusing Diffusivity: A Model for Anomalous, yet Brownian, Diffusion, Phys. Rev. Lett. 113, 098302 (2014).
  * (79) R. Jain and K.L. Sebastian, Diffusion in a Crowded, Rearranging Environment, J. Phys. Chem. B 120, 3988 (2016).
  * (80) N. Tyagi and B. J. Cherayil, Non-Gaussian Brownian diffusion in dynamically disordered thermal environments, J. Phys. Chem. B 121, 7204 (2017).
  * (81) Y. Lanoiselee and D. S. Grebenkov, A model of non-Gaussian diffusion in heterogeneous media, J. Phys. A: Math. Theor. 51, 145602 (2018).
  * (82) Y. Lanoiselee, N. Moutal, and D.S. Grebenkov, Diffusion-limited reactions in dynamic heterogeneous media, Nature Comm. 9, 4398 (2018).
  * (83) S. Song, S. J. Park, M. Kim, J. S. Kim, B. J. Sung, S. Lee, J.-H. Kim, and J. Sung, Transport dynamics of complex fluids, PNAS 116 (26), 12733-12742 (2019).
  * (84) J. M. Miotto, S. Pigolotti, A. V. Chechkin, and S. Roldan-Vargas, Length scales in Brownian yet non-Gaussian dynamics, arXiv:1911.07761.
  * (85) R. Metzler, Superstatistics and non-Gaussian diffusion, Eur. Phys. J.: Spec. Top. 229, 711-728 (2020).
  * (86) I. M. Sokolov and J. Klafter, Field-induced dispersion in subdiffusion, Phys. Rev. Lett. 97, 140602 (2006).
  * (87) A. V. Chechkin, R. Gorenflo, and I. M. Sokolov, Fractional diffusion in inhomogeneous media, J. Phys. A: Math. and Gen. 38 L679 (2005).
  * (88) E. Barkai and R. J. Silbey, Fractional Kramers Equation, J. Phys. Chem. B 104, 3866-3874 (2000).
  * (89) I. M. Sokolov, Thermodynamics and fractional Fokker-Planck equations, Phys. Rev. E 63, 056111, (2001).
  * (90) I. M. Sokolov, Lévy flights from a continuous-time process, Phys. Rev. E 63, 011104 (2000).
  * (91) Equation (30) is always fulfilled for independent random variables. However, there exist the examples when it is also valid for dependent variables, see e.g. G. G. Hamedani, Sub-Independence - An explanatory Perspective, Commun. in Stat. Theory and Methods, 42, 3615-3638 (2013).
  * (92) M. Abramovitz and I.A. Stegun, Handbook of mathematical functions, Dover, N.Y., 1972