# Fractional Dissipative PDEs

Franz Achleitner Goro Akagi Christian Kuehn
Jens Markus Melenk Jens D.M. Rademacher Cinzia Soresina Jichen Yang TU Wien,
Institute for Analysis and Scientific Computing, 1040 Vienna, Austria,
<EMAIL_ADDRESS>University, Mathematical Institute and
Graduate School of Science, 980-8578 Sendai, Japan,
<EMAIL_ADDRESS>University of Munich, Department of
Mathematics, Research Unit “Multiscale and Stochastic Dynamics”, 85748
Garching b. München, Germany<EMAIL_ADDRESS>University of Munich,
Munich Data Science Institute (MDSI), 85748 Garching b. München,
GermanyExternal Faculty, Complexity Science Hub Vienna, 1080 Vienna, AustriaTU
Wien, Institute for Analysis and Scientific Computing, 1040 Vienna, Austria,
<EMAIL_ADDRESS>Hamburg, MIN faculty, Department of
Mathematics, Applied Dynamical Systems, 20146 Hamburg, Germany,
<EMAIL_ADDRESS>Graz, Department of Mathematics and
Scientific Computing, 8010 Graz, Austria<EMAIL_ADDRESS>Engineering University, College of Mathematical Sciences, 150001 Harbin,
China<EMAIL_ADDRESS>
###### Abstract

In this chapter we provide an introduction to fractional dissipative partial
differential equations (PDEs) with a focus on trying to understand their
dynamics. The class of PDEs we focus on are reaction-diffusion equations but
we also provide an outlook on closely related classes of PDEs. To simplify the
exposition, we only discuss the cases of fractional time derivatives and
fractional space derivatives in the PDE separately. As our main tools, we
describe analytical as well as numerical methods, which are generically
necessary to study nonlinear dynamics. We start with the analytical study of
steady states and local linear stability for fractional time derivatives. Then
we extend this view to a global perspective and consider time-fractional PDEs
and gradient flows. Next, we continue to steady states, linear stability
analysis and bifurcations for space-fractional PDEs. As a final analytical
consideration we discuss existence and stability of traveling waves for space-
fractional PDEs. In the last parts, we provide numerical discretization
schemes for fractional (dissipative) PDEs and we utilize these techniques
within numerical continuation in applied examples of fractional reaction-
diffusion PDEs. We conclude with a brief summary and outlook on open questions
in the field.

##  1 Introduction

Reaction-diffusion partial differential equations (PDEs) are a crucial
modeling tool appearing in all the natural sciences, engineering and beyond
[1, 2, 3]. In particular, reaction-diffusion PDEs are relevant at very small
length scales up to very large macroscopic length scales. The existence,
uniqueness and regularity theory for large classes of reaction-diffusion PDEs
is quite classical by now and many aspects are well-understood [4, 5, 6, 7,
8]. Yet, asking for a more precise description of the qualitative and
quantitative dynamics of quite small classes of reaction-diffusion PDEs, a lot
of questions appear that have led to very active research areas. There are few
absolutely complete and comprehensive results about the dynamics of reaction-
diffusion PDEs. This is expected as such a dynamical classification does not
even exist for nonlinear ODEs formed by the reaction terms. Yet, for the case
with enough dissipativity, there is a well-established strategy that can be
applied. First, one uses dissipativity to establish boundedness of solutions,
often obtaining the existence of a (finite-dimensional) attractor as a
consequence of global dissipation [9, 10]. Unfortunately, just having the
existence of an attractor yields relatively little information on the precise
dynamics, so one then proceeds to study finer information via invariant sets
and patterns [11, 12, 13]. A first step in this direction is the analysis of
steady states. Often one starts with homogeneous steady states and then
proceeds to more complex steady states using complementary methods, e.g.,
bifurcation theory [14, 15, 16] and variational methods [17]. Once steady
states are found, one continues to investigate their local linear stability
via classical methods such as spectral theory [18]. A next step is the
investigation of non-stationary patterns such as traveling waves [12], which
are effectively invariant steady sets in a moving frame. Again, one then
studies the stability of traveling waves [19, 20]. To cross-validate and
globally extend dynamical results, it is often imperative to combine
analytical approaches with numerical techniques, which are based upon spatial
discretization combined with temporal simulation or numerical continuation
methods [21, 22, 23].

To illustrate this classical paradigm, let us consider a typical example given
by the PDE

$\partial_{t}u=\partial_{x}^{2}u+u(p_{1}-u)(u-p_{2})=:\partial_{x}^{2}u+f(u;p),\quad
u_{0}(x)=u(0,x),$ (1)

where $u=u(t,x)$, $x\in\Omega\subseteq\mathbb{R}$, $t\geq 0$, and
$p=(p_{1},p_{2})\in\mathbb{R}^{2}$ are parameters. The reaction-diffusion PDE
(1) has many different names in different fields: Fisher equation with Allee
effect [3], Allen-Cahn equation [24], Chafee-Infante equation [25], Nagumo
equation [26], real Ginzburg-Landau equation [27], Schlögl model [28], and
(overdamped) $\Phi^{4}$ model [29]. This clearly demonstrates the importance
of (1) as a key example, or even as a normal form [30], for many dynamical
phenomena. One easily checks that (1) is dissipative due to the $-u^{3}$ term
in the nonlinearity [9, 10]. For the case of a bounded domain $\Omega=(a,b)$,
$a<b$, one obtains a finite-dimensional attractor [31, 32, 10, 9]. For
concreteness, let us consider homogeneous Neumann boundary data

$\frac{\partial u}{\partial x}(a)=0=\frac{\partial u}{\partial x}(b).$ (2)

In this case, we immediately spot the three homogeneous equilibria
$u_{*,0}(x)\equiv 0$, $u_{*,1}=p_{1}$, $u_{*,2}=p_{2}$. The standard
configuration for the Allen-Cahn equation is to fix $p_{2}=-p_{1}$ and set
$p_{1}^{2}=:p_{AC}$ so that the nonlinearity becomes $p_{AC}u-u^{3}$. Even in
this simple case, the homogeneous state $u_{*,0}$ undergoes bifurcations as
$p_{AC}$ is increased generating non-homogenous/patterned steady states in the
process [11, 15, 12]. Note that linear stability of any steady state on a
$\Omega=(a,b)$ reduces to the study of the eigenvalue problem

$({\textnormal{D}}_{u}f)(u_{*,\cdot};p)U=\lambda U,\qquad U=U(x),$

where ${\textnormal{D}}_{u}$ denotes the Fréchet derivative of $f$ on a
suitable function space $\mathcal{X}$, e.g., one could consider
$\mathcal{X}=H^{1}(\Omega)=W^{1,2}(\Omega)$. Computation and visualization of
the non-homogeneous steady states is possible locally near bifurcation points.
Globally, the situation is already very difficult and it helps to employ
numerical tools [21, 22, 23]. A natural step is to spatially discretize the
steady state problem, e.g., using finite elements, finite differences or a
spectral Galerkin method. Then one can apply path-following techniques to
generate global bifurcation diagrams [21, 22, 23].

For the case of an unbounded domain $\Omega=\mathbb{R}$ and the standard
Nagumo parameter choices $p_{1}=1$ and $p_{2}=p_{Nag}\in(0,1)$, it is well-
known that the PDE (1) has traveling wave solutions [33, 19, 12] of the form

$u(t,x)=u(x-st)=u(\zeta)$ (3)

for a wave speed $s=s(p_{Nag})$. In fact, the ansatz (3) just results in a
second-order ODE. One can prove that this ODE has a parametric family of
heteroclinic orbits corresponding to traveling front solutions and a
parametric family of homoclinic orbits corresponding to traveling pulse
solutions. One can then prove using linearization and Sturm-Liouville theory
that the traveling front is linearly (and in fact also nonlinearly) stable,
while the traveling pulse is already linearly unstable [33, 19, 12]. To
visualize the traveling wave dynamics, one option is to just combine spatial
discretization on a very large domain $(a,b)$ with $a\ll-1$ and $b\gg 1$ with
a time discretization. Alternatively, one could freeze the wave by using a
suitable reference frame [34].

Although the theory of (1) is very detailed, an important question for
applications is its robustness. This not only means proving stability results
in phase space but also refers to structural stability of the model itself,
i.e., if the PDE class is changed and/or extended, do dynamical phenomena
persist or do new effects appear? Before answering this question, we have to
consider the input from applications. Which extensions to reaction-diffusion
PDEs do occur in practice? One answer is provided by going back to the
microscopic derivation of diffusion. In particular, the Laplacian $\Delta$ as
well as the associated heat operator $\partial_{t}-\Delta$ can arise from
taking a macroscopic limit of a random walk described by a stochastic process
$z(t)$. If the stochastic process satisfies a power law for the mean-squared
displacement

$\mathbb{E}[z(t)^{2}]\sim 2\kappa_{\alpha}t^{\alpha},$

then one easily calculates that a standard random walk has the exponent
$\alpha=1$. The case $\alpha=1$ could also be called “classical” diffusion.
Classical diffusion is a Markovian and a Gaussian/Wiener process. Note that
$\kappa_{1}=:\sigma$ becomes the diffusion constant of the mean-field limit
diffusion equation

$\partial_{t}u=\sigma\Delta u$

for the probability density $u(t,x)$ of the random walker to be at position
$x$ at time $t$. Yet, if $\alpha>1$, then the behavior is super-diffusive (or
“ballistic”), while $\alpha\in(0,1)$ gives subdiffusive dynamics. The cases
$\alpha\neq 1$ are also called anomalous diffusion [35] and lead to different
macroscopic operators beyond the standard heat equation; we will define these
fractional differential operators, as well as several additional ones, in
Section 2. Once all these additional fractional operators arising from various
applications have been collected, we can try to proceed to study dissipative
fractional PDEs along the same lines as the classical example (1) presented
above.

The remaining part of this review is structured as follows. In Section 2 we
fix the notation we use for this chapter. In Section 3 we study the extension
of steady states and linear stability with a focus on fractional time
derivatives. In Section 4 we take a global perspective and consider time-
fractional reaction-diffusion equations and their associated gradient flows.
In Section 5 we briefly recap linear stability analysis and bifurcations for
space-fractional derivatives. In Section 6, we survey results for traveling
waves for space-fractional PDEs. In Section 7, we give an overview for
numerical discretization schemes for fractional (dissipative) PDEs, while in
Section 8, we show how these techniques can be embedded into numerical
continuation in practical examples of fractional reaction-diffusion PDEs
arising in applications. In Section 9, we provide a brief summary and outlook
towards interesting open problems.

We stress that our review provides only a partial snapshot of the field of
dynamics of fractional dissipative PDEs. It is probably already impossible to
cover the field comprehensively as it has grown so quickly in recent years.
Yet we hope that this chapter makes it much easier to find an entry point to
the results and techniques in the field, even beyond the results we mention
here. Furthermore, we believe that the snapshot of results discussed here
clearly indicates that the general approach to study dissipative reaction-
diffusion PDEs can be lifted from classical to fractional derivatives. Yet,
the techniques and phenomena are often substantially different, which makes it
an exciting area to explore.

##  2 Definitions and Notation

Let $\Omega$ now denote a domain in $\mathbb{R}^{d}$. Consider
$x\in\Omega\subseteq\mathbb{R}^{d}$, $t\geq 0$ and write
$u=u(t,x)\in\mathbb{R}^{N}$ for the (unknown) solution of the PDEs we study
here; often we will just take $d=1$ and $N=1$ as the main basic case. $U$ will
be used to highlight a linearized problem $U=U(t,x)\in\mathbb{R}^{N}$.
$\Delta$ denotes the classical Laplacian, $\partial_{x}^{2}$ the classical
Laplacian in one dimension, and $\partial_{t}$ the classical partial time
derivative.

###  2.1 Time-Fractional Derivative(s)

First, we define the Riemann-Liouville fractional integral and derivative for
the subdiffusion range $\alpha\in(0,1)$. Let $f(t)\in\mathnormal{L}^{1}(0,T)$
for any $T>0$. The fractional Riemann-Liouville integral of order $\alpha$ is
defined as

$(\mathfrak{D}_{0,t}^{-\alpha}f)(t):=(k_{\alpha}\ast
f)(t)=\frac{1}{\Gamma(\alpha)}\int_{0}^{t}(t-s)^{\alpha-1}f(s)\,\mathrm{d}s,\quad\alpha\in(0,1),$
(4)

where the (Laplace) convolution kernel is
$k_{\alpha}(t):=t^{\alpha-1}/\Gamma(\alpha)$, and $\Gamma(\alpha)$ is the
Gamma function. We can also generalize the time dependence of the integrals
yielding the natural definitions

$\displaystyle\mathfrak{D}^{-\alpha}_{t_{0},t}f(t):=\frac{1}{\Gamma(\alpha)}\int_{t_{0}}^{t}(t-s)^{\alpha-1}f(s)\,\mathrm{d}s,\qquad\mathfrak{D}^{-\alpha}_{t,T}f(t):=\frac{1}{\Gamma(\alpha)}\int_{t}^{T}(s-t)^{\alpha-1}f(s)\,\mathrm{d}s.$

Notably, for $\alpha=1$ we get
$(\mathfrak{D}_{0,t}^{-\alpha}f)(t)=\int_{0}^{t}f(s)\mathrm{d}s$. For
$f:[0,T]\to\mathbb{R}$ the Riemann-Liouville fractional derivative of order
$\alpha$ is defined as

$(\mathfrak{D}_{0,t}^{\alpha}f)(t):=\frac{\mathrm{d}}{\mathrm{d}t}(k_{1-\alpha}\ast
f)(t)=\frac{1}{\Gamma(1-\alpha)}\frac{\mathrm{d}}{\mathrm{d}t}\int_{0}^{t}(t-s)^{-\alpha}f(s)\,\mathrm{d}s,\quad\alpha\in(0,1),$
(5)

with the obvious generalizations to other endpoints $t_{0}$, $T$, which is
then reflected in the notation by adding suitable indices. Notably, the
Riemann-Liouville fractional derivative is nonzero on constants,

$(\mathfrak{D}_{0,t}^{\alpha}1)(t)=k_{1-\alpha}(t),$

which tends to zero as $t^{-\alpha}$ for $t\to\infty$ and is unbounded for
$t\to 0$. At $\alpha=0$, $\mathfrak{D}_{0,t}^{\alpha}$ is the identity
operator, i.e., $(\mathfrak{D}_{0,t}^{\alpha}f)(t)=f(t)$, while for
$\alpha=1$, $(\mathfrak{D}_{0,t}^{\alpha}f)(t)$ formally yields
$(\mathrm{d}/\mathrm{d}t)\left(\delta\ast f\right)(t)=f^{\prime}(t)$ with
Dirac delta distribution $\delta(t)$ and
${}^{\prime}:=\mathrm{d}/\mathrm{d}t$. A simple sufficient condition for the
existence of the Riemann-Liouville fractional derivative is as follows.

######  Lemma 2.1 ([36, Lemma 2.2]).

Let $f(t)\in\mathnormal{AC}([0,T])$ and $\alpha\in(0,1)$. Then
$(\mathfrak{D}_{0,t}^{\alpha}f)(t)$ exists almost everywhere. Moreover
$(\mathfrak{D}_{0,t}^{\alpha}f)(t)\in\mathnormal{L}^{p}(0,T)$, $1\leq
p<1/\alpha$, and

$\displaystyle(\mathfrak{D}_{0,t}^{\alpha}f)(t)=(k_{1-\alpha}\ast
f^{\prime})(t)+k_{1-\alpha}(t)f(0).$

Here $\mathnormal{AC}([0,T])$ is the space of functions $f$ which are
absolutely continuous on $[0,T]$, i.e., $f(t)=c+\int_{0}^{t}g(s)\mathrm{d}s$
for some $g\in\mathnormal{L}^{1}(0,T)$ and constant $c$. This lemma motivates
the definition of the (Djrbashian-)Caputo fractional derivative. For
$f\in\mathnormal{AC}([0,T])$ the _(Djrbashian-)Caputo fractional derivative_
of order $\alpha$ is defined as

$(\tilde{\mathfrak{D}}_{0,t}^{\alpha}f)(t):=(k_{1-\alpha}\ast
f^{\prime})(t)=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-s)^{-\alpha}f^{\prime}(s)\,\mathrm{d}s,\quad\alpha\in(0,1).$
(6)

Hence, for sufficiently smooth $f:[0,T]\to\mathbb{R}$, Lemma 2.1 gives a well-
known relation between Riemann-Liouville and Caputo derivatives

$\displaystyle\tilde{\mathfrak{D}}_{0,t}^{\alpha}f(t)=\mathfrak{D}_{0,t}^{\alpha}(f(t)-f(0)).$

One can extend these definitions also beyond the case $\alpha\in(0,1)$ to any
$\alpha>0$. Define $m:=\lceil\alpha\rceil\in{\mathbb{N}}_{0}$ for
$\alpha\not\in{\mathbb{N}}_{0}$ and $m:=\alpha+1$ for
$\alpha\in{\mathbb{N}}_{0}$. Then, the Riemann-Liouville derivatives are given
by

$\displaystyle\mathfrak{D}^{\alpha}_{0,t}f(t)$
$\displaystyle:=\frac{\mathrm{d}^{m}}{\mathrm{d}t^{m}}\mathfrak{D}^{-(m-\alpha)}_{0,t}f(t),$
$\displaystyle\mathfrak{D}^{\alpha}_{t,T}f(t)$
$\displaystyle:=(-1)^{m}\frac{\mathrm{d}^{m}}{\mathrm{d}t^{m}}\mathfrak{D}^{-(m-\alpha)}_{t,T}f(t).$

The left and right Caputo derivatives are given by

$\displaystyle\mathfrak{\tilde{D}}^{\alpha}_{0,t}f(t)$
$\displaystyle:=\mathfrak{D}^{-(m-\alpha)}_{0,t}f^{(m)}(t),$
$\displaystyle\mathfrak{\tilde{D}}^{\alpha}_{t,T}f(t)$
$\displaystyle:=(-1)^{m}\mathfrak{D}^{-(m-\alpha)}_{t,T}f^{(m)}(t),$

where one requires sufficient regularity, in terms of a generalization of
Lemma 2.1 to ensure existence of the derivatives. The Riemann-Liouville and
the Caputo fractional derivatives are related to each other: For sufficiently
smooth $f$, one has

$\displaystyle\mathfrak{D}^{\alpha}_{0,t}f(t)$
$\displaystyle=\mathfrak{\tilde{D}}^{\alpha}_{0,t}u+\sum_{k=0}^{m-1}\frac{f^{(k)}(0)(t-0)^{k-\alpha}}{\Gamma(k+1-\alpha)}.$

###  2.2 Space-Fractional Derivative(s)

We start with the case of the spatial domain $\Omega=\mathbb{R}$. In
experiments asymmetric superdiffusive behavior has been observed in many
applications [37, 38, 39, 40, 41] (see also the references in [42, 43] and
[44, Part 3]). Using continuous time random walks (CTRW) with an asymmetric
distribution of jump lengths allows one to derive evolution equations with
asymmetric Lévy operators [45, §3.3.3]. In this context, a very general class
of space-fractional derivatives is derived microscopically from Lévy
operators. A subclass of Lévy operators are the Riesz–Feller operators, which
are the infinitesimal generators of stable Lévy processes; see [46, §14] and
[47]. In the one-dimensional setting, Riesz–Feller operators can be
characterized by two parameters: the index of stability $\gamma$ and another
parameter $\theta$ measuring the asymmetry. Riesz–Feller operators can be
defined as Fourier multiplier operators

$\mathcal{F}[\operatorname{D}_{\theta}^{\gamma}u](\xi)=\psi^{\gamma}_{\theta}(\xi)\
\mathcal{F}[u](\xi),\qquad\xi\in\mathbb{R},$

with $\mathcal{F}$ denoting the Fourier transform and symbol

$\psi^{\gamma}_{\theta}(\xi)=-|\xi|^{\gamma}\exp\left[{\textnormal{i}}(\mathrm{sgn}(\xi))\theta\tfrac{\pi}{2}\right]\quad\text{
for }\,0<\gamma\leq 2\,\text{ and }\,|\theta|\leq\min\\{\gamma,2-\gamma\\}.$

We will require enough regularity for $u=u(x)$ with $x\in\mathbb{R}$ so that
we can define its Fourier transform. Special cases for Riesz-Feller operators
are the second order derivative $\operatorname{D}_{0}^{2}=\partial^{2}_{x}{}$,
the (negative) fractional Laplacians
$\operatorname{D}_{0}^{\gamma}=-(-\partial^{2}_{x})^{\gamma/2}$ for
$0<\gamma\leq 2$ and $\theta=0$ (which are the only symmetric Riesz–Feller
operators), and Weyl fractional derivatives
$\operatorname{D}_{2-\gamma}^{\gamma}$ for $0<\gamma<2$ and $\theta=2-\gamma$
(which are extremely asymmetric Riesz–Feller operators).

For the higher-dimensional setting $\Omega=\mathbb{R}^{d}$, we will only need
the (positive) fractional Laplacian $(-\Delta_{x})^{\gamma/2}$,
$\gamma\in(0,2)$. In full space fractional Laplacians
$(-\Delta_{x})^{\gamma/2}$ can be defined in various equivalent ways, see,
e.g., [48]: For example, again as Fourier multiplier operators

$\mathcal{F}[(-\Delta_{x})^{\gamma/2}u](\xi)=|\xi|^{\gamma}\
\mathcal{F}[u](\xi),\qquad\xi\in\mathbb{R}^{d},$

with symbol $|\xi|^{\gamma}$. Or, we can define them as singular integral
operators in two different ways:

$(-\Delta_{x})^{\gamma/2}u(x):=c_{d,\gamma}\int_{\mathbb{R}^{d}}\,{\frac{-u(x+y)-u(x-y)+2u(x)}{2|y|^{d+\gamma}}}\,\,\operatorname{d}\\!y\,=c_{d,\gamma}P.V.\int_{\mathbb{R}^{d}}\frac{u(x)-u(y)}{|x-y|^{d+\gamma}}\,dy,$
(7)

where P.V. indicates that a principal value integral is considered, and the
normalizing constant

$c_{d,\gamma}:=\frac{2^{\gamma}\Gamma\big{(}\tfrac{d+\gamma}{2}\big{)}}{\pi^{d/2}\big{|}\Gamma\big{(}-\tfrac{\gamma}{2}\big{)}\big{|}}$
(8)

is such that the singular integral operators have the Fourier symbol
$|\xi|^{\gamma}$. On bounded domains, the situation becomes more complicated
and one has to carefully distinguish different definitions for the fractional
Laplacian, e.g., the spectral fractional Laplacian (SFL) and the integral
fractional Laplacian (IFL). These operators coincide on ${\mathbb{R}}^{d}$,
but they differ when restricted to bounded domains $\Omega$. To define the
Dirichlet spectral fractional Laplacian let $\Omega\subset{\mathbb{R}}^{d}$ be
a bounded Lipschitz domain and let
$(\phi_{i},\lambda_{i})_{i\in{\mathbb{N}}}\subset
H^{1}_{0}(\Omega)\setminus\\{0\\}\times{\mathbb{R}}$ be the eigenpairs of

$\Delta\phi=\lambda\phi\quad\mbox{ in $\Omega$},\quad\phi|_{\partial\Omega}=0$
(9)

with the eigenvalues $\lambda_{i}\rightarrow-\infty$ sorted in descending
order and normalized such that $\|\phi_{i}\|_{L^{2}(\Omega)}=1$,
$i\in{\mathbb{N}}$. Note that $(\phi_{i})_{i\in{\mathbb{N}}}$ forms an
orthonormal basis of $L^{2}(\Omega)$. The (positive) Dirichlet spectral
fractional Laplacian $(-\Delta)^{\gamma/2}$ is defined as the operator
$u\mapsto(-\Delta)^{\gamma/2}u(x):=\sum_{i}(-\lambda_{i})^{\gamma/2}(u,\phi_{i})_{L^{2}(\Omega)}\phi_{i}(x)$.
Note that $\gamma=0$ reproduces the identity operator and $\gamma=2$ the
standard Laplacian. The operator is an isomorphism
$\widetilde{H}^{\gamma/2}(\Omega)\rightarrow
H^{-\gamma/2}(\Omega):=(\widetilde{H}^{\gamma/2}(\Omega))^{\prime}$, where for
$\beta\in[0,1]$ the fractional Sobolev space $\widetilde{H}^{\beta}(\Omega)$
is the closure of $C^{\infty}_{0}(\Omega)$ under the norm

$\|u\|_{\widetilde{H}^{\beta}(\Omega)}:=\left(\sum_{i}(-\lambda_{i})^{\beta}|(u,\phi_{i})_{L^{2}(\Omega)}|^{2}\right)^{1/2}.$
(10)

(See [49, Chap. 3] for alternative, equivalent norms on
$\widetilde{H}^{\beta}(\Omega)$.) The Dirichlet spectral fractional Laplacian
is based on Dirichlet boundary conditions in (9). Replacing the Dirichlet
conditions in (9) by Neumann boundary conditions leads to a similar operator,
the Neumann spectral fractional Laplacian. Throughout the text we will
abbreviate

$\Delta^{\gamma/2}:=-(-\Delta)^{\gamma/2}$

for the integral and the spectral fractional (both Dirichlet and Neumann)
Laplacian. Note that this is consistent with the 1D operator
$\operatorname{D}_{0}^{\gamma}$ introduced earlier.

######  Remark.

In the literature there is sometimes confusion regarding the signs of “the”
fractional Laplacian. This issue already arises for the classical Laplacian as
some authors prefer to have it as “negative” operator while others prefer a
“positive” operator. This is just a matter of sign conventions. In dynamical
systems it is more common to have the classical Laplacian with nonpositive
eigenvalues to indicate that there is decay for the heat semigroup in
accordance with the principle that stability is associated with negative real
parts for elements in the spectrum. So we will adopt this convention here,
which explains our sign conventions.

For further general references regarding fractional calculus, we refer to [50,
36], which cover the univariate case in significant detail. We refer to [44]
for a good survey of further theoretical and numerical aspects of fractional
differential operators.

##  3 Steady States & Linear Stability: Time Fractional PDEs

In case of subdiffusion $\alpha\in(0,1)$, the waiting time probability density
function from the continuous time random walk modelling is non-exponential and
yields a non-Markovian process with memory. Arguing via Fourier and Laplace
transforms, e.g., [35, 51, 52], one derives what we refer to here as the
subdiffusion equation for $u=u(t,x)\in\mathbb{R}$, $x\in\mathbb{R}$, $t>0$,

$\partial_{t}u=\mathfrak{D}_{0,t}^{1-\alpha}\sigma\partial_{x}^{2}u,$ (11)

with diffusion coefficient $\sigma$ and time-fractional Riemann-Liouville
derivative $\mathfrak{D}_{0,t}^{1-\alpha}$. This non-local convolution
operator in time with power law kernel entails that the system has memory.
Acting by the inverse of $\mathfrak{D}_{0,t}^{1-\alpha}$ on (11) yields the
distinct but mathematically equivalent form [53]

$\displaystyle\tilde{\mathfrak{D}}_{0,t}^{\alpha}u=\sigma\partial_{x}^{2}u,$
(12)

with the Caputo fractional derivative $\tilde{\mathfrak{D}}_{0,t}^{\alpha}$.
However, care has to be taken in the functional analytic interpretation of
(11) and (12), in particular regarding continuity of solutions at $t=0$ [52].
Indeed, replacing the left hand side of (12) with
$\mathfrak{D}_{0,t}^{\alpha}(u-u_{0})$ gives an alternative form, which
requires less smoothness of $u$ with $u(0,x)=u_{0}(x)$ [35, 54]. The Cauchy
problems of (11) and (12) are also equivalent to an integro-differential
equation that was first investigated via Mellin transforms [55], which we do
not discuss further here; more details on the comparisons can be found in
[56]. In all cases, the integral memory means that equations with subdiffusion
cannot be viewed as dynamical systems on the phase space of the natural
initial condition $u(0,x)$. Starting at time $t>0$ requires suitably
prescribing $u(s,x)$ on the elapsed time interval $0\leq s\leq t$, so that
such equations may be viewed as equations with monotonically increasing delay.
Moreover, as in general for non-autonomous equations, we expect that solutions
can cross the initial state and stability properties are not necessarily
determined by the spectrum of a time-independent linear operator.

###  3.1 The Linear Subdiffusion Equation

In order to highlight essential differences between classical diffusion and
subdiffusion we consider the basic subdiffusion equation in the Riemann-
Liouville form (11). Compared to the heat equation, it is more subtle to see
that solutions of (11) remain positive, and to determine the growth and decay
properties. We briefly outline some aspects and refer to, e.g., [35] for
further details and different perspectives involving the Fox-function. The
Green’s function of (11) has the Fourier transform [35, eq. 49]

$\displaystyle\mathcal{F}\Phi(t,q)=E_{\alpha}(-\sigma q^{2}t^{\alpha}),$ (13)

where $E_{\alpha}(z)=\sum_{n=0}^{\infty}z^{n}/\Gamma(1+n\alpha)$ is the
Mittag-Leffler function, and $q$ the Fourier-wavenumber. The function (13) has
the asymptotic behavior [57, eq. 20]

$\displaystyle E_{\alpha}(-\sigma
q^{2}t^{\alpha})\sim\begin{cases}\exp(-\frac{\sigma
q^{2}t^{\alpha}}{\Gamma(1+\alpha)}),&t\ll(\sigma q^{2})^{1/\alpha},\\\ (\sigma
q^{2}t^{\alpha}\Gamma(1-\alpha))^{-1},&t\gg(\sigma
q^{2})^{1/\alpha}.\end{cases}$ (14)

Here we see an effect of memory for nonzero wavenumber: exponential decay for
relatively short times, and algebraic decay for long times, here with power
$-\alpha$. The separation of these depends in particular on the wavenumber.
The algebraic decay effectively replaces the exponential decay of Fourier
modes for the classical heat equation, and illustrates that linear stability
for processes with subdiffusion has a different character from that for
classical diffusion. On the one hand, this behavior is also foreshadowed in
Fourier-Laplace space, for which we view the right-hand side of (11) realized
as an operator with space-time domain
$\mathnormal{AC}([0,T];\mathnormal{H}^{2}(\mathbb{R}))$ into
$\mathnormal{L}^{1}(0,T;\mathnormal{L}^{2}(\mathbb{R}))$. With the Laplace
transformation $\mathfrak{L}$ we obtain

$\displaystyle(\mathfrak{L}\mathcal{F}{\Phi})(s,q)=(s+\sigma
q^{2}s^{1-\alpha})^{-1},\quad{\textnormal{Re}}(s)>0,$ (15)

where $s\in\mathbb{C}$ is the Laplace space variable. Indeed, one expects that
the branch point of (15) at the origin creates algebraic decay in the inverse
Laplace transform. This analysis will be the topic of Section 3.3 for reaction
subdiffusion models and directly relates to decay in fractional linear ODE.
This is always algebraic with rate $-\alpha$, in contrast to classical ODE,
where decay is exponential for eigenvalues with negative real parts, cf. [58,
59]. On the other hand, in $x$-space a scaling analysis shows that the
parabolic scaling of the classical heat equation is replaced in (11) by a
fractional self-similar scaling: if $u(t,x)$ solves (11) then so does
$u(\varepsilon^{2/\alpha}t,\varepsilon x)$ for $\varepsilon\in\mathbb{R}$.
Based on the similarity variable $x/t^{\alpha/2}$ a series expansion of the
Green’s function can be found (e.g. [35, eq. 46], [60, eq. 4.23]),

$\Phi(t,x)=\frac{1}{\sqrt{4\sigma
t^{\alpha}}}\sum_{n=0}^{\infty}\frac{(-1)^{n}}{n!\Gamma\left(1-\alpha(n+1)/2\right)}\left(\frac{|x|}{\sqrt{\sigma
t^{\alpha}}}\right)^{n},\quad t>0.$ (16)

Via this form one can find that $\Phi$ has an expression in terms of the
Wright function [61] and is therefore positive [60] and decays algebraically
locally uniformly in $x$ for $t\gg 1$ with power $-\alpha/2$ [62, 63, 64]. In
addition, we remark that, using the method of separation of variables, it was
shown in [65] that solutions to the Cauchy problem (11) with homogeneous
Dirichlet boundary condition decay pointwise algebraically with power
$-\alpha$ for $t\gg 1$, i.e., faster than on the unbounded domain, which was
also found to be optimal in the more general case (18) below. A more detailed
discussion of the Mittag-Leffler and Wright functions, their estimates and
relation to fractional calculus can be found in [52]; we also refer to [66,
60, 67] for further details of the Wright functions serving as the Green’s
functions to various (space-)time fractional diffusion equations.

We also briefly outline some aspects of the subdiffusion equation in the
Caputo form (12). The formulation via $\tilde{\mathfrak{D}}_{0,t}^{\alpha}$
motivates to consider a more general form of time fractional evolution
equations with a convolution $(k\ast v)(t)=\int_{0}^{t}k(t-s)v(s)ds$, $t\geq
0$ for a kernel $k$, generalizing (12) to

$k\ast\partial_{t}u=\sigma\partial_{x}^{2}u.$ (17)

The theory for Cauchy-problems in (12) and (17) is by now quite well developed
and an excellent introduction is provided by [52]. For (17) we particularly
refer to [68] and [69], which provide decay estimates on bounded and unbounded
domains, respectively. For instance, in the subdiffusion case with homogeneous
Dirichlet boundary condition the decay estimate of solutions has a different
character, and exponential decay of Fourier modes is replaced by algebraic
decay with rate $-\alpha$, consistent with the separable solutions to (11);
moreover, [68] proves $\mathnormal{L}^{2}$-decay with a fixed rate $-\alpha$
for any space dimension. On unbounded domains, the character of algebraic
decay persists, but the $\mathnormal{L}^{2}$-decay depends on the spatial
dimension $d$: the rate is slower than that on the bounded domain for $d<4$
and is fixed at $-\alpha$ for $d\geq 4$, cf. [69]. Hence, the decay rate on
the bounded domain provides a lower bound for the case of the unbounded domain
in any space dimension. In contrast to the classical diffusion, the
subdiffusion is so slow that increasing the spatial dimensions (particularly
for $n>4$) and restricting to bounded domains does not increase the decay
rate. We also refer to [54] for various analytical aspects of time-fractional
diffusion equations. In addition, to the subdiffusion equation, various
examples of kernels $k$ can be found in [68]. Lastly, we remark that the
subdiffusion equations of the forms (12) and (17) can be extended as
“parabolic-type” equations by adding uniformly elliptic operators $L$ to the
fractional derivatives: $\tilde{\mathfrak{D}}_{0,t}^{\alpha}+L$ or
$k\ast\partial_{t}+L$; we refer to [63, 68, 69] for details.

###  3.2 Subdiffusion-Reaction Models

Modelling additional reactions and also the extension to systems is a subtle
task that can lead to rather different results. Most approaches come as
subdiffusion-limited (S) or activation-limited (A) reactions.

##### S: Subdiffusion limited.

In the more straightforward and broadly studied class (S), either the non-
local operator acts on the diffusion terms as well as the reaction terms [70,
71, 72, 73, 74, 75], or one simply adds an extra source or sink to the right
hand side of (12) [76], so that in the scalar case we obtain

$\displaystyle\partial_{t}u=\mathfrak{D}_{0,t}^{1-\alpha}[\sigma\partial_{x}^{2}u+f(u)]\quad\text{or}\quad\tilde{\mathfrak{D}}_{0,t}^{\alpha}u=\sigma\partial_{x}^{2}u+f(u),$
(18)

with $f:\mathbb{R}\to\mathbb{R}$. Here the reaction is “anomalous” even
without subdiffusive transport $\sigma=0$ and for spatially uniform states.
Analogously, adding a nonlinear source term to (17) gives rise to a time
fractional evolution equation

$k\ast\partial_{t}u=\sigma\partial_{x}^{2}u+f(u).$ (19)

The theory for Cauchy-problems in (18) is well developed (mostly on bounded
domains) and an excellent introduction is provided by [52], for which we also
refer to the numerous references in this active field. The initial-boundary
value problem of (19) and even more general semilinear diffusion cases have
been considered in [76, 77], which provide well-posedness and stability
theory. The existence problem of steady states in (18) and (19) is exactly the
same as for classical reaction-diffusion systems, although the stability
properties may differ. For instance, the linear stability principle has a
different character, and exponential decay in the stable case is replaced by
algebraic decay with rate $-\alpha$. Indeed, (18) for linear $f(u)=bu$ and for
a Fourier mode $\tilde{q}$ such that $b<\sigma\tilde{q}^{2}$ has the same form
as the subdiffusion equation (11) for the Fourier mode
$q=\sqrt{\tilde{q}^{2}-b/\sigma}$, consistent with the power $-\alpha$ in
(14). There are numerous studies of asymptotic stability for equations of the
form (19); we refer to [76, 77] as abstract studies and the references
therein. In addition, some recent studies concern the asymptotic stability for
modified forms of (19), e.g., with fractional Laplacian [78] or nonlinearity
involving infinite delays [79]. Another difference in stability properties is
that spatially constant steady states tend to be spectrally more stable
against non-steady perturbations in the presence of subdiffusion of class (S)
than classical diffusion, e.g. [58, 59, 80]. We refer to [59] for a recent
account of stability in fractional ODE.

Both aspects, i.e., non-expoential decay and (de)stabilizing effect, also
arise for the less considered class (A) on which we focus in the following.
Here the reaction kinetics is not affected by subdiffusive transport, i.e.,
any nonlocal time-dependent operator acts on the diffusion terms only [81, 70,
82, 83, 75]. Within this class of models we distinguish the two cases of added
source/sink (Aadd) and linear creation/annihilation (Alin) as detailed next.

##### Aadd: subdiffusion with added source/sink.

Akin to classical reaction-diffusion models, following [81], in this case one
simply adds an extra source or sink to the right-hand side of (11). More
generally, for a possibly nonlinear reaction term $f(u)$ one thus obtains

$\partial_{t}u=\sigma\mathfrak{D}_{0,t}^{1-\alpha}\partial_{x}^{2}u+f(u).$
(20)

While for the classical diffusion at $\alpha=1$ the usual maximum principle
holds, the Green’s function for (20) with negative linear reaction dynamics,
e.g., $f(u)=-u$, partially takes negative values, cf. [70]. Heuristically, the
reason is that the sink can remove particles much faster than they can be
transported from other positions. Hence, (20) is not a model for an absolute
density $u$, see also [84]. However, when $u$ represents a density
perturbation from a saturated strictly positive state this phenomenon is
acceptable, at least for short times. Within (20) this would concern $u\approx
u_{*}$, where $u_{*}>0$ solves $f(u_{*})=0$. A corresponding ad hoc model, cf.
[85], for multiple species is

$\partial_{t}{U}=\Sigma\mathfrak{D}_{0,t}^{1-\alpha}\partial_{x}^{2}{U}+{F}({U}),$
(21)

where ${U}={U}(t,x)\in\mathbb{R}^{N}$ is the vector of the perturbation
densities,
$\Sigma=\mathrm{diag}(\sigma_{1},\ldots,\sigma_{N})\in\mathbb{R}^{N\times N}$
is the diffusion matrix with positive entries, and
${F}:\mathbb{R}^{N}\to\mathbb{R}^{N}$ is a possibly nonlinear function that
models the reaction kinetics.

##### Alin: subdiffusion with linear creation/annihilation.

In this case, the addition or removal of particles through (linear) reactions
occurs only during the waiting time [70] and the resulting scalar equation
with reaction rate $a$ reads

$\partial_{t}u=\sigma
e^{at}\mathfrak{D}_{0,t}^{1-\alpha}(e^{-at}\partial_{x}^{2}u)+au.$

From the corresponding integral form one can infer that this model preserves
positivity of solutions since the amount of removed particles is less than the
amount of existing particles. For multiple species and reaction rate matrix
${A}\in\mathbb{R}^{N\times N}$ the model

$\partial_{t}{U}=\Sigma
e^{{A}t}\mathfrak{D}_{0,t}^{1-\alpha}(e^{-{A}t}\partial_{x}^{2}{U})+{A}{U},$
(22)

has been proposed in [86, 84], and rigorously derived in [87]. Since the
matrix exponentials are suitable for linear reactions only, the form of an
appropriate model for nonlinear reactions is not clear. We also refer to [84]
for a discussion of the modelling of reaction subdiffusion equations.

###  3.3 Spectra, Growth and Decay

The Fourier-Laplace transform can be employed also for systems to derive
dispersion relations of temporal and spatial modes, analogous to the classical
reaction diffusion systems. We illustrate decay and growth in subdiffusion
from this viewpoint. In contrast to classical diffusion, the fractional powers
in the dispersion relation yield branch points, and we already alluded to the
related algebraic decay in the subdiffusion equation in Section 3.1. Also
generally, the branch points strongly impact the decay laws determined by an
inverse Laplace transform, and moreover the decay properties of models (Aadd)
and (Alin) differ significantly, although the models are equal in absence of
reactions. In addition, the solutions to the dispersion relation to the left
of the most unstable branch point in the complex plane depend on the choice of
branch cut, and following [64] we therefore refer to these solutions as
pseudo-spectrum; when there is no ambiguity in this section we refer to all
solutions of dispersion relations (in the admissible domain) as pseudo-
spectrum. Non-canonical branch cuts admit to relate pseudo-spectrum to regular
spectrum of the classical diffusion equations through the limit $\alpha\to 1$.
For illustration, we consider two-component systems and start with classical
diffusion, i.e., a system

$\partial_{t}{U}=\Sigma\partial_{x}^{2}{U}+{F}({U}),$ (23)

where ${U}={U}(t,x)\in\mathbb{R}^{2}$ is the vector of density of the species
with $x\in\mathbb{R}$, $t>0$, ${F}:\mathbb{R}^{2}\to\mathbb{R}^{2}$ represents
the reaction kinetics, and $\Sigma=\mathrm{diag}(1,\sigma)$ is the diagonal
matrix of positive constant diffusion coefficients. Here the first diffusion
coefficient is normalized to $1$ by a suitable scaling of $x$ so that $\sigma$
is in fact the ratio of diffusion coefficients. Suppose that (23) possesses a
spatially homogeneous steady state ${U}_{*}\in\mathbb{R}^{2}$, i.e.,
${F}({U}_{*})=0$. The linearization of (23) about ${U}_{*}$ is given by

$\partial_{t}{U}=\Sigma\partial_{x}^{2}{U}+{A}{U},\quad{A}:=\left({\textnormal{D}}_{U}{F}({U}_{*})\right)=\begin{pmatrix}a_{1}&a_{2}\\\
a_{3}&a_{4}\end{pmatrix}.$ (24)

The right-hand side of (24) defines the linear operator
$\mathcal{L}:=\Sigma\partial_{x}^{2}+{A}$, with eigenvalue problem

$\mathcal{L}{U}=\Sigma\partial_{x}^{2}{U}+{A}{U}=s{U}.$

Upon Fourier transforming $x\in\mathbb{R}$ to Fourier-wavenumbers
$q\in\mathbb{R}$, the eigenvalue problem becomes

$\mathcal{L}(\mathcal{F}{U})=-q^{2}\Sigma(\mathcal{F}{U})+{A}(\mathcal{F}{U})=s(\mathcal{F}{U}),$

whose solvability condition as a matrix eigenvalue problem yields the
dispersion relation

$D_{{\rm reg}}(s,q^{2}):=\det(s{I}+q^{2}\Sigma-{A})=(s+q^{2}-a_{1})(s+\sigma
q^{2}-a_{4})-a_{2}a_{3}=0.$ (25)

The solution set $\Lambda_{\rm reg}:=\\{s\in\mathbb{C}:D_{{\rm
reg}}(s,q^{2})=0\text{ for a }q\in\mathbb{R}\\}$ defines in particular the
$\mathnormal{L}^{2}$-spectrum of $\mathcal{L}$ with domain
$(\mathnormal{H}^{2}(\mathbb{R}))^{2}$, e.g., [88], which we refer to as
regular spectrum. The homogeneous steady state ${U}_{*}$ is called strictly
spectrally stable (unstable) if $\max({\textnormal{Re}}(\Lambda_{\rm reg}))<0$
($>0$). It is then also linearly and nonlinearly stable (unstable) for (24)
and (23), respectively [88].

##### Dispersion relation for Aadd:

Proceeding in the same way for (21) gives

$(s{I}+s^{1-\alpha}q^{2}\Sigma-{A})\,\mathfrak{L}\mathcal{F}{U}=\mathcal{F}{U}_{0},\quad{\textnormal{Re}}(s)>0,\
q\in\mathbb{R},$

with the constraint on $s$ due to the branch point of $s^{1-\alpha}$ at the
origin. Hence, the analogue of (25) is the dispersion relation on
$\Omega_{0}^{+}:=\\{s\in\mathbb{C}\setminus\\{0\\}:\arg(s)\in(-\pi/2,\pi/2)\\}$
given by

$\displaystyle D_{\rm add}(s,q^{2})$
$\displaystyle:=\det(s{I}+s^{1-\alpha}q^{2}\Sigma-{A})=(s+s^{1-\alpha}q^{2}-a_{1})(s+s^{1-\alpha}\sigma
q^{2}-a_{4})-a_{2}a_{3}=0.$

For any given $q\in\mathbb{R}$ we can extend the domain of $D_{\rm add}$ on
the principal branch such that the values on positive reals are preserved. Due
to the discrete solution set, we can choose $\theta_{1}(q)\in(0,\pi/2)$ such
that there is no solution of the dispersion relation on the branch cut

$\mathfrak{B}_{0}^{\theta_{1}(q)}:=\\{s\in\mathbb{C}:{\textnormal{Im}}(s)/{\textnormal{Re}}(s)=\tan(\theta_{1}(q)),\,{\textnormal{Re}}(s)<0\\},$

i.e., $D_{\rm add}(s,q^{2})\neq 0$ for $s\in\mathfrak{B}_{0}^{\theta_{1}(q)}$.
For simplicity we suppress the dependence of $\theta_{1}$ on $q$, and obtain
the extended domain of $D_{\rm add}$ as

$\Omega_{0}:=\\{s\in\mathbb{C}\setminus\\{0\\}:\arg(s)\in(-\pi+\theta_{1},\pi+\theta_{1}),\,\theta_{1}\in(0,\pi/2)\\}.$

The pseudo-spectrum $\Lambda_{\rm add}$ is the corresponding set of roots of
$D_{\rm add}$ on $\Omega_{0}$.

##### Dispersion relation for Alin:

For (22) we follow [86, 64, 87] and make the simplifying assumption that ${A}$
is diagonalizable, so that with the eigenvalues $\mu_{1},\mu_{2}$ of ${A}$,
the coordinate change

${P}=\begin{pmatrix}a_{2}&a_{2}\\\
\mu_{1}-a_{1}&\mu_{2}-a_{1}\end{pmatrix},\quad{P}^{-1}=\begin{pmatrix}\frac{a_{1}-\mu_{2}}{a_{2}(\mu_{1}-\mu_{2})}&\frac{1}{\mu_{1}-\mu_{2}}\\\
-\frac{a_{1}-\mu_{1}}{a_{2}(\mu_{1}-\mu_{2})}&-\frac{1}{\mu_{1}-\mu_{2}}\end{pmatrix}$

gives $\bar{{A}}:={P}^{-1}{A}{P}=\mathrm{diag}(\mu_{1},\mu_{2})$. We order the
eigenvalues so that
${\textnormal{Re}}(\mu_{1})\geq{\textnormal{Re}}(\mu_{2})$. Setting
${U}={P}{W}$ turns (22) into

$\partial_{t}{W}=\bar{\Sigma}e^{\bar{{A}}t}\mathfrak{D}_{0,t}^{1-\alpha}(e^{-\bar{{A}}t}\partial_{x}^{2}{W})+\bar{{A}}{W},\quad\bar{\Sigma}=\begin{pmatrix}d_{1}&d_{2}\\\
d_{3}&d_{4}\end{pmatrix}:={P}^{-1}\Sigma{P}.$

Fourier-Laplace transforming then gives the dispersion relation

$\begin{split}D_{\rm lin}(s,q^{2}):=\
&\det(s{I}-\bar{{A}}+(s{I}-\bar{{A}})^{1-\alpha}\bar{\Sigma}q^{2})=\det((s{I}-\bar{{A}})^{1-\alpha})\det((s{I}-\bar{{A}})^{\alpha}+\bar{\Sigma}q^{2})\\\
=\
&(s-\mu_{1})^{1-\alpha}(s-\mu_{2})^{1-\alpha}\big{(}((s-\mu_{1})^{\alpha}+d_{1}q^{2})((s-\mu_{2})^{\alpha}+d_{4}q^{2})-d_{2}d_{3}q^{4}\big{)}=0,\end{split}$

where ${\textnormal{Re}}(s-\mu_{1}),\,{\textnormal{Re}}(s-\mu_{2})>0$. As
above, we may extend the domain to $s\in\Omega_{\rm lin}\subset\mathbb{C}$ by
choosing branch cuts that connect $\infty$ to the finite branch points
$\mu_{1}$ and $\mu_{2}$, respectively. It turns out, cf. [64], that $s=\infty$
is a branch point of $(s-\mu_{1})^{\alpha}(s-\mu_{2})^{\alpha}$ for
$\alpha\in(0,1/2)\cup(1/2,1)$; for $\alpha=1/2$ we can think of $\infty$ as an
intermediate point on the branch cut connecting $\mu_{1}$ and $\mu_{2}$. For
simplicity, we assume real $\mu_{1}>\mu_{2}$; the general case can be found in
[64]. In order to keep non-integer powers of positive reals on the positive
real axis, we then choose branch cuts $\mathfrak{B}_{\mu_{j}}^{\theta_{j}}$,
$j=1,2$, such that

$\Omega_{\rm
lin}:=\\{s\in\mathbb{C}\setminus\\{\mu_{1},\mu_{2}\\}:\arg(s-\mu_{j})\in(-\pi+\theta_{j},\pi+\theta_{j}),\
\theta_{j}\in(0,\pi/2),\ j=1,2\\}.$

As for (Aadd), given $q\in\mathbb{R}$ we select $\theta_{1},\theta_{2}$ such
that the solutions of $D_{\rm lin}(\cdot,q^{2})=0$ do not lie on the branch
cuts. The pseudo-spectrum $\Lambda_{\rm lin}$ is the corresponding set of
roots of $D_{\rm lin}$ on $\Omega_{\rm lin}$.

######  Remark.

We illustrate the special nature of rational $1-\alpha=n/m\in\mathbb{Q}$ for
the non-trivial factor of $D_{\rm lin}$,

$\displaystyle D_{{\rm lin}2}(s,q^{2})$
$\displaystyle:=\big{(}(s-\mu_{1})^{\alpha}+d_{1}q^{2}\big{)}\big{(}(s-\mu_{2})^{\alpha}+d_{4}q^{2}\big{)}-d_{2}d_{3}q^{4}=0,\quad
s\in\Omega_{\rm lin}.$

In the rational case it can be recast as a polynomial with respect to the two
variables $z_{j}=(s-\mu_{j})^{1/m}$, $j=1,2$, of degree $m-n$ in each
variable, given by

$\tilde{D}_{{\rm
lin}2}(z_{1},z_{2},q^{2})=(z_{1}^{m-n}+d_{1}q^{2})(z_{2}^{m-n}+d_{4}q^{2})-d_{2}d_{3}q^{4}.$

By resultant theory, the polynomial system $\tilde{D}_{{\rm
lin}2}(z_{1},z_{2},q^{2})=0$, $z_{1}^{m}+\mu_{1}=z_{2}^{m}+\mu_{2}$ has a
finite number of roots $(z_{1},z_{2})$, which gives a finite number of roots
$s$ for $D_{{\rm lin}2}(s,q^{2})$, including multiplicities.

##### Growth and decay:

Via an inverse Laplace transform, the location of pseudo-spectrum translates
into growth and decay properties in the time-domain for each Fourier-
wavenumber. For case (Aadd), we quote a simplified form of [64, Theorem 4.2],
where the following set of solutions with maximal real part is relevant:

$S^{+}:=\\{(s,q)\in\overline{\Omega_{0}^{+}}\times\mathbb{R}:D_{\rm
add}(s,q^{2})=0\text{ and }{\textnormal{Re}}(s)\mbox{ is maximal for }q\\}.$

######  Theorem 3.1 (Case Aadd).

Let $\alpha=\tfrac{n}{m}\in(0,1),\,n,m\in\mathbb{N}$, be a reduced fraction,
$\lambda:=\sup({\textnormal{Re}}(\Lambda_{\rm add}))$, and
$\mathcal{F}{U}(t,q)$ a solution to the Fourier transformed linearization of
(21) in $U_{*}$ with $\mathcal{F}{U}(0,q)=\mathcal{F}{U}_{0}(q)$.

  * (1)

If $\lambda\geq 0$ then for any $(s_{0},q_{0})\in S^{+}$ we have
$\mathcal{F}{U}(t,q_{0})=C(q_{0},\mathcal{F}{U}_{0}(q_{0}))t^{k-1}e^{s_{0}t}+o(t^{k-1}e^{{\textnormal{Re}}(s_{0})t})$
with $C(q_{0},\mathcal{F}{U}_{0}(q_{0}))\in\mathbb{C}^{2}$ nonzero for almost
all $\mathcal{F}{U}_{0}$ and $k$ the multiplicity of $s_{0}$ as a root of
$D_{\rm add}(\cdot,q_{0}^{2})=0$.

  * (2)

Assume $\det({A})\neq 0$. If (i) $\lambda=0$ and $S^{+}=\emptyset$ or (ii)
$\lambda<0$ or (iii) $\Lambda_{\rm add}=\emptyset$, then for any $q\neq 0$
there is a $q$-independent invertible ${M}\in\mathbb{R}^{2\times 2}$ such that
$\mathcal{F}{U}(t,q)=q^{2}{M}\mathcal{F}{U}_{0}(q)t^{\alpha-2}+\mathcal{O}(t^{\alpha-2-1/m})$.
In addition, the entries of ${M}$ approach $0$ as $\alpha\to 1$ and $\infty$
as $\det({A})\to 0$.

Concerning item (1) we note that $S^{+}\neq\emptyset$ if $\lambda>0$. In item
(2) the condition $\det({A})\neq 0$ has been added compared to [64, Theorem
4.2], where this was assumed a priori, cf. Section 3.4. In case $q=0$ we
obtain $\partial_{t}(\mathcal{F}{U})={A}(\mathcal{F}{U})$, whose solutions are
either constant, grow algebraically as in item (1) for $s_{0}=0$, or decay
exponentially. Theorem 3.1 means that the growth of Fourier modes is the same
as for classical diffusion, $\alpha=1$, but the decay properties differ: in
the stable case each Fourier mode decays algebraically as $t^{\alpha-2}$.
Notably, this is in contrast to the decay law $t^{-\alpha}$ for the
subdiffusion-limited class (S) discussed in Section 3.2, see also [58, 74,
59]. For the case (Alin) with $\mu_{1}>\mu_{2}$ we quote a simplified form of
[64, Theorem 5.3]. Here we denote $\Omega_{\rm lin}^{-}:=\\{s\in\Omega_{\rm
lin}:{\textnormal{Re}}(s)<{\textnormal{Re}}(\mu_{1})\\}$, $\Omega_{\rm
lin}^{0+}:=\Omega_{\rm lin}\setminus\Omega_{\rm lin}^{-}$ and

$\displaystyle S^{+}$ $\displaystyle:=\\{(s,q)\in\Omega_{\rm
lin}^{0+}\times(\mathbb{R}\setminus\\{0\\}):D_{\rm lin}(s,q^{2})=0,\text{ and
}{\textnormal{Re}}(s)\text{ maximal for }q\\},$ $\displaystyle S^{-}$
$\displaystyle:=\\{(s,q)\in\Omega_{\rm
lin}^{-}\times(\mathbb{R}\setminus\\{0\\}):D_{\rm lin}(s,q^{2})=0,\text{ and
}{\textnormal{Re}}(s)\text{ maximal for }q\\}.$

######  Theorem 3.2 (Case Alin).

Assume $\mu_{1}>\mu_{2}$. Let
$\alpha=\tfrac{n}{m}\in(0,1),\,n,m\in\mathbb{N}$, be a reduced fraction,
$\lambda:=\sup({\textnormal{Re}}(\Lambda_{\rm lin}))$, and
$\mathcal{F}{U}(t,q)$ a solution to the Fourier transform of (22) with
$\mathcal{F}{U}(0,q)=\mathcal{F}{U}_{0}(q)$.

  * (1)

If $\lambda\geq\mu_{1}$, then for any $(s_{0},q_{0})\in S^{+}$ we have
$\mathcal{F}{U}(t,q_{0})=C(q_{0},\mathcal{F}{U}_{0}(q_{0}))t^{k-1}e^{s_{0}t}+o(t^{k-1}e^{{\textnormal{Re}}(s_{0})t})$
with $C(q_{0},\mathcal{F}{U}_{0}(q_{0}))\in\mathbb{C}^{2}$ nonzero for almost
all initial data and $k$ the multiplicity of $s_{0}$ as a root of $D_{{\rm
lin}2}(\cdot,q_{0}^{2})=0$.

  * (2)

If (i) $\lambda=\mu_{1}$ and $S^{+}=\emptyset$ or (ii)
$\lambda\in(\mu_{2},\mu_{1})$, then $S^{-}\neq\emptyset$ and for any
$(s_{0},q_{0})\in S^{-}$ there exist $q$-independent nonzero
$c_{j}\in\mathbb{C},\,j=1,2$, and a singular ${M}_{1}(q)\in\mathbb{C}^{2\times
2}$ with one-dimensional kernel such that for any initial data satisfying
${P}^{-1}\mathcal{F}{U}_{0}(q_{0})\notin\ker({M}_{1}(q_{0}))$ we have

$\displaystyle\mathcal{F}{U}(t,q_{0})=$ $\displaystyle\
q_{0}^{-2}{P}\begin{pmatrix}c_{1}&0\\\
0&c_{2}t^{-1}\end{pmatrix}{M}_{1}(q_{0}){P}^{-1}\mathcal{F}{U}_{0}(q_{0})\,e^{\mu_{1}t}t^{-\alpha}+\mathcal{O}(t^{-\alpha-1/m}e^{\mu_{1}t}).$

  * (3)

If (i) $\lambda\leq\mu_{2}$ or (ii) $\Lambda_{\rm lin}=\emptyset$, then for
any $q\neq 0$ and any initial data satisfying
${P}^{-1}\mathcal{F}{U}_{0}(q)\notin\ker({M}_{1}(q))$, with $c_{j},{M}_{1}(q)$
as in item (2), it holds that

$\displaystyle\mathcal{F}{U}(t,q)=$ $\displaystyle\
q^{-2}{P}\begin{pmatrix}c_{1}&0\\\
0&c_{2}t^{-1}\end{pmatrix}{M}_{1}(q){P}^{-1}\mathcal{F}{U}_{0}(q)\,e^{\mu_{1}t}t^{-\alpha}+\mathcal{O}(t^{-\alpha-1/m}e^{\mu_{1}t}).$

In addition, $c_{j}\to 0$ as $\alpha\to 1$ for $j=1,2$, and the entries of
${M}_{1}(q)$ are bounded for any $q\neq 0$.

Concerning item (1), if $\lambda>\mu_{1}$, then $S^{+}$ is guaranteed to be
non-empty. Items (2,3) can be understood as thresholds for decay rates. In
particular, exponential decay for any Fourier mode cannot be faster than that
of the kinetics. For initial data in the kernel of ${M}_{1}(q)$, the leading
order terms of the corresponding solution is rather involved; we refer to [64,
Theorem 5.3] for further details. The statement also holds for $\mu_{1}=0$ and
scalar equations, and thus includes the subdiffusion equation (11). The decay
law $t^{-\alpha}$ (for rational $\alpha$) indeed coincides with that derived
in Section 3.1 via (14), and with that for (18), e.g., in [58, 68]. As for
Theorem 3.1 the method of proof only allows for rational $\alpha$. In contrast
to Theorem 3.1, the coefficients decay as $q^{-2}$ in Theorem 3.2 so that some
smoothing of initial data is conceivable.

In summary, a pseudo-spectrum with positive real parts readily implies
exponential instability with the same spectral rate that would arise in
classical diffusion with the same real parts. Negative real parts, however, do
not imply the exponential decay of the classical case: In case (Aadd) there is
only algebraic decay (in systems at least for rational $\alpha$) of power
$\alpha-2$, and indeed the origin is always a branch point of the dispersion
relation. In case (Alin) the branch points lie at the eigenvalues of ${A}$,
which admits exponential decay if the pseudo-spectra have negative real parts.
However, the decay cannot be faster than the rate of the most unstable
eigenvalue of ${A}$, and there is an algebraic correction with power
$-\alpha$. From a technical point of view, the essence of two distinct decays
is that in case (Aadd) the branch point does not solve the dispersion
relation, while it does in case (Alin). As mentioned, also in model class (S)
strictly stable spectrum is expected to give only algebraic decay, which is
proven by energy methods for Dirichlet boundary conditions (together with
further results) in [68, 76]; we also refer to Section 4 for more details on
global/gradient-flow methods based upon energy dissipation functionals.

In the results for subdiffusion models of class (A), $\alpha$ is taken
rational since then the inverse Laplace transform can be computed explicitly
via residues. Although the results suggest uniformity, there seems to be no
proof. Moreover, the estimates for a selected Fourier wavenumber, and although
naturally locally uniform and with asymptotic rates, there appears to be no
additional results in $x$-space. Both issues may be approachable by energy
methods. For the case (Aadd), i.e., (21), constants grow as $q^{2}$, thus
requiring sufficient smoothness in $x$-space. In contrast, in the case (Alin),
i.e., (22), the analogous coefficients decay as $q^{-2}$, which indicates some
smoothing. However, results on well-posedness for these models seem not to be
available in the literature.

###  3.4 Turing-type instability

Concerning the onset of instability, models (21) and (22) differ significantly
from each other and from the case of normal diffusion $\alpha=1$. We quote
some results for which we refer to [64] and the references therein, in
particular [83, 75, 85, 89]. These concern the case of two species systems,
$N=2$, and when parameters are such that $\alpha=1$ admits a so-called Turing
instability for a critical diffusion ratio $\sigma>\sigma_{\rm reg}$. For this
instability the onset at $\sigma=\sigma_{\rm reg}$ occurs with a unique finite
wavenumber $q=q_{\rm reg}$, while ${A}$ has eigenvalues with strictly negative
real parts, i.e., $\mathrm{tr}({A})=a_{1}+a_{4}<0$ and
$\det({A})=a_{1}a_{4}-a_{2}a_{3}>0$. Without loss of generality we may assume
$a_{1}>0$, $a_{4}<0$, $a_{2}a_{3}<0$, which gives $\sigma_{\rm reg}>1$ via the
well known instability condition $\sigma a_{1}+a_{4}>0$, e.g., [3]. Turing
instability for class (S) is well studied: roughly speaking, the subdiffusion
does not have an effect on the steady instability threshold; we refer to [74,
75, 80] for further details. However, this is not the case for models of class
(A).

##### Aadd: destabilizing subdiffusion.

For (21), if the spectrum for $\alpha=1$ is Turing unstable, then the spectrum
for all $\alpha\in(0,1)$ is unstable. Somewhat counterintuitively, this means
that, in terms of the diffusion ratio $\sigma$, the reaction-subdiffusion in
case (Aadd) is always less stable than classical diffusion. In particular, the
diffusion instability threshold $\sigma_{\text{add}}$ for (Aadd) is smaller
than that for classical diffusion. A more significant deviation from parabolic
equations is that the spectrum becomes unstable through infinite wavenumbers,
in the following sense: for any $\sigma$ the grow rate at $q=\pm\infty$ is
zero, and for $\sigma>\sigma_{\text{add}}$ the unstable modes have wavenumber
$q$ with $|q|\in[q(\sigma),\infty)$, where $q(\sigma)\to\infty$ as
$\sigma\searrow\sigma_{\text{add}}$. In particular, the signature finite
wavenumber selection at the onset of a Turing instability is absent. At the
same time, in the limit $\alpha\to 1$ the pseudo-spectrum converges locally
uniformly as a set in $\mathbb{C}$ to the regular spectrum. Hence, there is
also an instability threshold in terms of $\alpha$. We refer to [64] for some
numerical explorations.

##### Alin: stabilizing subdiffusion.

For (22) the analytical comparison of regular and pseudo-spectrum is
constrained to $s\approx 0$. In this region, the pseudo-spectrum is real and
more stable than normal diffusion, and pseudo-spectrum becomes unstable via
finite wavenumber at a diffusion ratio $\sigma_{\text{lin}}>\sigma_{\rm reg}$,
if $\alpha\geq\alpha_{A}$ with a certain threshold $\alpha_{A}$. For
$\alpha<\alpha_{A}$ the instability is completely suppressed. In the limit
$\alpha\to 1$ the pseudo-spectrum again converges locally uniformly as a set
in $\mathbb{C}$ to the regular spectrum. This means that for $\alpha\approx 1$
the critical pseudo-spectrum is real and therefore it is globally more stable
than regular spectrum. Numerical computations suggest that this relation holds
broadly in parameter space, cf. [64].

##  4 Gradient Flows: Time Fractional PDEs

In the last section, we have seen how to analyze dissipative subdiffusion
equations locally by focusing on spectra and bifurcations. In this section, we
describe the natural complementary view on this class of differential
equations via global gradient flow techniques. Gradient flows have been
vigorously studied in the field of nonlinear PDEs. An abstract form of
gradient flows in an appropriate space ${\mathcal{X}}$ (e.g., Hilbert and
Banach spaces) reads,

$\partial_{t}u(t)=-\nabla E(u(t))\ \mbox{ for }\ t>0,\quad u(0)=u_{0},$ (26)

where $\partial_{t}=\text{\rm d}/\text{\rm d}t$,
$u:[0,\infty)\to{\mathcal{X}}$ is a gradient flow emanating from an initial
datum $u_{0}$ and $\nabla E$ denotes a functional derivative in a proper sense
(e.g., Fréchet derivative and subdifferential) of an energy (or entropy)
functional $E:{\mathcal{X}}\to\bar{\mathbb{R}}:=[-\infty,\infty]$. There are
many parabolic equations (systems), e.g., Cahn-Hilliard system, Allen-Cahn
equation, nonlinear diffusion equations such as porous medium and fast
diffusion equations, many of phase transition models, which admit variational
structures of gradient flow type. One of crucial features of gradient flows is
the decrease (in time) of the energy $t\mapsto E(u(t))$. Moreover, any
elements of the $\omega$-limit set of each gradient flow are equilibria, i.e.,
critical points of the functional $E$. On the other hand, the non-emptiness of
the $\omega$-limit set is nontrivial, and moreover, even if it is nonempty, it
is rather delicate whether or not the $\omega$-limit set is singleton.

This section is concerned with a fractional gradient flow, which is a variant
of gradient flows with time-fractional derivatives instead of the classical
one, say

$\partial_{t}^{\alpha}(u-u_{0})(t)=-\nabla E(u(t))\ \mbox{ for }\ t>0,$ (27)

where $\alpha\in(0,1)$ and
$\partial_{t}^{\alpha}:=\mathfrak{D}_{0,t}^{\alpha}$ is a shorthand notation
for the Riemann-Liouville derivative throughout this section.

###  4.1 Brief History

The study of fractional gradient flows may date back to the 1970s, when the
following nonlinear Volterra equation was studied extensively:

$u(t)+(a*\xi)(t)=g(t),\quad\xi(t)\in A(u(t)),\quad t>0,$ (28)

where $A:{\mathcal{X}}\to{\mathcal{X}}$ is a (possibly multivalued) nonlinear
operator in a Hilbert or Banach space ${\mathcal{X}}$,
$g:(0,T)\to{\mathcal{X}}$ is given and $a:[0,+\infty)\to[0,+\infty]$ is a
kernel function of an appropriate class. Actually, the convolution of both
sides of (27) with the conjugate kernel $k_{\alpha}$ yields (28) with the
choice $a=k_{\alpha}$, $A=-\nabla E$ and $g(t)\equiv u_{0}$, since
$k_{\alpha}*k_{1-\alpha}$ is identically equal to $1$. In those days, equation
(28) was often studied in a spirit of the nonlinear semigroup theory (see,
e.g., [90, 91, 92]).

Barbu [93, 94, 95] studied an abstract nonlinear Volterra equation (28), which
arises in the study of mechanical systems with memory effects, under the
assumptions that $a$ is of class $W^{1,1}_{\rm loc}([0,+\infty))$ (hence it is
differentiable almost everywhere and finite at $t=0$) and positive, $g\in
W^{1,2}(0,T;{\mathcal{X}})$ and $A=\partial\varphi$ in a Hilbert space
${\mathcal{X}}$. Differentiating both sides of (28) in time, we see that

$u^{\prime}(t)+a(0)\xi(t)+\big{(}a^{\prime}*\xi\big{)}(t)=g^{\prime}(t),\quad\xi(t)\in
A(u(t)),\quad 0<t<T,$

which can also be regarded as a nonlocal (in time) perturbation problem of a
(local) nonlinear evolution equation. We also refer the reader to [96, 97, 98,
99, 100, 101]. However, these results cannot be applied to the Riemann-
Liouville kernel $k_{\beta}$, for $k_{\beta}^{\prime}(0+)$ is always infinite.
As for singular kernels such as the Riemann-Liouville kernel $k_{\beta}$,
Clément and Nohel [102, Theorem 3.1] studied (28) for _completely positive
kernels_ $a\in L^{1}_{\rm loc}([0,+\infty))$ and proved existence of a
generalized solution, which is a limit of a certain class of approximate
solutions. Moreover, the abstract result is also applied to a couple of
nonlocal (in time) PDEs (see also [103]). We further refer the reader to [104]
and [105, Theorems 2 and 3], which are closely related to the main focus of
this section. On the other hand, in these works, no application to the time-
fractional problem is explicitly mentioned. Evolution equations including
nonlocal derivatives (e.g., Riemann-Liouville derivative) are studied in [106,
107, 108] by finding out that nonlocal differential operators are
$m$-accretive in Bochner spaces (see Section 4.2 below). On the other hand,
most existence results are established for generalized solutions and, to the
best of our knowledge, there had been no result which can guarantee the
existence of strong (in time) solutions for time-fractional abstract evolution
equations of gradient flow type prior to Zacher’s pioneer work [109], which is
still concerned with only some linear evolution equations and will be more
detailed below. One of the reasons for such a stagnation might reside in the
lack of useful chain-rule formulae, which play an indispensable role in the
analysis of (classical) gradient flows (26).

Let $V$ and $V^{*}$ be a Hilbert space and its dual space, respectively, and
let $H$ be a pivot Hilbert space such that the Gelfand triplet holds,

$V\hookrightarrow H\simeq H^{*}\hookrightarrow V^{*}.$

Based on the Galerkin approximation, Zacher [109] proved existence and
uniqueness of strong (in time) solutions to the linear evolution equation,

$\dfrac{\text{\rm d}}{\text{\rm
d}t}\big{(}\left[k*(u-u_{0})\right](t),v\big{)}_{H}+b(t,u(t),v)=\langle
f(t),v\rangle_{V}\ \mbox{ for all }v\in V,$ (29)

where $(\cdot,\cdot)_{H}$ and $\langle\cdot,\cdot\rangle_{V}$ stand for an
inner product of $H$ and a duality pairing between $V$ and $V^{*}$,
respectively, $k$ is a completely positive kernel, $f:(0,T)\to V^{*}$ and
$u_{0}\in H$ are given and $b(t,\cdot,\cdot)$ is a (time-dependent) bounded
coercive bilinear form defined on $V$. Key devices for proving the existence
of strong (in time) solutions to (29) consist not only of the $m$-accretivity
of the nonlocal time-differential operator but also of a nonlocal energy
identity (which is a sort of the chain-rule formula),

$\displaystyle{\left(\dfrac{\text{\rm d}}{\text{\rm
d}t}\left(h*u\right)(t),u(t)\right)_{H}=\dfrac{1}{2}\dfrac{\text{\rm
d}}{\text{\rm
d}t}\left(h*\|u\|_{H}^{2}\right)(t)+\dfrac{1}{2}h(t)\|u(t)\|_{H}^{2}}$
$\displaystyle-\dfrac{1}{2}\int^{t}_{0}h^{\prime}(s)\|u(t)-u(t-s)\|_{H}^{2}\,\text{\rm
d}s\quad\mbox{ for a.e. }t\in(0,T)$ (30)

for any $h\in W^{1,1}(0,T)$ and $u\in L^{2}(0,T;H)$ (cf. Proposition 4.3
below). Moreover, in [110], convergence of nonlocal (in time) gradient flows
in Euclidean spaces $\mathbb{R}^{d}$ to equilibria is proved under certain
assumptions, which seem to exclude Riemann-Liouville and Caputo derivatives
but can treat their regularizations. Let us also briefly review other
contributions of R. Zacher and his collaborators (e.g., V. Vergara) without
any claim of completeness. In [111], boundedness (both in space and time) of
solutions to the time-fractional PDE,

$\partial_{t}^{\alpha}(u-u_{0})-\mathrm{div}\left(A(x,t)\nabla
u\right)=g(x,t),$ (31)

where $A(x,t)$ is a uniformly elliptic but possibly discontinuous matrix field
and $u_{0}$ and $g$ are prescribed, as well as its quasilinear variants (e.g.,
$A$ also depends on $u$) is proved by means of De Giorgi’s iteration technique
(see also [112]). This result was extended to degenerate and singular
diffusion equations of $p$-Laplacian type in [113], and moreover, a weak
Harnack inequality is also proved for (31) in [114]. Furthermore, an optimal
decay estimate is established for some time-fractional diffusion equations in
[68], where a chain-rule formula (cf. Proposition 4.3) as well as an
asymptotic analysis of solutions to time-fractional ODEs play a crucial role
(see also [115]). On the other hand, existence of solutions was assumed and is
still open as a question in these papers. We shall discuss the existence of
solutions in Section 4.3. In [76], the method developed in [68] is also
extended to a finite-time blow-up analysis of solutions to time-fractional
semilinear diffusion equation. Furthermore, a series of results on the time-
fractional (linear) diffusion equation is conveniently summarized in [116],
where an optimal decay estimate (see also [69]) and an $L^{p}$-maximal
regularity result (see also [117, 118]) are particularly interesting.

We close this section by pointing the reader to two monographs on time-
fractional evolution equations. A theory based on an integral equation, which
is standard for the case $\alpha=1$, rather than energy methods is well
established by Gal and Warma in [119], where various time-fractional
semilinear parabolic problems are comprehensively studied. Furthermore, we
also quote [120] as a nice summary of abstract differential equations with
fractional time-derivatives. In particular, a novel Sobolev framework for
fractional time-derivatives is developed, and moreover, decay properties of
solutions to time-fractional PDEs are discussed.

###  4.2 Abstract Theory for Fractional Gradient Flows

In this subsection, we consider a fractional gradient flow of subdifferential
type in a Hilbert space. To this end, we set up a couple of notions: Let $H$
be a real Hilbert space and let $\phi:H\to(-\infty,+\infty]$ be a proper
(i.e., $\phi\not\equiv+\infty$) lower-semicontinuous convex functional with
_effective domain_ $D(\phi):=\\{w\in H\colon\phi(w)<+\infty\\}$. Then the
subdifferential operator $\partial\phi:H\to 2^{H}$ (here $2^{H}$ denotes the
power set of $H$) is defined by

$\displaystyle\partial\phi(w)=\\{\xi\in
H\colon\phi(v)-\phi(w)\geq(\xi,v-w)_{H}\ \mbox{ for }v\in D(\phi)\\}\ \mbox{
for }w\in D(\phi)$

with domain $D(\partial\phi):=\\{w\in
D(\phi)\colon\partial\phi(w)\neq\emptyset\\}$. The subdifferential is a
generalized notion of the Fréchet derivative for convex functionals, and it is
necessary to establish an abstract theory which is applicable to PDEs having
gradient flow structures. The Cauchy problem for a fractional gradient flow
reads,

$-\partial_{t}^{\alpha}(u-u_{0})(t)\in\partial\phi(u(t))\ \mbox{ for }\ t>0,$
(32)

where $u_{0}\in D(\phi)$ is a prescribed initial datum. We recall the Riemann-
Liouville kernel

$k_{\beta}(t):=\frac{t^{\beta-1}}{\Gamma(\beta)}$

of order $\beta\in(0,1)$ and
$\partial_{t}^{\alpha}(u-u_{0})=\mathfrak{D}_{0,t}^{\alpha}(u-u_{0})=\partial_{t}[k_{1-\alpha}*(u-u_{0})]$
for $\alpha\in(0,1)$ (see also (4)). We are then concerned with strong
solutions to (32) defined as

######  Definition 4.1 (Strong solution).

For $0<T<\infty$, a function $u\in L^{2}(0,T;H)$ is called a strong solution
to (32) on $[0,T]$ if the following conditions are all satisfied:

  1. (i)

$k_{1-\alpha}*(u-u_{0})\in W^{1,2}(0,T;H)$ and
$[k_{1-\alpha}*(u-u_{0})](0)=0$,

  2. (ii)

$u(t)\in D(\partial\phi)$ for a.e. $t\in(0,T)$ and

$-\partial_{t}^{\alpha}(u-u_{0})(t)\in\partial\phi(u(t))$

for a.e. $t\in(0,T)$.

Moreover, a function $u\in L^{2}_{\rm loc}([0,\infty);H)$ is called a global
strong solution to (32) if for every $0<T<\infty$ the restriction $u|_{[0,T]}$
of $u$ to $[0,T]$ is a strong solution to (32) on $[0,T]$.

Define the space $H^{\alpha,2}_{0}(0,T;H)$ by

$H^{\alpha,2}_{0}(0,T;H)=\\{w|_{[0,T]}\colon w\in H^{\alpha,2}(\mathbb{R};H)\
\mbox{ and }\ \mathrm{supp}\,w\subset[0,\infty)\\},$

where $H^{\alpha,2}(\mathbb{R};H)\subset L^{2}(\mathbb{R};H)$ denotes the
Bessel potential function space of order $\alpha$. Then the main result of
this section is stated as follows:

######  Theorem 4.2 ([121, Theorem 2.3, Corollary 2.7]).

For every $u_{0}\in D(\phi)$, the Cauchy problem admits a unique global strong
solution $u\in L^{2}_{\rm loc}([0,\infty);H)$ such that the following holds
true:

  1. (i)

Set
$\mathcal{F}(t):=\left[k_{1-\alpha}*\left(\phi(u(\cdot))-\phi(u_{0})\right)\right](t)$
for $t\geq 0$. Then the function $t\mapsto\mathcal{F}(t)$ is continuous on
$[0,\infty)$, and it holds that

$\displaystyle\int^{t}_{s}\left\|\partial_{t}^{\alpha}(u-u_{0})(\tau)\right\|_{H}^{2}\,\text{\rm
d}\tau+\mathcal{F}(t)\leq\mathcal{F}(s)$

for all $0\leq s\leq t<\infty$.

  2. (ii)

It holds that $u\in H^{\alpha,2}_{0}(0,T;H)$ for $0<T<\infty$, and moreover,

$u\in\begin{cases}L^{\infty}(0,T;H)&\mbox{ if }\ \alpha\in(0,1/2],\\\
C([0,T];H)&\mbox{ if }\ \alpha>1/2\end{cases}$ (33)

for $0<T<\infty$. In particular, if $\alpha>1/2$, then the initial condition
$u(0)=u_{0}$ holds true in a classical sense.

  3. (iii)

It holds that $\phi(u(\cdot))\in L^{\infty}(0,\infty)$ and
$k_{\alpha}*\|\partial_{t}^{\alpha}(u-u_{0})\|_{H}^{2}\in
L^{\infty}(0,\infty)$, and moreover,

$\displaystyle\frac{1}{2}\left(k_{\alpha}*\|\partial_{t}^{\alpha}(u-u_{0})\|_{H}^{2}\right)(t)+\phi(u(t))\leq\phi(u_{0})$

for a.e. $t>0$.

  4. (iv)

Let $u_{j}$ be a global strong solution to (32) with an initial datum
$u_{0,j}$ for $j=1,2$. Then it holds that

$\displaystyle\|u_{1}(t)-u_{2}(t)\|_{H}\leq\|u_{0,1}-u_{0,2}\|_{H}$ (34)

for a.e. $t>0$.

######  Remark.

  1. (i)

Even for $u_{0}\in\overline{D(\phi)}^{H}$, one can still construct a function
satisfying all the requirements except $k_{1-\alpha}*(u-u_{0})\in
W^{1,2}(0,T;H)$ for strong solutions to (32) in Definition 4.1. We can still
check $k_{1-\alpha}*(u-u_{0})\in W^{1,2}_{\rm loc}((0,T];H)$, which may not,
however, imply the uniqueness of such “solutions” immediately.

  2. (ii)

From the theorem above, we know that $\phi(u(t))\leq\phi(u_{0})$ for a.e.
$t>0$; however, it may still be open to question whether $t\mapsto\phi(u(t))$
is nonincreasing or not. On the other hand, the nonlocal quantity
$\mathcal{F}(t)$ is surely continuous and nonincreasing in $t$.

  3. (iii)

The assertions except (iv) of Theorem 4.2 can also be extended to _semiconvex_
functionals $\phi:H\to(-\infty,+\infty]$, that is, the case where

$\phi(w)+\frac{\lambda}{2}\|w\|_{H}^{2}$

is convex for some $\lambda\in\mathbb{R}$. In particular, $\phi$ may not be
convex when $\lambda>0$. Then the contraction property like (34) no longer
holds, but a Lipschitz property still holds true, i.e., for each $T>0$ there
exists a constant $C_{T}>0$ depending only on $T$, $\alpha$ such that

$\|u_{1}(t)-u_{2}(t)\|_{H}\leq C_{T}\|u_{0,1}-u_{0,2}\|_{H}$

for a.e. $t\in(0,T)$. See [121, Section 5] for more details.

We briefly review two key ingredients of the proof for Theorem 4.2 in [121],
that is, an operator theoretic frame for fractional derivatives and a
fractional chain-rule inequality. Let $\mathcal{H}=L^{2}(0,T;H)$ and define an
operator $\mathcal{B}:D(\mathcal{B})\subset\mathcal{H}\to\mathcal{H}$ by

$\displaystyle D(\mathcal{B})$ $\displaystyle=\\{w\in\mathcal{H}\colon
k_{1-\alpha}*w\in W^{1,2}(0,T;H),\ (k_{1-\alpha}*w)(0)=0\\},$
$\displaystyle\mathcal{B}(u)$
$\displaystyle=\partial_{t}\left(k_{1-\alpha}*w\right)\ \mbox{ for }\ w\in
D(\mathcal{B}).$

Then $\mathcal{B}$ is maximal monotone in $\mathcal{H}\times\mathcal{H}$, and
moreover, the Yosida approximation $\mathcal{B}_{\lambda}$ of $\mathcal{B}$ is
represented as

$\mathcal{B}_{\lambda}(w)=\frac{\text{\rm d}}{\text{\rm
d}t}\left(k_{1-\alpha,\lambda}*w\right)\quad\mbox{ for }\ w\in\mathcal{H},$

where $k_{1-\alpha,\lambda}$ is a nonnegative nonincreasing function of class
$W^{1,1}_{\rm loc}([0,\infty))$ such that

$k_{1-\alpha,\lambda}\to k_{1-\alpha}\quad\mbox{ strongly in }L^{1}(0,T)\
\mbox{ as }\ \lambda\to 0_{+}\ \mbox{ for any }\ T>0.$

We now move on to a fractional chain-rule inequality. Let us start with

######  Proposition 4.3 ([121, Proposition 3.4]).

Let $T>0$ and let $h$ be a nonnegative function of class $W^{1,1}(0,T)$. Let
$\psi:H\to(-\infty,+\infty]$ be a proper lower-semicontinuous convex
functional. Let $u\in L^{1}(0,T;H)$ be such that $\psi(u(\cdot))$ is
integrable over $(0,T)$. Then for any $t\in(0,T)$ satisfying $u(t)\in
D(\partial\psi)$ it holds that

$\displaystyle{\left(\dfrac{\text{\rm d}}{\text{\rm
d}t}[h*(u-u_{0})](t),g\right)_{H}}$ $\displaystyle=\dfrac{\text{\rm
d}}{\text{\rm
d}t}[h*\psi(u(\cdot))](t)+h(t)\left[(u(t)-u_{0},g)_{H}-\psi(u(t))\right]$
$\displaystyle\quad+\int^{t}_{0}h^{\prime}(\tau)\left[(u(t-\tau)-u(t),g)_{H}+\psi(u(t))-\psi(u(t-\tau))\right]\,\text{\rm
d}\tau$

for any $u_{0}\in H$ and $g\in\partial\psi(u(t))$. In addition, if $u_{0}\in
D(\psi)$, $h\geq 0$ and $h^{\prime}\leq 0$, then it holds that

$\left(\dfrac{\text{\rm d}}{\text{\rm
d}t}[h*(u-u_{0})](t),g\right)_{H}\geq\dfrac{\text{\rm d}}{\text{\rm
d}t}\left[h*\big{(}\psi(u(\cdot))-\psi(u_{0})\big{)}\right](t).$

The particular choice $\psi(w)=(1/2)\|w\|_{H}^{2}$ yields (30) immediately. We
usually apply the proposition mentioned above for the regularized kernel
$h=k_{1-\alpha,\lambda}\in W^{1,1}(0,T)$ for $T>0$. Noting that
$\mathcal{B}_{\lambda}(w)\to\mathcal{B}(w)$ in $\mathcal{H}$ as $\lambda\to
0_{+}$ for $w\in D(\mathcal{B})$, which is a well-known property of the Yosida
approximation, we find that

$\displaystyle{\left(\dfrac{\text{\rm d}}{\text{\rm
d}t}[k_{1-\alpha}*(u-u_{0})](t),g\right)_{H}}$
$\displaystyle=\left(\mathcal{B}_{\lambda}(u-u_{0})(t),g\right)_{H}+\left(\mathcal{B}(u-u_{0})(t)-\mathcal{B}_{\lambda}(u-u_{0})(t),g\right)_{H}$
$\displaystyle\geq\frac{\text{\rm d}}{\text{\rm
d}t}\left[k_{1-\alpha,\lambda}*\big{(}\psi(u(\cdot))-\psi(u_{0})\big{)}\right](t)$
$\displaystyle\quad+\left(\mathcal{B}(u-u_{0})(t)-\mathcal{B}_{\lambda}(u-u_{0})(t),g\right)_{H},$

provided that $u-u_{0}\in D(\mathcal{B})$. Let $u\in D(\mathcal{B})$ and
$\xi\in\mathcal{H}$ be such that $u(t)\in D(\partial\psi)$ and
$\xi(t)\in\partial\psi(u(t))$ for a.e. $t\in(0,T)$. Substitute $g=\xi(t)$ in
the above. The integration of both sides over $(s,t)$ yields

$\displaystyle{\int^{t}_{s}\left(\dfrac{\text{\rm d}}{\text{\rm
d}t}[k_{1-\alpha}*(u-u_{0})](\tau),\xi(\tau)\right)_{H}\,\text{\rm d}\tau}$
$\displaystyle\geq\left[k_{1-\alpha,\lambda}*\big{(}\psi(u(\cdot))-\psi(u_{0})\big{)}\right](t)-\left[k_{1-\alpha,\lambda}*\big{(}\psi(u(\cdot))-\psi(u_{0})\big{)}\right](s)$
$\displaystyle\quad-\left\|\mathcal{B}(u-u_{0})-\mathcal{B}_{\lambda}(u-u_{0})\right\|_{\mathcal{H}}\|g\|_{\mathcal{H}}$

for $0\leq s\leq t\leq T$. Passing to the limit as $\lambda\to 0_{+}$, we
obtain the following fractional chain-rule inequality,

$\displaystyle{\int^{t}_{s}\left(\partial_{t}^{\alpha}(u-u_{0})(\tau),\xi(\tau)\right)_{H}\,\text{\rm
d}\tau}$
$\displaystyle\geq\left[k_{1-\alpha}*\big{(}\psi(u(\cdot))-\psi(u_{0})\big{)}\right](t)-\left[k_{1-\alpha}*\big{(}\psi(u(\cdot))-\psi(u_{0})\big{)}\right](s)$
(35)

for a.e. $0<s\leq t<T$, whenever $u-u_{0}\in D(\mathcal{B})$, $\xi\in
L^{2}(0,T;H)$ and $[u(t),\xi(t)]\in G(\partial\psi)$ for a.e. $t\in(0,T)$. In
particular, if $\psi(u)\in L^{\infty}(0,T)$, then (35) holds for all $0\leq
s\leq t\leq T$.

###  4.3 Applications to Nonlinear Time-Fractional PDEs

In this subsection, we discuss applications of Theorem 4.2 to a couple of
nonlinear PDEs. In what follows, we let $\Omega$ be an open subset of
$\mathbb{R}^{d}$ with (smooth) boundary $\partial\Omega$.

(1) Parabolic $p$-Laplace equation. We consider the following Cauchy-Dirichlet
problem:

$\displaystyle\partial_{t}^{\alpha}(u-u_{0})$ $\displaystyle=\Delta_{p}u\quad$
$\displaystyle\mbox{ in }\Omega\times(0,\infty),$ (36) $\displaystyle u$
$\displaystyle=0$ $\displaystyle\mbox{ on }\partial\Omega\times(0,\infty),$
(37)

where $1<p<\infty$, $0<\alpha<1$ and $\Delta_{p}$ denotes the so-called
$p$-Laplacian given by

$\Delta_{p}u=\mathrm{div}\,\left(|\nabla u|^{p-2}\nabla u\right).$

Here we remark that the initial condition is already included in the PDE (36)
due to the formulation of strong solutions for time-fractional gradient flows
(see also Definition 4.1) (this will also be applied to the following two
other problems). We set $H=L^{2}(\Omega)$ and $\phi:H\to[0,\infty]$ as
follows:

$\phi(w)=\begin{cases}\frac{1}{p}\int_{\Omega}|\nabla w|^{p}\,\text{\rm
d}x&\mbox{ if }\ w\in W^{1,p}_{0}(\Omega),\\\ \infty&\mbox{
otherwise.}\end{cases}$

Then $\phi$ is proper, lower-semicontinuous, and convex in $H$. Hence, thanks
to Theorem 4.2, the Cauchy-Dirichlet problem (36)–(37) admits a unique global
$L^{2}$-solution $u\in L^{\infty}(0,\infty;W^{1,p}_{0}(\Omega))\cap
H^{\alpha,2}_{0,{\rm loc}}([0,\infty);L^{2}(\Omega))$. We also refer the
reader to [122] for the $1$-Laplacian (i.e., $p=1$) case.

(2) Nonlinear diffusion equation. Let us consider the following Cauchy-
Dirichlet problem:

$\displaystyle\partial_{t}^{\alpha}(u-u_{0})$
$\displaystyle=\Delta\left(|u|^{m-1}u\right)\quad$ $\displaystyle\mbox{ in
}\Omega\times(0,\infty),$ (38) $\displaystyle u$ $\displaystyle=0$
$\displaystyle\mbox{ on }\partial\Omega\times(0,\infty),$ (39)

where $0<m<\infty$ and $0<\alpha<1$. In order to reduce the Cauchy-Dirichlet
problem to the abstract gradient flow (32), we set $H=H^{-1}(\Omega)$ and
$\phi:H\to[0,\infty]$ as follows:

$\phi(w)=\begin{cases}\frac{1}{m+1}\int_{\Omega}|w|^{m+1}\,\text{\rm
d}x&\mbox{ if }\ w\in L^{m+1}(\Omega),\\\ \infty&\mbox{
otherwise,}\end{cases}$

which turns out to be a proper lower-semicontinuous convex functional in $H$.
Then using Theorem 4.2, the Cauchy-Dirichlet problem (38)–(39) admits a unique
global $H^{-1}$-solution $u\in L^{\infty}(0,\infty;L^{m+1}(\Omega))\cap
H^{\alpha,2}_{0,{\rm loc}}([0,\infty);H^{-1}(\Omega))$ (cf. see also [123]).

(3) Allen-Cahn equation. The Allen-Cahn equation reads,

$\displaystyle\partial_{t}^{\alpha}(u-u_{0})$ $\displaystyle=\Delta
u-\beta(u)+\lambda u\quad$ $\displaystyle\mbox{ in }\Omega\times(0,\infty),$
(40) $\displaystyle u$ $\displaystyle=0\quad$ $\displaystyle\mbox{ on
}\partial\Omega\times(0,\infty),$ (41)

where $0<\alpha<1$, $\beta:\mathbb{R}\to\mathbb{R}$ is a maximal monotone
function and $\lambda\in\mathbb{R}$. We set $H=L^{2}(\Omega)$ again and

$\phi(w)=\begin{cases}\frac{1}{2}\int_{\Omega}|\nabla u|^{2}\,\text{\rm
d}x+\int_{\Omega}\hat{\beta}(u)\,\text{\rm
d}x-\frac{\lambda}{2}\int_{\Omega}|u|^{2}\,\text{\rm d}x&\mbox{ if }\ u\in
D(\phi),\\\ \infty&\mbox{ otherwise,}\end{cases}$

where $\hat{\beta}$ is a primitive function of $\beta$, that is,
$\beta=\partial\hat{\beta}$ and $D(\phi)=\\{u\in
H^{1}_{0}(\Omega)\colon\hat{\beta}(u)\in L^{1}(\Omega)\\}$. Then $\phi$ is
(not convex but still) semiconvex and lower-semicontinuous in $H$. Hence
thanks to Remark Remark, the Cauchy-Dirichlet problem (40)–(41) admits a
unique global $L^{2}$-solution $u\in
L^{\infty}(0,\infty;H^{1}_{0}(\Omega))\cap H^{\alpha,2}_{0,{\rm
loc}}([0,\infty);L^{2}(\Omega))$.

###  4.4 Towards Dynamics of Fractional Gradient Flows

Recall that due to the nonlocal nature of fractional derivatives in time, one
cannot expect a classical semigroup property of the solution operator
$S(t):u_{0}\mapsto u(t)$ associated with the fractional gradient flow (32). On
the other hand, we may think of a trajectory approach as follows: We denote

$I:=\left\\{[0,t]\colon t\geq 0\right\\}$

and

$\displaystyle\mathfrak{Sol}=D(\phi)\cup\\{u:[0,t]\to D(\phi)\colon u\mbox{ is
a strong solution of \eqref{cp} on }[0,t]$ $\displaystyle\mbox{for some
initial datum }u_{0}\in D(\phi)\mbox{ and }t>0\\}.$

Moreover, we define an operator $U:I\times D(\phi)\to\mathfrak{Sol}$ by

$U([0,t],u_{0}):=\begin{cases}u_{0}&\mbox{ if }\ t=0,\\\ u&\mbox{ if }\
t>0,\end{cases}$

where $u$ denotes the strong solution of (32) on $[0,t]$ with the initial
datum $u_{0}$, and furthermore, we generate a map $\overline{U([0,t],u_{0})}$
from $(-\infty,0]$ into $D(\phi)$ by

$\overline{U([0,t],u_{0})}(\cdot):=\begin{cases}U([0,t],u_{0})(\cdot+t)&\mbox{
if }\ \cdot\in[-t,0],\\\ u_{0}&\mbox{ if }\ \cdot\not\in[-t,0]\end{cases}$

for $t>0$ and $\overline{U(\\{0\\},u_{0})}(\cdot)\equiv u_{0}$ for $t=0$.
Define a phase space $\mathfrak{X}$ by

$\displaystyle\mathfrak{X}=\left\\{\overline{U([0,t],u_{0})}\colon u_{0}\in
D(\phi),\ t\geq 0\right\\}.$

For each $\tau\geq 0$, we define an evolution operator
$T_{\tau}:\mathfrak{X}\to\mathfrak{X}$ by

$T_{\tau}(\overline{U([0,t],u_{0})})=\overline{U([0,t+\tau],u_{0})}$

for each $\overline{U([0,t],u_{0})}\in\mathfrak{X}$. Then we immediately
observe that

$T_{0}(\overline{U([0,t],u_{0})})=\overline{U([0,t],u_{0})},$

and moreover, for $\tau_{1},\tau_{2}\geq 0$, we find that

$\displaystyle{T_{\tau_{2}}\circ T_{\tau_{1}}(\overline{U([0,t],u_{0})})}$
$\displaystyle=T_{\tau_{2}}(\overline{U([0,t+\tau_{1}],u_{0})})$
$\displaystyle=\overline{U([0,t+\tau_{1}+\tau_{2}],u_{0})}=T_{\tau_{1}+\tau_{2}}(\overline{U([0,t],u_{0})}).$

Hence the semigroup property $T_{\tau_{1}}\circ
T_{\tau_{2}}=T_{\tau_{1}+\tau_{2}}$ holds true. Furthermore, we can also
define a functional $\mathcal{E}:\mathfrak{X}\to\mathbb{R}$ by

$\displaystyle\mathcal{E}(\overline{U([0,t],u_{0})})$
$\displaystyle:=\left[k_{1-\alpha}*\big{(}\phi(U([0,t],u_{0})(\cdot))-\phi(u_{0})\big{)}\right](t)$
$\displaystyle=\left[k_{1-\alpha}*\big{(}\phi(u(\cdot))-\phi(u_{0})\big{)}\right](t),$

which is decreasing in $t$ due to Theorem 4.2. Therefore we obtain

$\mathcal{E}(\overline{U([0,t],u_{0})})\leq\mathcal{E}(\overline{U([0,s],u_{0})})$

if $s\leq t$. That is, $\mathcal{E}$ is a Lyapunov functional for the
dynamical system $(\mathfrak{X},(T_{\tau})_{\tau\geq 0})$.

##  5 Steady States & Linear Stability: Space Fractional PDEs

In Section 3, we have presented spectra and bifurcations for dissipative
subdiffusion equations. In this section, we focus on one-dimensional space-
fractional reaction–diffusion problems. The space-fractional case is
technically a bit easier to deal with as one still studies Markovian dynamics
and the spectral theory for the space-fractional Laplacian is a bit more
classical, or it is even encoded immediately in its definition for the
spectral fractional Laplacian (SFL) as introduced in Section 2.2. We proceed
to study the effect of the fractional derivative order on the spectra, steady
states and bifurcations, also considering three well-known PDEs in which the
standard Laplacian is replaced by a spectral fractional Laplacian. A generic
fractional reaction–diffusion system can be formulated as

$\partial_{t}u=-D(-\Delta)^{\gamma/2}{u}+F({u};p)=:D\Delta^{\gamma/2}{u}+F({u};p),$
(42)

where we recall the shorthand notation

$\Delta^{\gamma/2}:=-(-\Delta)^{\gamma/2}$

for the fractional diffusion operator with fractional order $\gamma\in(0,2)$,
considered on a bounded domain $\Omega\subset\mathbb{R}^{d}$ and with suitable
boundary conditions specified below for each example. As for classical
diffusion, a key role in this context is played by the eigenvalues of the
fractional diffusion operator. A natural goal, similar to Section 3, is to
study patterns and we are going to give a brief summary of the spectral
properties for different definitions and analytical stability results for
bifurcations and Turing instability.

###  5.1 Spectra for Space-fractional PDEs

As discussed in Section 2, the several definitions of the fractional Laplacian
in $\mathbb{R}^{d}$ are all equivalent in that setting. This is not the case
on a bounded domain, so different definitions also have different spectral
properties.

Regarding its spectral properties, the spectral fractional Laplacian of
fractional order $\gamma$ turns out to be the simplest case, as the
eigenvalues and the eigenfunctions are the $(\gamma/2)$-power of the
eigenvalues of the Laplacian and the very same eigenfunctions of the
Laplacian, respectively. In particular, we recall the spectral definition of
the fractional Laplacian $\Delta^{\gamma/2}$ [124] (cf. Section 2.2)

$\Delta^{\gamma/2}u(x)=-\sum_{j=1}^{\infty}(-\lambda_{j})^{\gamma/2}u_{j}\phi_{j}(x),\qquad
u_{j}:=\langle
u,\phi_{j}\rangle_{L^{2}_{\Omega}}=\int_{\Omega}u\phi_{j}~{}\textnormal{d}x$
(43)

where $\lambda_{j}$ and $\phi_{j}$ the eigenvalues and the normalized
eigenfunctions, respectively, of the standard Laplacian $\Delta$ on $\Omega$,
i.e., the solutions of the eigenvalue problem
$\Delta\phi_{j}=\lambda_{j}\phi_{j}$ with specific boundary conditions [125].
For the one-dimensional case $\Omega=(0,L)$ the Dirichlet spectral fractional
Laplacian has eigenpairs

$(\phi_{j},\tilde{\lambda}_{j})=(\sin\left({j\pi
x}/{L}\right),-\left({j\pi}/{L}\right)^{\gamma}),\quad j\geq 1,$ (44)

while the Neumann spectral fractional Laplacian has eigenpairs

$(\phi_{j},\tilde{\lambda}_{j})=(\cos\left({j\pi
x}/{L}\right),-\left({j\pi}/{L}\right)^{\gamma}),\quad j\geq 0.$ (45)

For other definitions, as for instance the integral fractional Laplacian,
obtaining results on the spectral properties is not as straightforward as for
the spectral definition. In [126], the differences between these two
operators, namely the spectral and the integral definitions, have been
emphasized by studying their spectral properties. In particular, with respect
to Dirichlet boundary conditions the first eigenvalue of the integral
fractional Laplacian is strictly less than that of the spectral one. We refer
also to [127, 128, 129] for more detailed studies.

###  5.2 Analytical Results for Turing Bifurcations

Pattern formation in reaction–diffusion systems with fractional diffusion has
been investigated in different contexts: predator–prey interactions [130], a
semiarid vegetation model [131], a model for coral reefs [132], and chemical
systems [133, 134], just to name a few. As for classical diffusion problems,
also in the fractional setting we expect to obtain conditions, depending on
the fractional exponent $\gamma$, that can lead to the appearance of non-
homogeneous solutions and how the fractional order influence the possibility
of pattern formation (and eventually the type) is an interesting topic.
However, the interplay between the fractional order and the other parameters
can be nontrivial. In particular it is not obvious if the presence of the
fractional operator leads to a common effect in this class of problems or it
depends on the particular structure of the equation (nonlinear part).
Considering the spectral definition, the modification of the classical Turing
instability analysis is quite straightforward, since the fractional order
enters as exponent of the eigenvalues in the diffusion matrix (hence in the
dispersion relation). With other definitions, the first step is to obtain the
corresponding linearized problem. To this end, the Fourier transformed
$k-$space is the right framework, since it clarifies the relation between
properties of the spectrum and the structure of the solutions; cf. Section 3.

We now consider three important benchmark problems: the Allen–Cahn equation
[24], the Swift–Hohenberg equation [27], and the Schnakenberg system [135] in
a non-local setting. In particular, we consider the spectral definition.

The fractional Allen–Cahn equation. The model with the standard Laplacian,
proposed in [24], is a well-studied equation for various polynomial
nonlinearities [136, 137, 138, 139]. In particular, the classical Allen–Cahn
PDE ($\gamma=2$) with a cubic nonlinearity is also known as the (real)
Ginzburg–Landau PDE [140], an amplitude equation or normal form for
bifurcations from homogeneous states [141, 142]. Moreover, it is also common
to consider an additional quintic term added to the Ginzburg–Landau PDE [143,
144].

Different versions of the fractional Allen–Cahn have been recently
investigated, see for instance [145, 146, 147]. Following [148], we take with
the fractional Allen–Cahn equation with cubic–quintic nonlinearities

$\partial_{t}u=\Delta^{\gamma/2}u+pu+u^{3}-qu^{5},\quad\quad u\in\mathbb{R},\
\ p,q\in\mathbb{R},$ (46)

on a bounded interval $\Omega\subset\mathbb{R}$ of length $L$ with homogeneous
Dirichlet boundary conditions and we consider $p$ as the main bifurcation
parameter. As for the standard Allen–Cahn equation ($\gamma=2$), the
homogeneous solution $u_{*}=0$ is a stationary steady state of equation (46).
Furthermore, using the Fourier property of the fractional Laplacian, the
bifurcation values from the homogeneous states can be analytically computed
[148] and are located at

$p_{j}=\left(\frac{j\pi}{L}\right)^{\gamma},\quad j\geq 1.$ (47)

These values depend on the fractional order $\gamma$, and in particular all
bifurcation points move towards $p=~{}1$ as $\gamma$ decreases. The
homogeneous steady state $u_{*}=0$ is stable for $p<p_{1}$ and unstable
otherwise. The destabilization of the homogeneous steady state then happens at
$p_{1}$. Thus, decreasing the fractional order $\gamma$ pushes the stability
region of the homogeneous state to greater values of the bifurcation parameter
$p$. As a last observation, from equation (46) we can see that the fractional
order has no influence on the order of the bifurcation points, as shown in
Figure 1(a).

The Fractional Swift–Hohenberg Equation. The classical Swift–Hohenberg
equation ($\gamma=2$) [149] is a widely studied model in pattern dynamics [27,
150, 151], and also a standard test case for deriving reduced
amplitude/modulation equations [12, 13, 30, 152, 153].

Recently, the Swift–Hohenberg equation with nonlocal reaction terms has been
investigated [154, 155, 148]. Here, as in [148], we consider the fractional
Swift–Hohenberg equation with competing cubic–quintic nonlinearities

$\partial_{t}u=-(1+\Delta^{\gamma/2})^{2}u+pu+qu^{3}-u^{5},\quad\quad
u\in\mathbb{R},\ \ p,q\in\mathbb{R},\ \ q>0,$ (48)

on a bounded interval $\Omega\subset\mathbb{R}$ of length $L$ with homogeneous
Dirichlet boundary conditions. The parameter $p$ is taken as the main
bifurcation parameter. As in the previous model, the homogeneous state $u_{*}$
is a stationary solution to the fractional equation (48). According to (44),
the bifurcations from the homogeneous states occur at

$p_{j}=\left(1-\left(\frac{j\pi}{L}\right)^{\gamma}\right)^{2},\qquad j\geq
1.$ (49)

For simplicity, we fix $L=m\pi$, $m\in\mathbb{N}_{>0}$. Thus, the homogeneous
state $u_{*}$ is stable for $p<0$ and first becomes unstable at $p_{c}=0$,
independently of the fractional order $\gamma$ and with wavenumber $m$. In
addition, note that as $\gamma$ tends to zero, the bifurcation points on the
homogeneous branch $p_{j}$ accumulate at $p_{c}=0$. Interestingly and in
contrast to the previous case, as the fractional order decreases, it may
happen that bifurcation points shift their position, as shown in Figure 1(b)
for different wavenumbers.

\begin{overpic}[width=433.62pt]{AC-nsc}
\put(-8.0,45.0){\rotatebox{90.0}{$p_{j}$}} \put(85.0,-2.0){$\gamma/2$}
\put(30.0,20.0){$p_{1}$} \end{overpic} (a) fractional Allen–Cahn

\begin{overpic}[width=433.62pt]{SH-nsc}
\put(-8.0,45.0){\rotatebox{90.0}{$p_{j}$}} \put(85.0,-2.0){$\gamma/2$}
\put(90.0,43.0){$p_{9}$} \put(90.0,65.0){$p_{11}$} \put(90.0,7.0){$p_{10}$}
\end{overpic} (b) fractional Swift–Hohenberg

Figure 1: Effect of the fractional order $\gamma$ on the bifurcation points
$p_{j}$ on the homogeneous branch for the fractional Allen–Cahn equation
(left) and the Swift–Hohenberg equation (right). Each line corresponds to a
different wavenumber, using the analytical formula (47) for Allen–Cahn with
$L=10$ and formula (49) for Swift–Hohenberg with $L=10\pi$.

The Fractional Schnakenberg System. In its classical formulation, the
Schnakenberg system [135] is one of the prototype reaction–diffusion systems
exhibiting Turing patterns [3]. The fractional Schnakenberg system has been
investigated in [148] with a focus on the bifurcation structure and on the
snaking branch between periodic branches appearing with additional
nonlinearities in the reaction [156, 23]. As in [148], we focus on the
fractional Schnakenberg system with a simple reaction part

$\partial_{t}\begin{pmatrix}u_{1}\\\ u_{2}\end{pmatrix}=\begin{pmatrix}1&0\\\
0&d\end{pmatrix}\begin{pmatrix}\Delta^{\gamma/2}u_{1}\\\
\Delta^{\gamma/2}u_{2}\end{pmatrix}+\begin{pmatrix}-u_{1}+u_{1}^{2}u_{2}\\\
p-u_{1}^{2}u_{2}\end{pmatrix},\quad\quad u_{1},u_{2}\in\mathbb{R},\ \
p,d\in\mathbb{R},\ \ d>1.$ (50)

on a bounded interval $\Omega\in\mathbb{R}$ of length $L$ together with
homogeneous Neumann boundary conditions. The usual bifurcation parameter is
$p$ but we also want to understand the effect of the fractional order
$\gamma$. As for the Schnakenberg system with standard Laplacian ($\gamma=2$),
the spatially homogeneous state $u_{*}=(p,{1}/p)$ is a stationary solution to
(50). Moreover, according to (45), we can compute analytically the bifurcation
values on the homogeneous branch, which are located at

$p_{j}=k_{j}^{\gamma/2}\sqrt{\frac{d(1-k_{j}^{\gamma})}{1+k_{j}^{\gamma}}},\quad
k_{j}=\frac{j\pi}{L},\quad j\geq 0,$ (51)

provided that $1-k_{j}^{\gamma}>0$ so that the associated wavenumber lies in
the Turing instability region. Moreover, the critical parameter $p_{c}$ at
which the first instability occurs is

$p_{c}=\sqrt{d(3-\sqrt{8})},$ (52)

thus not affected by the fractional diffusion. However, in contrast to the
previous cases, the critical value $k_{c}$ at the destabilisation of the
homogeneous state depends on the fractional order $\gamma$

$k_{c}(\gamma/2)=\left(\sqrt{\frac{d-\mu_{c}^{2}}{2d}}\right)^{\frac{2}{\gamma}}=\left(\sqrt{\sqrt{2}-1}\,\right)^{\frac{2}{\gamma}},$
(53)

which is plotted in Figure 1(a). Note that we scaled here by a factor of two,
as above, to have the horizontal axis of the bifurcation parameter normalized
to $[0,1]$. We observe that the wavelength of non-homogeneous stationary
solutions bifurcating from the first bifurcation point increases as $\gamma$
decreases. To keep the same wavenumber at the primary bifurcation point, the
domain can be “tuned” (enlarged) to the primary bifurcation for each
fractional order $\gamma$. In Figure 2(b), we show the effect of the
fractional order $\gamma$ on the bifurcation points $p_{j}$ on the homogeneous
branch. As for the Swift–Hohenberg equation, two bifurcation points
corresponding to $j=1$ and $j=9$ shift their position as $\gamma$ decreases.
Note that additional nonlinearties as in [156, 23] do not affect $p_{c}$ and
$k_{c}$.

\begin{overpic}[width=433.62pt]{s-kc} \put(-8.0,45.0){$k_{c}$}
\put(85.0,-2.0){$\gamma/2$} \end{overpic} (a) Critical value $k_{c}(\gamma/2)$

\begin{overpic}[width=433.62pt]{s-nsc}
\put(-8.0,45.0){\rotatebox{90.0}{$p_{j}$}} \put(85.0,-2.0){$\gamma/2$}
\put(92.0,18.0){$p_{1}$} \put(92.0,55.0){$p_{9}$} \put(92.0,7.0){$p_{0}$}
\end{overpic} (b) Bifurcation values $p_{j}(\gamma/2)$

Figure 2: Effect of the fractional order $s$ on critical value $k_{s}$ from
equation (53) (left) and on the bifurcation points $\mu_{j}$ on the
homogeneous branch from equation (51) (right) for the Schnakenberg system. In
the left panel, each line corresponds to a different wavenumber, using the
analytical formula (51) with $L=10\pi$ and $d=60$.

##  6 Traveling Waves: Space Fractional PDEs

For earlier reviews on front propagation in reaction-diffusion equations with
anomalous diffusion (including subdiffusive systems), see e.g. [157] and the
references therein.

######  Remark.

Subdiffusion limited reaction-diffusion equations like (18), as any typical
non-autonomous equation, do not support (non-trivial) traveling wave solutions
of the form (3) with nonzero speed. For example, consider the subdiffusion
limited reaction-diffusion equation in (18) with _(Djrbashian-)Caputo
fractional derivative_ $\tilde{\mathfrak{D}}_{0,t}^{\alpha}$ of order
$\alpha\in(0,1)$. Computing the (Djrbashian-)Caputo fractional derivative with
respect to $t\geq 0$ of the traveling wave ansatz $u(t,x)=u(x-st)=u(\zeta)$
with $\zeta=x-st$ and $s<0$ yields

$\displaystyle\tilde{\mathfrak{D}}_{0,t}^{\alpha}\big{[}u(x-st)\big{]}$
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-\sigma)^{-\alpha}\big{(}\partial_{\sigma}u(x-s\sigma)\big{)}\,\mathrm{d}\sigma$
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{0}^{t}(t-\sigma)^{-\alpha}\big{(}-su^{\prime}(x-s\sigma)\big{)}\,\mathrm{d}\sigma$
and using $\xi:=x-s\sigma$, $\zeta=x-st$, $\delta:=\zeta-\xi$,
$\displaystyle=\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{x}^{x-st}\frac{u^{\prime}(\xi)}{(x-st-\xi)^{\alpha}}\,\mathrm{d}\xi=\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{x}^{\zeta}\frac{u^{\prime}(\xi)}{(\zeta-\xi)^{\alpha}}\,\mathrm{d}\xi$
$\displaystyle=\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{0}^{-st}\frac{u^{\prime}(\zeta-\delta)}{\delta^{\alpha}}\,\mathrm{d}\delta.$

The integral is rewritten using the traveling wave variables $\zeta$ and
$\xi$, but still depends either on the spatial variable $x$ or on the time
variable $t$. Therefore, a subdiffusion limited reaction-diffusion equation
(18) with Caputo fractional derivative can only support traveling wave
solutions of the form $u(t,x)=u(x-st)$ that are either stationary ($s=0$) or
steady states (if $s\neq 0$). However, in the natural long-term memory limit
(for fixed $\zeta$),

$\displaystyle\lim_{t\to\infty}\tilde{\mathfrak{D}}_{0,t}^{\alpha}\big{[}u(x-st)\big{]}$
$\displaystyle=\lim_{t\to\infty}\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{0}^{-st}\frac{u^{\prime}(\zeta-\delta)}{\delta^{\alpha}}\,\mathrm{d}\delta=\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{0}^{\infty}\frac{u^{\prime}(\zeta-\delta)}{\delta^{\alpha}}\,\mathrm{d}\delta$
$\displaystyle=\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{-\infty}^{\zeta}\frac{u^{\prime}(\delta)}{(\zeta-\delta)^{\alpha}}\,\mathrm{d}\delta,$

the traveling wave equation of (18) turns into an integro-differential
equation of the form,

$\frac{(-s)^{\alpha}}{\Gamma(1-\alpha)}\int_{-\infty}^{\zeta}\frac{u^{\prime}(\delta)}{(\zeta-\delta)^{\alpha}}\,\mathrm{d}\delta=\sigma\partial_{\zeta}^{2}u+f(u),\qquad\zeta\in\mathbb{R},$
(54)

cf. [158, 157]. For the subdiffusive Allen-Cahn equation, a formal analysis of
bistable fronts in [158] predicts algebraic decay in one of the tails due to
memory effects. For the time-fractional Fisher–KPP equation (18) with Caputo
derivative $\tilde{\mathfrak{D}}_{0,t}^{\alpha}$, the existence of invading
(i.e., monostable) fronts satisfying (54) has been proven in [159] (for
complementary results see [160, 161]). Some numerical simulations (see e.g.,
[159, §2]) suggest that solutions of the initial value problem for (18) with
initial datum given by a Heaviside step function converge to traveling wave
solutions of (54). However, we are not aware of an analytical proof of the
convergence of solutions of (18) to traveling wave fronts of (54) for
$\alpha\in(0,1)$.

In the remaining part of this section, we are going to focus on the existence
of traveling wave solutions for reaction-diffusion equations with anomalous
diffusion of the form

$\partial_{t}{u}=\mathcal{L}_{Levy}u+f(u)\,,\quad x\in\mathbb{R}^{d}\,,\quad
t\in(0,\infty)\,,$ (55)

where $f\in C^{\infty}(\mathbb{R})$ is a nonlinear function and
$\mathcal{L}_{Levy}$ is a Lévy operator (i.e., the generator of a Lévy
process).

###  6.1 Traveling Wave Solutions in One-Dimensional Settings

A particularly important subclass arises in the Riesz-Feller case leading to
evolution equations

$\partial_{t}{u}=\operatorname{D}_{\theta}^{\gamma}u+f(u)\,,\quad
x\in\mathbb{R}\,,\quad t\in(0,\infty)\,,$ (56)

where $f\in C^{\infty}(\mathbb{R})$ is a nonlinear function of bistable type,
i.e., $f$ has precisely three roots $u_{1}<u_{2}<u_{3}$ in the interval
$[u_{1},u_{3}]$ such that

$\displaystyle f(u_{1})=f(u_{2})=f(u_{3})=0\,,\quad
f^{\prime}(u_{1})<0\,,\quad f^{\prime}(u_{3})<0,$ (57a) $\displaystyle
f<0\text{ on }(u_{1},u_{2})\text{ and }f>0\text{ on }(u_{2},u_{3}),$ (57b)

and we recall that $\operatorname{D}_{\theta}^{\gamma}$ denotes a Riesz–Feller
operator for some fixed parameters $1<\gamma\leq 2$ and
$|\theta|\leq\min\\{\gamma,2-\gamma\\}$; see Section 2.2. To study the
traveling wave problem, it is necessary to extend the nonlocal operator to
$C^{2}_{b}(\mathbb{R})$. The following integral representations may be used to
accomplish this task:

######  Theorem 6.1.

If $0<\gamma<1$ or $1<\gamma<2$ and $|\theta|\leq\min\\{\gamma,2-\gamma\\}$,
then for all $f\in\mathcal{S}(\mathbb{R})$ and $x\in\mathbb{R}$

$\displaystyle\operatorname{D}_{\theta}^{\gamma}f(x)=\tfrac{c_{1}-c_{2}}{1-\gamma}f^{\prime}(x)$
$\displaystyle+c_{1}\int_{0}^{\infty}\,{\tfrac{f(x+y)-f(x)-f^{\prime}(x)\,y{\bf
1}_{(-1,1)}(y)}{y^{1+\gamma}}}\,\,\operatorname{d}\\!y\,$
$\displaystyle+c_{2}\int_{0}^{\infty}\,{\tfrac{f(x-y)-f(x)+f^{\prime}(x)\,y{\bf
1}_{(-1,1)}(y)}{y^{1+\gamma}}}\,\,\operatorname{d}\\!y\,$ (58)

where ${\bf 1}_{(-1,1)}(\cdot)$ is an indicator function and some constants
$c_{1},c_{2}\geq 0$ with $c_{1}+c_{2}>0$.

The result follows from [46, Theorem 31.7]; see also [162, Theorem 2.4]. A
traveling front solution of (56) is a traveling wave solution (3) connecting
distinct endstates $u_{\pm}\in\\{u_{1},u_{3}\\}$ in the sense that
$\lim_{\zeta\to\pm\infty}\phi(\zeta)=u_{\pm}$. The profile $u=\phi(\zeta)$ of
a traveling wave solution satisfies (the traveling wave equation)

$-s\phi^{\prime}(\zeta)=\operatorname{D}_{\theta}^{\gamma}{}[\phi](\zeta)+f(\phi),\qquad\zeta\in\mathbb{R},$
(59)

where $\operatorname{D}_{\theta}^{\gamma}$ has to be understood as its
extension (6.1) to $C^{2}_{b}$-functions. The existence of traveling front
solutions of (56) has been proved in case of
$\operatorname{D}_{0}^{2}=\partial_{x}^{2}$ [163, 164], fractional Laplacians
$\operatorname{D}_{0}^{\gamma}=-(-\partial^{2}_{x})^{\gamma/2}$,
$\gamma\in(0,2)$ modeling symmetric super-diffusion [165, 166, 167, 168], and
Riesz–Feller operators $\operatorname{D}_{\theta}^{\gamma}$ with $1<\gamma<2$
and $|\theta|\leq\min\\{\gamma,2-\gamma\\}$ in [162]. If a traveling front
solution of (1) with profile $\phi$ exists then its wave speed $s$ satisfies

$s=-\tfrac{\int_{u_{-}}^{u_{+}}\,{f(w)}\,\,\operatorname{d}\\!w\,}{\int_{-\infty}^{\infty}\,{(\phi^{\prime}(x))^{2}}\,\,\operatorname{d}\\!x\,}\,.$
(60)

Thus the potential
$F(u):=F(u_{-})+\int_{u_{-}}^{u}\,{f(v)}\,\,\operatorname{d}\\!v\,$ indicates
which stable state, either $u_{-}$ or $u_{+}$, will replace the other one. In
case of a balanced potential, i.e., $F(u_{-})=F(u_{+})$ or
$\int_{u_{-}}^{u_{+}}\,{f(v)}\,\,\operatorname{d}\\!v\,=0$, a stationary
traveling front will exist, i.e., both stable states will co-exist. In
contrast, in case of an unbalanced potential,
$\int_{u_{-}}^{u_{+}}\,{f(v)}\,\,\operatorname{d}\\!v\,\neq 0$, the stable
state with smaller potential value will replace the one with larger potential
value, also called the metastable state. It is important to note that many of
the following results are restricted to balanced bistable functions $f$. For
bistable reaction-diffusion equations (56) with fractional Laplacian
$\operatorname{D}_{0}^{\gamma}$ where $0<\gamma<2$,

$\partial_{t}{u}=\operatorname{D}_{0}^{\gamma}u+f(u)\,,\quad
x\in\mathbb{R}\,,\quad t\in(0,\infty)\,,$ (61)

the first preliminary investigations of traveling wave solutions highlighted
some similarities and differences to the classical case (1), see [169, 170,
171]: Using the simple reaction function $f$ of McKean [172], traveling wave
solutions of (61) can be constructed explicitly [171]. This model may serve as
a benchmark example. For (61) with $\gamma\in(1,2)$ and bistable reaction
function $f(u)=u(1-u^{2})$, a variational formulation can be derived such that
(61) is the associated Euler-Lagrange equation [169]. For general reaction
functions, traveling front solutions (3) with profile $\phi$ and wave speed
$s$ satisfy

$-s\int_{-\infty}^{\infty}\,{(\phi^{\prime}(\zeta))^{2}}\,\,\operatorname{d}\\!\zeta\,=\int_{-\infty}^{\infty}\,{\operatorname{D}_{0}^{\gamma}{}[\phi](\zeta)\
u^{\prime}(\zeta)}\,\,\operatorname{d}\\!\zeta\,+\int_{u_{-}}^{u_{+}}\,{f(w)}\,\,\operatorname{d}\\!w\,.$
(62)

For monotone profiles $u$, the first term on the right-hand side vanishes such
that the wave speed $s$ satisfies again (60), see [170]. However, the tail
behavior of profiles is different: traveling wave profiles $u$ approach the
endstates at an algebraic rate $1/|\zeta|^{\gamma}$ for $1\leq\gamma<2$, in
contrast, to an exponential rate for $\gamma=2$, see [170]. The existence and
uniqueness of traveling front solutions of fractional Allen–Cahn equation (61)
was proven via different methods [165, 166, 168, 167] (which we discuss
below). The stability of traveling front solutions has been proven via the
squeezing method using sub- and supersolutions [162, 173], and via spectral
analysis [174]. Next, we discuss the different approaches to prove the
existence of traveling front solutions.

####  6.1.1 Approximation of the Diffusion Operator

Fractional Laplacians $\operatorname{D}_{0}^{\gamma}$ can be approximated by a
family of convolution operators $J_{\epsilon}\ast\phi-\phi$ with
$J_{\epsilon}\in L^{1}(\mathbb{R})$ such that
$J_{\epsilon}\ast\phi-\phi\longrightarrow\operatorname{D}_{0}^{\gamma}\phi$
for $\epsilon\to 0$ in a suitable sense. The associated (traveling wave)
equations $-s\phi^{\prime}=J_{\epsilon}\ast\phi-\phi+f(\phi)$ exhibit for all
sufficiently small $\epsilon>0$ a unique monotone solution
$(\phi_{\epsilon},s_{\epsilon})$ with $\phi_{\epsilon}^{\prime}>0$ see [175,
176]. The limit $\phi=\lim_{\epsilon\to 0}\phi_{\epsilon}$ exists and yields
the profile $\phi$ of a traveling front solution of (61).

####  6.1.2 Using an Extension Problem

The operator square root of the Laplacian $\operatorname{D}_{0}^{1/2}$ can be
obtained from the harmonic extension problem to the upper half space as the
operator that maps the Dirichlet boundary condition to the Neumann condition.
For fractional Laplacians and other integro-differential operators, a similar
characterization is available [177, 178] (see also Section 7.2.2 and the
related discussion of the SFL in (89)). From those characterizations,
properties of these integro-differential equations can be derived from purely
local arguments in the extension problems [177]. For example, the stationary
traveling wave problem

$0=\operatorname{D}_{0}^{\gamma}w+f(w)\quad\text{ in }\,\mathbb{R}$ (63)

for $\gamma\in(0,2)$ and a bistable function $f$ can be related to a nonlinear
boundary value problem for a partial differential equation. This allows one to
prove that stationary traveling front solutions $w(x)=\phi(x)$ of (63) exist
if and only if $f$ is bistable with a balanced potential; and are unique (up
to translations) under the additional assumption $f^{\prime}(u_{\pm})<0$, see
[165, Theorem 2.4]. Moreover, the asymptotic behavior of the profile of
traveling front solutions can be analyzed, see [165, Theorem 2.7].

####  6.1.3 Minimizers of a Nonlocal Energy Functional

Consider the nonlocal energy functional

$\mathcal{E}(w,\Omega):=\mathcal{K}(w,\Omega)+\int_{\Omega}\,{F(w(x))}\,\,\operatorname{d}\\!x\,$
(64)

where $\mathcal{K}(w,\Omega)$ is (related to) the squared
Gagliardo–Sobolev–Slobodeckij $H^{s}$ semi-norm of $w$, and $F$ is a double-
well potential with $F(u_{\pm})=0$. In the one-dimensional setting, stationary
traveling front solutions $w(x)=\phi(x)$ of (63) are local minimizers of the
functional $\mathcal{E}(w,\mathbb{R})$. For a bistable function $f\in
C^{1}(\mathbb{R})$ with balanced potential, there exists a unique (up to
translations) nontrivial global minimizer $w$ of (64). This minimizer $w$ is a
traveling front solution of (63), and is unique (up to translations) in the
class of monotone solutions of (63), see [168, Theorem 2].

####  6.1.4 Homotopy Arguments

In [165, 168], the existence of _stationary_ traveling front solutions of (61)
(i.e.) with _balanced_ potentials have been shown. For general bistable
reaction functions $f$ (with possibly unbalanced potential), the existence of
(unique) traveling front solutions of (61) can be proven via a homotopy from
the balanced case [167].

####  6.1.5 Semigroup Approach and Squeezing Method using Sub- and Super-
solutions

For reaction-diffusion equations supporting a comparison principle, the
squeezing method using sub- and supersolutions allows one to prove the
existence, uniqueness and local asymptotic stability with exponential decay of
traveling wave solutions. The squeezing method was formulated—first—for a
class of nonlinear (nonlocal) evolution equations [176] including (1) and

$\partial_{t}{u}=J\ast u-u+f(u)=\int_{\mathbb{R}}\,{J(x-y)\
u(y,t)}\,\,\operatorname{d}\\!y\,-u(t,x)+f(u(t,x))\,,\quad
x\in\mathbb{R}\,,\quad t>0\,,$ (65)

for non-negative functions $J\in L^{1}(\mathbb{R})$ and general bistable
functions $f\in C^{2}(\mathbb{R})$; and—later—extended to other classes of
reaction-diffusion equations, such as (56) see [162].

######  Remark (Metastability).

The Allen–Cahn equation (1) exhibits metastable behavior, i.e., some solutions
appear to be stationary on extremely long—yet still transient—time scales, see
e.g., [179] and references therein. The metastable behavior of the fractional
Allen–Cahn equation (61) has been studied via asymptotic analysis and
numerical simulations, see [145].

For an earlier review on numerical methods to construct traveling front
solutions for the fractional Allen–Cahn equation (61), see e.g. [180].

###  6.2 Traveling Wave Solutions in Multi-Dimensional Settings

Next, we discuss the existence of traveling wave solutions for reaction-
diffusion equations (55) with bistable function $f$ in the multi-dimensional
setting. Consider the classical bistable reaction-diffusion equation

$\partial_{t}{u}=\Delta u+f(u)\,,\quad x\in\mathbb{R}^{d}\,,\quad
t\in(0,\infty)\,,$ (66)

and bistable reaction function $f$ satisfying (57). Here, we briefly discuss
the theory of traveling front solutions to (66) and refer to [181, 182, 183]
for further information. Planar traveling waves in $\mathbb{R}^{d}$ are of the
form

$u(t,x):=\phi(x\cdot e-st)$ (67)

for some vector $e\in\mathbb{R}^{d}$ with $\|e\|=1$ and $s\in\mathbb{R}$. If a
planar traveling wave (67) is a solution of (66) then its profile $\phi$
satisfies a second-order ODE

$\phi^{\prime\prime}+s\phi^{\prime}(\zeta)+f(\phi(\zeta))=0\qquad\text{for
}\zeta\in\mathbb{R},$ (68)

just like in the one-dimensional setting. A planar traveling front is a planar
traveling wave (67) which satisfies
$\lim_{\zeta\to\pm\infty}\phi(\zeta)=u_{\pm}$ for some $u_{\pm}\in\mathbb{R}$.
The existence of planar traveling fronts of (66) and the properties of its
profile $\phi$ and speed $s$ follow again from the analysis of ODE (68), see
[163, 164]. For bistable reaction function $f$ satisfying (57), $u_{-}=u_{3}$
and $u_{+}=u_{1}$, there exists a monotone profile
$\phi\in\mathbb{C}^{1}(\mathbb{R})$ and a unique speed $s$ satisfying (60).
The level sets of a planar traveling front at level $c\in\mathbb{R}$ and time
$t\geq 0$ are

$\\{x\in\mathbb{R}^{d}:\ u(t,x)=c\\}=\\{x\in\mathbb{R}^{d}:\ x\cdot
e=\phi^{-1}(c)+st\\},$ (69)

i.e., the level sets are parallel hyperplanes which are orthogonal to the
direction of propagation $e$. Hence, they are called planar. Moreover, planar
traveling fronts are only unique up to shifts. In particular, these traveling
fronts are invariant in the moving coordinate frame with speed $s$ in the
direction $e$. In the one-dimensional setting, bistable traveling fronts are
asymptotically stable. In the multi-dimensional setting, the stability w.r.t.
perturbations of the initial data is more subtle, but see, e.g., [184, 185,
186, 187, 183, 188]. Robustness/stability of traveling front solutions w.r.t.
perturbations of the nonlinear bistable reaction function has been studied
also [189].

For classical bistable reaction-diffusion equations (66) in the multi-
dimensional setting, other kinds of traveling front solutions exist, see the
reviews in e.g. [182, 183]: For example, traveling fronts with non-planar
level sets exist such as (i) (axisymmetric) conical-shape, (ii) (non-
axisymmetric) pyramidal shape, and (iii) exponential/parabolic shape [190], to
name a few examples. Again, these traveling wave solutions (resp. their
profiles) are invariant in a moving coordinate frame with constant speed.

The notion of transition fronts has been introduced to encompass all traveling
front solutions which converge to $u_{\pm}$ far away from their level sets
uniformly in time [191, 182]. There exist transition fronts of the classical
bistable reaction-diffusion equation (66) which are not invariant in any
moving coordinate frame with constant speed [182, Theorem 2.9]. Thus,
traveling fronts are a (strict) subset of transition fronts. Nonetheless, for
(bistable) transition fronts a  _global mean speed_ can be defined which is
independent of the shape of level sets of the transition fronts, see [182].

######  Remark (Rigidity results for the Allen-Cahn equation).

In case of bistable reaction functions $f$ with balanced potential $F$,
Gibbons’ conjecture may be stated as follows [192, 193]: Let $u$ be a bounded
solution of the stationary Allen-Cahn equation

$0=\Delta u+u-u^{3},\quad x\in\mathbb{R}^{d},$ (70)

such that, say w.r.t. $x_{d}$,

$\lim_{x_{d}\to\pm\infty}u(x^{\prime},x_{d})=u_{\pm}\qquad\text{uniformly
w.r.t. }x^{\prime}:=(x_{1},\ldots,x_{d-1})\in\mathbb{R}^{d-1}$ (71)

then $u(x)=\phi(x_{d})$ for some $\phi:\mathbb{R}\to\mathbb{R}$. Gibbons’
conjecture for (70) is true for any $d\in\mathbb{N}$, see the references given
in [192, 193]. De Giorgi’s conjecture about the rigidity of the Allen-Cahn
equation (70) assumes the monotonicity of the bounded solution $u$ in one
direction, say $\partial u/\partial x_{d}>0$; instead of the uniform limits in
(71). Variants of De Giorgi’s conjecture have been proven for $d\leq 8$,
whereas for $d\geq 9$ counterexamples have been constructed, see the
references given in [190, 194, 195, 196]. The analogous conjecture for the
parabolic equation (66) does not hold, due to the existence of traveling wave
solutions with exponential/parabolic shaped level sets [190].

Consider the fractional bistable reaction-diffusion equation

$\partial_{t}{u}=-(-\Delta)^{\gamma/2}u+f(u)\,,\quad
x\in\mathbb{R}^{d}\,,\quad t\in(0,\infty)\,,$ (72)

with fractional Laplacians $(-\Delta_{x})^{\gamma/2}$, $\gamma\in(0,2)$ in
$\mathbb{R}^{d}$. If a planar traveling wave (67) with $\phi\in
C^{2}_{b}(\mathbb{R})$ is a solution of (72) then its profile $\phi$ satisfies
a one-dimensional traveling wave equation (59) with $\theta=0$, i.e.

$-s\phi^{\prime}(\zeta)=-(-\partial_{\zeta}^{2})^{\gamma/2}[\phi](\zeta)+f(\phi),\qquad\zeta\in\mathbb{R}.$
(73)

For (67) with $\phi\in C^{2}_{b}(\mathbb{R})$, the traveling wave equation
(73) is derived by using the identity

$\displaystyle-(-\Delta)^{\gamma/2}\big{[}\phi(x\cdot e-st)\big{]}$
$\displaystyle=c_{d,\gamma}\int_{\mathbb{R}^{d}}\,{\frac{\phi((x+y)\cdot
e-st)+\phi((x-y)\cdot e-st)-2\phi(x\cdot
e-st)}{2|y|^{d+\gamma}}}\,\,\operatorname{d}\\!y\,$
$\displaystyle=c_{1,\gamma}\Bigg{(}\int_{\mathbb{R}}\,{\frac{\phi(\zeta+y_{e})+\phi(\zeta-
y_{e})-2\phi(\zeta)}{2|y_{e}|^{1+\gamma}}}\,\,\operatorname{d}\\!y_{e}\,\Bigg{)}\Bigg{|}_{\zeta=x\cdot
e-st}$ $\displaystyle=-(-\partial_{\zeta}^{2})^{\gamma/2}[\phi](\zeta),$

where $y_{e}:=y\cdot e$, see also [197]. Thus, if there exist traveling front
solutions of the fractional Allen-Cahn equation (61) in the one-dimensional
setting, then these solutions of (61)—with profile $\phi$ and speed $s$—can be
used to construct families of planar traveling front solutions of (72) in the
multi-dimensional setting: For $e\in\mathbb{R}^{d}$ with $\|e\|=1$, the
function $u(x):=\phi(x\cdot e-st)$ is a planar traveling wave front solution
of (72).

For the fractional Allen–Cahn equation (72) with $\gamma\in(1,2)$, the
existence (and stability) of  _non-planar_ traveling wave solutions has been
proven in [198, 199, 200]. Stability of steady states for the fractional
bistable reaction-diffusion equation has been studied in [201].

######  Remark (Rigidity results for the fractional Allen-Cahn equation).

Following Gibbons and De Giorgi, analogous conjectures for

$0=-(-\Delta)^{\gamma/2}u+u-u^{3},\quad x\in\mathbb{R}^{d},$ (74)

have been studied: The fractional Gibbons’ conjecture holds for any
$d\in\mathbb{N}$, see [192] and references therein. The fractional De Giorgi’s
conjecture, i.e., De Giorgi’s conjecture for (74), has been partially proven,
see a review of results in [196, 193, 202, 203]. Notably, the parameter
regimes for $\gamma\in(0,1)$ and $\gamma\in[1,2)$ are very different, see the
discussion in [193, 202].

##  7 Discretization and Simulation: Fractional PDEs

Eventually, the analytical/pen-and-paper methods discussed in the previous
sections will reach their limits for fractional dissipative PDEs similarly to
the case of standard differential operators. Hence, we have to develop
numerical schemes to understand in more detail the dynamics we are interested
in. Two aspects make the numerical treatment of fractional (partial)
differential equations more challenging than that of standard (partial)
differential equations. First, fractional differential operators are
inherently nonlocal operators which brings about increased computational
complexity due to fully populated matrices or the necessity to track the
history in time-dependent problems. Second, solutions to fractional
differential equations typically exhibit singularities even for smooth input
data, which often lowers convergence rates, unless specialized techniques or
locally refined meshes are employed.

###  7.1 Numerical Methods for Fractional Operators in 1D

Several classes of numerical methods for fractional differential operators
have been proposed in the univariate setting. They can be broadly divided into
2 groups: the first group is related to discretizing the operator in the style
of a finite difference method, which leads to methods with a time-marching
flavor whereas the second class are Galerkin/collocation methods, which
typically involve solving a non-trivial linear system of equations. Detailed
discussions can be found in [204, 205] and in particular the excellent
monograph [206] to which we will frequently refer. The simplest model of a
time-fractional differential equation (FODE) takes the form

$D^{\alpha}_{0,t}u+\lambda u=f(t),\quad\mbox{ $t\in(0,T)$},$ (75)

where $\lambda\in{\mathbb{R}}$ and $D^{\alpha}_{0,t}$ is either the Riemann-
Liouville or the Caputo fractional differential operators introduced in
Section 2.1 and $\alpha\not\in{\mathbb{N}}$. More complex time-fractional PDEs
(FPDE) of the form

$D^{\alpha}_{0,t}u(t,x)+Au(x)=f(t,x)\quad\mbox{ $t\in(0,T),\quad
x\in\Omega\subset{\mathbb{R}}^{d}$}$ (76)

for spatial operators $A$ can successfully be treated numerically by combining
separate discretizations of $D^{\alpha}_{0,t}u$ and $Au$ in a way analogous to
the procedure in classical time-dependent problems, viz., the famous “method
of lines” in which first a spatial discretization is applied and then the
resulting (fractional) ODE system is discretized or the well-known “Rothe
method”/“method of transverse lines” in which first the temporal problem is
discretized and afterwards the resulting spatial problems are solved
numerically. In both approaches for (76), it is useful to have a good
understanding of discretizations for the scalar case (75). Indeed, a
diagonalization of the operator $A$ (or its discretization) leads to problems
of the form (75) so that (75) may be thought of as a simple model for (76). In
this section, we will not discuss the discretization of the spatial operator
$A$, and we will mostly restrict ourselves to the case $\alpha\in(0,1)$.

######  Remark.

A hallmark of solutions of (75) (and thus also (76)) is that they have a
(weak) singularity at the initial time $t=0$ [50], [207]; see also the remarks
in Section 4. Consider the case of Caputo derivatives and $\alpha\in(0,1)$.
For smooth $f$ and $\lambda=0$, solutions of (75) can be expected to have the
form $u(t)=t^{\alpha}\tilde{u}(t)+c$ for $c\in{\mathbb{R}}$ and smooth
$\tilde{u}$ since a particular solution of
$\mathfrak{\tilde{D}}^{\alpha}_{0,t}u=t^{\beta}$ is given by
$u(t)=\frac{\Gamma(\beta+1)}{\Gamma(\beta+\alpha+1)}t^{\alpha+\beta}$, [207,
Ex. 1]. For $\lambda\neq 0$, the solution structure is more involved in that
the behavior $O(t^{\alpha})$ is only the leading order behavior as
$t\rightarrow 0$. For example, the function $u(t)=E_{\alpha,1}(-\lambda
t^{\alpha})$ with the standard (entire) Mittag-Leffler function $E_{\alpha}$
solves the homogeneous equation $\mathfrak{\tilde{D}}^{\alpha}_{0,t}u+\lambda
u=0$, [50, Thm. 4.3]. Hence, for smooth $f$ one should expect solutions of
(75) to have singularity components $t^{j\alpha}$ for $j=1,2,\ldots$.

####  7.1.1 The L1 Scheme and Related Methods

The definitions of $\mathfrak{D}^{\alpha}_{0,t}$ and
$\mathfrak{\tilde{D}}^{\alpha}_{0,t}$ may motivate various approximations. For
example, for $\alpha\in(0,1)$ a discretization of
$\mathfrak{\tilde{D}}^{\alpha}_{0,t}u$ in the points
$0=t_{0}<t_{1}<\cdots<t_{n}$ may be taken to be

$\displaystyle\mathfrak{\tilde{D}}^{\alpha}_{0,t}u(t_{n})$
$\displaystyle=\frac{1}{\Gamma(1-\alpha)}\int_{\tau=0}^{t_{n}}(t_{n}-\tau)^{-\alpha}u^{\prime}(s)\,d\tau\approx\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}\int_{\tau=t_{k}}^{t_{k+1}}(t_{n}-\tau)^{-\alpha}\frac{u(t_{k+1})-u(t_{k})}{t_{k+1}-t_{k}}\,d\tau.$
(77)

This approximation is known as the L1 method. For the case of equidistant
points $t_{k}=k{\Delta t}$ the sum takes the form of a (discrete) convolution:

$\displaystyle\mathfrak{\tilde{D}}^{\alpha}_{0,t}u(t_{n})$
$\displaystyle\approx\frac{1}{\Gamma(1-\alpha)}\sum_{k=0}^{n-1}b_{n-k-1}\left[u(t_{k+1})-u(t_{k})\right],\qquad
b_{k}:=\frac{{\Delta
t}^{-\alpha}}{\Gamma(2-\alpha)}\left[(k+1)^{1-\alpha}-k^{1-\alpha}\right].$
(78)

A numerical method involving $\mathfrak{\tilde{D}}^{\alpha}_{0,t}$ is obtained
by replacing $u(t_{k})$ with approximations $u_{k}$. When solving numerically
a FODE of the form (75) with initial condition $u(0)=u_{0}$, the
discretization (78) of $\mathfrak{\tilde{D}}^{\alpha}_{0,t}$ leads to a
numerical method that has the form of a time-stepping scheme, i.e., the
approximations $u_{1}$, $u_{2},\ldots$ are determined in turn. In contrast to
classical ODEs, the cost of computing $u_{n}$ is not $O(n)$ but $O(n^{2})$ due
to the “memory” effect, i.e., the nonlocal nature of the fractional
derivative. It is possible to reduce the cost to from $O(n^{2})$ to $O(n\log
n)$ (see Section 7.1.3 below).

######  Remark.

The L1 scheme is obtained in (77) by replacing $u$ by its piecewise linear
interpolant $\Pi_{1}u$. For $\alpha\in(0,1)$, the L2-scheme is obtained by
approximating $u$ by the linear interpolant on the first interval $(0,t_{1})$
and on the interval $(t_{j-1},t_{j})$ by the quadratic polynomial $\Pi_{2}u$
interpolating in three consecutive nodes $t_{j-1}$, $t_{j}$, $t_{j+1}$, and by
the quadratic interpolant in the nodes $t_{n-2}$, $t_{n-1}$, $t_{n}$ for the
last interval $(t_{n-1},t_{n})$, [206, Chap. 4]. Likewise, an approximation of
the Caputo derivative for $\alpha>1$ can be based on locally higher order
polynomial interpolation, [205, Sec. 3]. A variant, which includes the
fractional form of the Crank-Nicolson scheme (if $\theta=1/2$), is known as
Alikhanov’s method (also known as $L2-L1_{\sigma}$) and approximates the
fractional derivative at intermediate points $t_{j+\theta}=\theta
t_{j}+(1-\theta)t_{j+1}$ using again polynomial approximations on the
intervals $(t_{j},t_{j+1})$ and the terminal interval $(t_{n},t_{n+\theta})$,
see [206, Sec. 4.1], [205, Sec. 3].

The L1 scheme has a truncation error $O({\Delta t}^{2-\alpha})$ for _smooth_
$u$, [206, Lem. 4.2]; however, since by Remark Remark the solution of FODEs
typically have at least a (weak) singularity near the initial point $t=0$,
lower convergence rates are observed, [208]. For example, for a semidiscrete
analysis of (76) with an elliptic operator $A$ on uniform temporal meshes
[206, Sec. 5.1] shows convergence $O({\Delta t})$ in the discrete
$\ell^{\infty}$ norm (see, e.g., [209, 208] for the fully discrete error
analysis). For smooth right-hand sides $f$ in (75) or (76), better convergence
rates can be obtained using one of the following two devices. The first one is
to use the “corrected L1 scheme”, which is based on a uniform mesh and
modifies the initial 2 steps of the time-marching formulation; for example,
errors $O(({\Delta t}/t_{n})^{2-\alpha})$ are then recovered for the case of
(76) and an elliptic operator $A$, [206, Thm. 5.3]. The second device is to
resort to graded meshes better to resolve the singularity at $t=0$. Currently
successful techniques to obtain error estimates for the L1 scheme and related
ones on uniform and non-uniform meshes include two strategies: The first one
uses new variants of Gronwall’s Lemma [206, Thm. 6.3] that are adapted to
FODEs and originate from earlier work on Volterra integral equations, [210].
The second technique exploits that the L1 scheme is associated with an
M-matrix and obtains estimates via suitably constructed barrier functions,
[206, Thm. 6.6].

####  7.1.2 Convolution Quadrature

The temporal fractional differential operators have a form that makes them
amenable to convolution quadrature (CQ), going back to [211, 212]; see also
[213] and [206, Chap. 3] for an overview. CQ realizes numerically operators of
the form (we follow the notation introduced by Lubich)

$(K(\partial_{t})u)(t):=\int_{s=0}^{t}k(t-s)u(\tau)\,d\tau$ (79)

in that good approximations at uniformly spaced points $t_{n}=nh$ are obtained
in the form

$K(\partial^{\Delta t}_{t})u(t_{n})=\sum_{j=0}^{n}w_{n-j}u(t_{j})$

for suitable weights $w_{j}$ (more below). Important properties of CQ are:  a)
determining the weights $w_{j}$ requires only the Laplace transform of $k$;
b) CQ include the standard (numerical) differentiation and integration
operators;  c) CQ has the structure of an algebra under composition, which
greatly simplifies the development of methods for operators that are the
composition of, e.g., an integral operator of the form (79) and a
differentiation operator as it appears in the definition of the operators
$\mathfrak{\tilde{D}}^{\alpha}_{0,t}$ and $\mathfrak{D}^{\alpha}_{0,t}$.

######  Example 7.1.

(Derivation of CQ) Several ways to motivate the form of CQ are given in [213,
Sec. 2]. One way, taken from [211], derives CQ from ODE solvers (here:
multistep methods), representing the kernel in terms of the inverse Laplace
transform, and using the so-called Z-transform. In more detail but proceeding
purely formally, we express $k$ in terms of the inverse Laplace transform and
write with the Laplace transform $K={\mathfrak{L}}k$ of the kernel $k$

$\displaystyle y(t)$
$\displaystyle:=(K(\partial_{t})u)(t)=\int_{\tau=0}^{t}k(t-\tau)u(\tau)\,d\tau=\frac{1}{2\pi{\rm
i}}\int_{s\in\sigma+{\rm
i}{\mathbb{R}}}K(s)\int_{\tau=0}^{t}e^{s(t-\tau)}u(s)\,d\tau\,ds$ (80)
$\displaystyle=\frac{1}{2\pi{\rm i}}\int_{s\in\sigma+{\rm
i}{\mathbb{R}}}K(s)v(t;s)\,ds,$

where the function $t\mapsto v(t;s)$ solves the initial value problem (IVP)
$v^{\prime}(t;s)=sv(t;s)+u(t)$, $v(0,s)=0$ (variation of constants formula).
Next, this IVP is solved with a numerical method to get an approximation
$v_{n}(s)\approx v(t_{n};s)$, which leads to the approximation
$y(t_{n})\approx y_{n}:=(2\pi{\rm i})^{-1}\int_{s\in\sigma+{\rm
i}{\mathbb{R}}}K(s)v_{n}(s)\,ds$. To compute $v_{n}(s)$, one may employ a
multistep method with $k$-steps (determined by the coefficients $\alpha_{i}$,
$\beta_{i}$), which takes the form

$\displaystyle\frac{1}{{\Delta t}}\sum_{j=n-k}^{n}\alpha_{n-j}v_{j}$
$\displaystyle=\sum_{j=n-k}^{n}\beta_{n-j}(sv_{j}+u(t_{j})),$

where the required additional initial data $v_{j}$, $j=-k+1,\ldots,-1$ are
taken to be zero. The values $v_{j}$ can be computed recursively. A closed