# Dissipative systems fractionally coupled to a bath

A. Vertessen1, R. C. Verstraten1, C. Morais Smith1 1Institute for Theoretical
Physics, Utrecht University, Princetonplein 5, 3584CC Utrecht, The Netherlands

(July 20, 2023)

###### Abstract

Quantum diffusion is a major topic in condensed-matter physics, and the
Caldeira-Leggett model has been one of the most successful approaches to study
this phenomenon. Here, we generalize this model by coupling the bath to the
system through a Weyl fractional derivative. The Weyl fractional Langevin
equation is then derived without imposing a non-Ohmic macroscopic spectral
function for the bath. By investigating the short- and long-time behavior of
the mean squared displacement (MSD), we show that this model is able to
describe a large variety of anomalous diffusion. Indeed, we find ballistic,
sub-ballistic, and super-ballistic behavior for short times, whereas for long
times we find saturation, and sub- and super-diffusion.

##  I Introduction

Quantum diffusion [1, 2, 3, 4, 5, 6] has been at the attention of researchers
for many years, posing the fundamental questions on the description of
dissipation in open quantum systems. A more complete understanding of these
mechanisms could lead to important breakthroughs in engineering devices that
have a better energy efficiency. Despite the existence of several models,
numerous unanswered questions remain. For instance, experimental evidence [7,
8, 9] indicates that diffusion does not always conform to the linear behavior
described by Einstein [10], but instead exhibits a power-law scaling in the
mean squared displacement (MSD), commonly known as anomalous diffusion [11,
12]. One such example is quantum diffusion on a fractal [13], where the
exponent scales with the fractal dimension. Typically, models that describe
anomalous diffusion assume a different environment compared to regular
diffusion. Our aim here is to achieve a description of anomalous diffusion
within the usual environment using instead a fractional interaction.

In order to build a dissipative model, the environment of a system should be
included in the Hamiltonian, since energy is conserved in quantum mechanics.
In the Caldeira-Leggett model [14, 15, 16, 17, 18, 19], the system of
interest, which we will take to be a single quantum particle, is coupled
linearly with a thermal bath that consists of harmonic oscillators. Although
this might seem arbitrary at first, it is a very good approximation for
coupling a system to an in-equilibrium environment, since each degree of
freedom of such an environment must oscillate around a local minimum. The
benefit of this approximation is that the harmonic oscillators can be
integrated out exactly, thus leading to an effective description of the
particle. In doing so, one needs to know the spectral function of the bath,
which is given by the imaginary part of the Fourier transform of the retarded
dynamical susceptibility. This spectral function is, however, often difficult
to measure experimentally, and is therefore often assumed to be linear (Ohmic)
for simplicity. This assumption on the reservoir leads to a Langevin equation,
which describes linear diffusion. However, more recent works [20, 21, 22, 23]
have shown that a power-law (non-Ohmic) bath leads to a fractional Langevin
equation, which will typically have colored noise due to the fluctuation-
dissipation theorem, but can also have white noise [21].

In the fractional Langevin equation, the first derivative friction term is
replaced by a fractional derivative [24, 25]. These are operators which
generalize the order of derivatives to be not just an integer number, like a
first or second derivative, but any real (or even complex) number. There are,
however, many different fractional derivatives, which are not all equivalent
[23, 26, 27], and this can often lead to different results. The first main
question then is: which fractional derivative should one use? For this, we
have to realize that fractional derivatives are non-local operators, which
means that it is important to specify the time domain. Refs. [20, 21] both use
a Caputo derivative, which is only possible if we have a clear initialization
time. Since we are interested in dynamics, the boundaries of the action are
taken from $-\infty$ to $+\infty$, which is an issue if we want to use Caputo
derivatives inside the action. Therefore, we must take the domain of the
fractional derivative to be the same as of the action. Furthermore, one would
like to preserve the Fourier transform, since this has proven to be a powerful
mathematical tool in physics. Hence, we would like a fractional derivative
that is compatible with it. Because of these two reasons, we thus select the
Weyl fractional derivative [28, 29, 30].

Here, we generalize the Caldeira-Leggett model by introducing a Weyl
fractional derivative coupling term and show that it can describe a whole
family of anomalous diffusion. Unlike previous works [23, 21, 20], we derive a
fractional Langevin equation without imposing a non-Ohmic macroscopic spectral
function. Motivated by the general velocity depended coupling in Ref. [31], we
consider a fractional Weyl derivative to couple the system to the bath. In
order to study the Lagrangian, we first derive the Weyl fractional Euler-
Lagrange equation in analogy with Ref. [32]. Then, we are able to derive the
Weyl fractional Langevin equation for an Ohmic bath, which we solve
analytically using Fox H-functions. We calculate the asymptotic behavior of
the MSD and show that it is comparable to previous results of a non-Ohmic
Caldeira-Leggett model [20, 33, 34, 35, 36, 37].

The outline of this paper is the following: In Sec. II we introduce the Weyl
fractional derivative because it is not well known. Then, in Sec. III we
derive the Weyl fractional Langevin equation. Its solution for the free case
is presented in Sec. IV. In Sec. V, we investigate the asymptotic behavior of
the MSD, and in Sec. VI we present our conclusions and outlook.

##  II The Weyl fractional derivative

Although the Weyl fractional derivative has been mentioned briefly in the
literature (see Refs. [38, 11, 21, 23, 39, 40]), its mathematical exploration
has been very limited, to the best of our knowledge. Therefore, we provide a
summary of its main properties. Let $\alpha\in\mathbb{R}$, then

$\displaystyle\prescript{W}{+}{\bm{D}_{t}^{\alpha}}f(t)=\mathcal{F}^{-1}\left\\{(i\omega)^{\alpha}\mathcal{F}\left[f(t);\omega\right];t\right\\};$
$\displaystyle\prescript{W}{-}{\bm{D}_{t}^{\alpha}}f(t)=\mathcal{F}^{-1}\left\\{(-i\omega)^{\alpha}\mathcal{F}\left[f(t);\omega\right];t\right\\}$
(1)

are called the Weyl and pseudo Weyl fractional derivative of order $\alpha$,
respectively111The used conventions for the Fourier transform and its inverse
are
$\displaystyle\mathcal{F}\left[f(t);\omega\right]=\frac{1}{2\pi}\int_{-\infty}^{\infty}dtf(t)e^{-i\omega
t},$
$\displaystyle\mathcal{F}^{-1}\left[f(\omega);t\right]=\int_{-\infty}^{\infty}d\omega
f(\omega)e^{i\omega t}.$ , where $(\pm i\omega)^{\alpha}$ is the principal
value of the complex power function [41]. The sign of the frequency determines
the rotational direction in complex space, which makes their combination of
particular interest, since only scaling remains: The combination

$\displaystyle\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}f(t)$
$\displaystyle=\mathcal{F}^{-1}\left\\{|\omega|^{2\alpha}\mathcal{F}\left[f(t);\omega\right];t\right\\}$
(2)

is known as the fractional Laplacian $\left(-\Box_{t}\right)^{2\alpha}$ [42,
43]. In particular, $\alpha=1/4$ yields the operator $\sqrt{-\Box_{t}}$ which
is encountered among others in the context of Pseudo QED [44, 45]. This theory
emerges from a projected QED [44] because in 2D systems like graphene, the
dynamics of the electrons is restricted to a 2D plane, whereas the photons
intermediate their interactions in 3D. Fractional Laplacians were shown to be
relevant to describe the quantum valley Hall effect in graphene [46], as well
as excitons in transition metal dichalcogenides [47, 48].

The Weyl and pseudo Weyl fractional derivatives are linear and real operators.
In addition, they commute
$\prescript{W}{-}{\bm{D}_{t}^{\vphantom{\beta}\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\vphantom{\alpha}\beta}}=\prescript{W}{+}{\bm{D}_{t}^{\vphantom{\alpha}\beta}}\prescript{W}{-}{\bm{D}_{t}^{\vphantom{\beta}\alpha}}$
and satisfy the semigroup property, meaning that their orders add up. Some
very important properties of the (pseudo) Weyl fractional derivative are that
their convolution satisfies

$\displaystyle
g(t)\ast\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}f(t)=\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}g(t)\ast
f(t)=\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}\left[g(t)\ast f(t)\right]$ (3)

and their partial integration formula reads

$\displaystyle\int_{-\infty}^{\infty}dt\;g(t)\prescript{W}{+}{\bm{D}^{\alpha}_{t}}f(t)$
$\displaystyle=\int_{-\infty}^{\infty}dt\;f(t)\prescript{W}{-}{\bm{D}^{\alpha}_{t}}g(t),$
(4)

where we call the attention to the change of sign, from
$\prescript{W}{+}{\bm{D}^{\alpha}_{t}}$ to
$\prescript{W}{-}{\bm{D}^{\alpha}_{t}}$. The proofs of these statements are
shown in the Appendix A. By making use of Eq. (4), one can derive the Euler-
Lagrange equations for Lagrangians of the form

$\displaystyle\mathcal{L}\left(t,y(t),\left\\{\prescript{W}{+}{\bm{D}_{t}^{\alpha_{k}}}y(t),\prescript{W}{-}{\bm{D}_{t}^{\beta_{j}}}y(t)\right\\}_{k,j}\right),$
(5)

where $k=1,...,n$, $j=1,...,m$, and all $\alpha_{k}$’s and $\beta_{j}$’s are
different from each other and non-zero to avoid over-counting. In this case,
the Weyl fractional Euler-Lagrange equation reads

$\displaystyle\frac{\partial\mathcal{L}}{\partial
y}+\sum_{k=1}^{n}\prescript{W}{-}{\bm{D}_{t}^{\alpha_{k}}}\frac{\partial\mathcal{L}}{\partial\prescript{W}{+}{\bm{D}_{t}^{\alpha_{k}}}y}+\sum_{j=1}^{m}\prescript{W}{+}{\bm{D}_{t}^{\beta_{j}}}\frac{\partial\mathcal{L}}{\partial\prescript{W}{-}{\bm{D}_{t}^{\beta_{j}}}y}=0.$
(6)

A detailed derivation of the above equation can be found in App. A. It is
partially inspired by Ref. [32], which however used a Riemann-Liouville
fractional derivative to a simpler version of Eq. (6). Note that some of the
derivatives could be of integer order. This is not excluded.

##  III Derivation of the Weyl fractional Langevin equation

The fractional Caldeira-Leggett model is an extension of the original model,
in which we modify the interaction term by introducing a Weyl fractional
derivative in the coupling between the system and bath. The rational behind
this modification is that the particle and the bath remain exactly the same,
however the interaction between them changes. Ref. [31] showed the effect of
modifications in this term. The alteration of the coupling could be attributed
to several external factors, such as the geometry in which the particle or
bath resides or the interaction mechanism. We consider the Lagrangian

$\displaystyle\mathcal{L}=\mathcal{L}_{p}+\mathcal{L}_{bath}+\mathcal{L}_{int}+\mathcal{L}_{CT},$
(7)

where

$\displaystyle\mathcal{L}_{p}$ $\displaystyle=\frac{1}{2}M\dot{Q}^{2}-V(Q),$
(8) $\displaystyle\mathcal{L}_{bath}$
$\displaystyle=\frac{1}{2}\sum_{k=1}^{N}m_{k}(\dot{q_{k}}^{2}-\omega_{k}^{2}q_{k}^{2}),$
(9) $\displaystyle\mathcal{L}_{int}$
$\displaystyle=\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q\sum_{k=1}^{N}C_{k}q_{k}.$
(10)

Here, $Q$ is the coordinate, $M$ is the mass of the system, and $V(Q)$ is the
potential that the system undergoes. The harmonic oscillators of the bath have
coordinates $q_{k}$, frequency $\omega_{k}$, and mass $m_{k}$.
$\alpha\in\mathbb{R}$ is the order of the fractional derivative and
$\mathcal{L}_{CT}$ is a counter-term Lagrangian. It is obtained by imposing
that the system of interest must be subject only to the potential $V(Q)$ and
not to some renormalized potential due to the coupling with the reservoir.
Therefore, we must impose that $\frac{\partial\mathcal{L}}{\partial
Q}\stackrel{{\scriptstyle!}}{{=}}-\frac{\partial V}{\partial Q}$. To do this,
we first set

$\displaystyle 0$ $\displaystyle=\frac{\partial\mathcal{L}}{\partial
q_{j}}=-m_{j}\omega_{j}^{2}q_{j}+C_{j}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q,$
(11)

because we want the bath to be non-interacting, meaning that each particle of
the bath should not feel any potential. This implies that

$\displaystyle q_{j}$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q.$
(12)

Now, we can find the minimum of $\mathcal{L}$ with respect to $Q$, which
should be the minimum of the potential, such that

$\displaystyle\frac{\partial\mathcal{L}}{\partial Q}$
$\displaystyle=-\frac{\partial V}{\partial Q}+\frac{\partial}{\partial
Q}\left(\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q\right)\sum_{k=1}^{N}C_{k}q_{k}+\frac{\partial\mathcal{L}_{CT}}{\partial
Q}$ $\displaystyle\stackrel{{\scriptstyle!}}{{=}}-\frac{\partial V}{\partial
Q}.$ (13)

Upon substituting Eq. (12) into Eq. (13) and solving for $\mathcal{L}_{CT}$,
we find

$\displaystyle\mathcal{L}_{CT}=-\frac{1}{2}\sum_{k=1}^{N}\frac{C_{k}^{2}}{m_{k}\omega_{k}^{2}}\left(\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q\right)^{2}.$
(14)

Now that we have the full Lagrangian, we can study the equations of motion.
The Euler-Lagrange equations for the bath are given by

$\displaystyle
m_{j}\ddot{q_{j}}=-m_{j}\omega_{j}^{2}q_{j}+C_{j}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q.$
(15)

Because we want to solve equations using Fourier transforms, we need to find
two solutions: a particular solution, which has no boundary conditions, and a
homogeneous solution, in which we set the dependence of external functions to
zero, yielding precisely the boundary conditions. Taking the Fourier transform
of Eq. (15) yields the particular solution, $q_{j}^{p}(t)$,

$\displaystyle q_{j}^{p}(\omega)$
$\displaystyle=-\frac{C_{j}}{m_{j}(\omega^{2}-\omega_{j}^{2})}(i\omega)^{\alpha}Q(\omega).$
(16)

In time, this can be written as

$\displaystyle q_{j}^{p}(t)=$
$\displaystyle\frac{C_{j}}{m_{j}\omega_{j}^{2}}\left\\{\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)\right.$
$\displaystyle-\left.\frac{d}{dt}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\left[\cos(\omega_{j}t)\Theta(t)\ast
Q(t)\right]\right\\}.$ (17)

The details of this derivation can be found in the Appendix B. The homogeneous
part of the solution, $q_{j}^{h}(t)$, is given by

$\displaystyle q_{j}^{h}(t)=A\cos(\omega_{j}t)+B\sin(\omega_{j}t),$ (18)

where $A$ and $B$ are constants determined by the initial conditions. Thus,
the full solution of Eq. (15) is given by

$\displaystyle q_{j}(t)$ $\displaystyle=q_{j}^{p}(t)+q_{j}^{h}(t).$ (19)

Now that we have found the solutions for the bath, we can use these in the
equation of motion of the system, which can be derived from the Weyl
fractional Euler-Lagrange equation. By applying Eq. (6) to Eq. (7), we obtain

$\displaystyle M\ddot{Q}(t)=$ $\displaystyle-\frac{\partial V(Q)}{\partial
Q}+\sum_{k=1}^{N}\Bigg{[}C_{k}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}q_{k}(t)$
$\displaystyle-\frac{C_{k}^{2}}{m_{k}\omega_{k}^{2}}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)\Bigg{]}.$
(20)

The bath coordinates can be eliminated by inserting Eqs. (17)-(19) into Eq.
(20), and we eventually obtain

$\displaystyle M\ddot{Q}(t)=-\frac{\partial V(Q)}{\partial Q}$
$\displaystyle-\sum_{k=1}^{N}\Bigg{\\{}\frac{C_{k}^{2}}{m_{k}\omega_{k}^{2}}\frac{d}{dt}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\left[\cos(\omega_{k}t)\Theta(t)\ast
Q(t)\right]$
$\displaystyle+C_{k}A\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\cos(\omega_{k}t)+C_{k}B\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\sin(\omega_{k}t)\Bigg{\\}}.$
(21)

Now, we continue by calculating the bath spectral function associated with the
Lagrangian (7), given by the imaginary part of the retarded dynamical
susceptibility of the bath

$\displaystyle J(\omega)$
$\displaystyle=\text{Im}\left[-i\theta(t-t^{\prime})\left\langle\left[F(t),F(t^{\prime})\right]\right\rangle\right]$
$\displaystyle=\text{Im}\left[\chi_{B}(\omega)\right],$ (22)

where $F(t)$ is the force acting on the bath by the system. The spectral
function can be interpreted as describing how much of each frequency is
available in the bath to interact with the system. The force coming from the
interaction term in Fourier space is given by
$F(\omega)=(i\omega)^{\alpha}\sum_{k=1}^{N}C_{k}q_{k}(\omega)$. Eq. (19) shows
that $q_{k}(\omega)$ is a function of $Q(\omega)$ only via the particular
solution, such that Eq. (16) can be directly used. The dynamical
susceptibility is thus given by

$\displaystyle\chi_{B}(\omega)$ $\displaystyle\equiv\frac{\partial
F(\omega)}{\partial
Q(\omega)}=-\sum_{k=1}^{N}\frac{C_{k}^{2}(i\omega)^{2\alpha}}{m_{k}(\omega^{2}-\omega_{k}^{2})}.$
(23)

Writing $(i\omega)^{2\alpha}=|\omega|^{2\alpha}e^{2\alpha
i\text{arg}(i\omega)}=|\omega|^{2\alpha}e^{\alpha
i\pi\text{sign}(\omega)}=|\omega|^{2\alpha}\left[\cos(\alpha\pi)+i\
\text{sign}(\omega)\sin(\alpha\pi)\right]$ yields

$\displaystyle\chi_{B}(\omega)$
$\displaystyle=-\sum_{k=1}^{N}|\omega|^{2\alpha}\frac{C_{k}^{2}}{m_{k}(\omega^{2}-\omega_{k}^{2})}$
$\displaystyle\times\left[\cos(\alpha\pi)+i\
\text{sign}(\omega)\sin(\alpha\pi)\right].$ (24)

Interpreting $1/(\omega^{2}-\omega_{k}^{2})$ as its principal value yields

$\displaystyle\frac{1}{\omega^{2}-\omega_{k}^{2}}$ $\displaystyle=\frac{\pi
i}{2\omega_{k}}\left[\delta(\omega+\omega_{k})-\delta(\omega-\omega_{k})\right].$
(25)

Substituting Eq. (25) into Eq. (24) and taking the imaginary part gives

$\displaystyle J(\omega)$
$\displaystyle=\text{Im}\left[\chi_{B}(\omega)\right]$
$\displaystyle=\frac{\pi}{2}|\omega|^{2\alpha}\cos(\pi\alpha)\sum_{k=1}^{N}\frac{C_{k}^{2}}{m_{k}\omega_{k}}\left[\delta(\omega-\omega_{k})-\delta(\omega+\omega_{k})\right]$
$\displaystyle=\frac{\pi}{2}\omega^{2\alpha}\cos(\pi\alpha)\sum_{k=1}^{N}\frac{C_{k}^{2}}{m_{k}\omega_{k}}\delta(\omega-\omega_{k})$
$\displaystyle=\omega^{2\alpha}\cos(\pi\alpha)J_{CL}(\omega),$ (26)

where in the second line we imposed that physical frequencies are positive,
$\omega>0$ and $\omega_{k}>0$. In the last line, we identified the usual
Caldeira-Leggett spectral function [14], which is given by

$\displaystyle
J_{CL}(\omega)=\frac{\pi}{2}\sum_{k=1}^{N}\frac{C_{k}^{2}}{m_{k}\omega_{k}}\delta(\omega-\omega_{k}).$
(27)

Eq. (26) suggests that the Lagrangian (7) can describe anomalous diffusion
without imposing a non-Ohmic macroscopic spectral function, because the
microscopic spectral function has a non-Ohmic prefactor $\omega^{2\alpha}$
that appears naturally. This is in agreement with investigations of spin-boson
models [49, 50, 51], where a spin-half particle is coupled to a bosonic
harmonic oscillator bath. Indeed, in Ref. [51], they show that the effect of a
sub-Ohmic bath can be fully characterized by a single effective interaction
term, which is quadratic in the interaction. This explains why we find the
modified spectral function with an extra $2\alpha$ power in the frequency. In
contrast to just assuming a non-linear spectral function, however, our method
provides a microscopic justification for the non-linearity from the coupling
term. In the limit of $\alpha\to 0$, we obtain the usual Caldeira Leggett
spectral function.

Returning back to Eq. (21), we identify the friction force $F_{fr}(t)$ and the
fluctuating force $f(t)$,

$\displaystyle
F_{fr}(t)=\sum_{k=1}^{N}\frac{C_{k}^{2}}{m_{k}\omega_{k}^{2}}\frac{d}{dt}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\left[\cos(\omega_{k}t)\Theta(t)\ast
Q(t)\right],$ (28) $\displaystyle
f(t)=\sum_{k=1}^{N}C_{k}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\left[A\cos(\omega_{k}t)+B\sin(\omega_{k}t)\right].$
(29)

Note that in the limit of $\alpha\to 0$, we get exactly the same friction and
fluctuating forces as in the usual Caldeira-Leggett model [14].

Starting with the friction force, we can insert Eq. (27) into Eq. (28) to
obtain

$\displaystyle F_{fr}(t)$
$\displaystyle=\frac{2}{\pi}\frac{d}{dt}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\int_{0}^{t}d\tau\int_{0}^{\infty}d\omega$
$\displaystyle\times\frac{J_{CL}(\omega)}{\omega}\cos[\omega(t-\tau)]Q(\tau).$
(30)

To verify whether the above expression actually is a friction term, we need to
acquire more knowledge of $J_{CL}(\omega)$. The spectral function tells us how
the bath reacts to a given frequency, which means that $J_{CL}(0)=0$ if the
bath is in equilibrium. Therefore, we can perform a Taylor approximation of
$J_{CL}(\omega)$ to get linear behavior for small frequencies. This is
effectively the same as taking the Caldeira-Leggett spectral function to be
Ohmic,

$\displaystyle\begin{split}J_{CL}(\omega)=\begin{cases}&\eta\omega\quad\text{
for }\omega<\Omega,\\\ &0\qquad\text{for
}\omega>\Omega,\end{cases}\end{split}$ (31)

where $\Omega$ is a high-frequency cut-off ($\Omega\to\infty$) and $\eta$ is a
macroscopic friction coefficient. Although the spectral function could in
principle have any shape, this is the most common assumption in literature
[14, 52]. With this assumption, we see that the frequency only remains in the
argument of the cosine in Eq. (30), which is a delta function representation.
Hence, we find that

$\displaystyle F_{fr}(t)$
$\displaystyle=2\eta\frac{d}{dt}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\left[\int_{0}^{t}d\tau\delta(t-\tau)Q(\tau)\right]$
$\displaystyle=\eta\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\dot{Q}(t),$
(32)

where, like in the calculation for the Caldeira-Leggett friction term, the
delta function is evaluated at the boundary, which gives an extra factor of
$1/2$. Notice that Eq. (32) reduces to a normal friction term in the limit of
$\alpha\to 0$, which is indeed when our model should reduce to the usual
Caldeira-Leggett model.

Next, we turn our focus to the fluctuating force given by Eq. (29). The first
step is to determine the constants $A$ and $B$ in Eq. (19). This is done by
only coupling the particle coordinate $Q(t)$ to the bath after $t=0$, such
that $A$ and $B$ can simply be found by inserting $t=0$ in Eq. (18). We thus
find

$\displaystyle A=q_{j}(0),\quad B=\dot{q_{j}}(0)/\omega_{j}.$ (33)

Inserting this into Eq. (28) yields

$\displaystyle
f(t)=\sum_{k=1}^{N}C_{k}\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\left[q_{k}(0)\cos(\omega_{k}t)+\frac{\dot{q_{k}}(0)}{\omega_{k}}\sin(\omega_{k}t)\right].$
(34)

Eq. (34) is the same as the fluctuating force in the usual Caldeira-Leggett
model up to $\prescript{W}{-}{\bm{D}_{t}^{\alpha}}$. Since this operator is
linear we ignore it during the calculation such that we can proceed
analogously to the calculation done for the Caldeira-Leggett model [21]. We
then find

$\displaystyle\langle f(t)\rangle=0,$ (35) $\displaystyle\langle
f(t)f(t^{\prime})\rangle=2\eta
k_{B}T\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{-}{\bm{D}_{t^{\prime}}^{\alpha}}\delta(t-t^{\prime}).$
(36)

We found that Eq. (36) actually describes colored noise. Indeed, using the
definition of the pseudo Weyl fractional derivative acting on the Fourier
representation of the delta function, we obtain

$\displaystyle\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{-}{\bm{D}_{t^{\prime}}^{\alpha}}\delta(t-t^{\prime})$
$\displaystyle=\frac{\sin(\pi\alpha)\Gamma(2\alpha+1)}{\pi}\left(t-t^{\prime}\right)^{-2\alpha-1}.$
(37)

This representation of colored noise is only valid for $(2\alpha+1)>0$.

Finally, substituting Eqs. (32) and (34) into Eq. (21) yields

$\displaystyle M\ddot{Q}(t)+\frac{\partial V(Q)}{\partial
Q}+\eta\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\dot{Q}(t)=f(t),$
(38)

with $\langle f(t)f(t^{\prime})\rangle=2\eta
k_{B}T\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{-}{\bm{D}_{t^{\prime}}^{\alpha}}\delta(t-t^{\prime})$
and $\langle f(t)\rangle=0$. Note that just by changing the coupling to the
bath, we obtained a modified version of the Langevin equation that we call the
Weyl fractional Langevin equation. There was no need to assume a non-Ohmic
bath explicitly. This is the first main result of of this paper.

##  IV Solving the Weyl fractional Langevin equation for the free case

When the particle is not subject to any external potential, Eq. (38) reduces
to

$\displaystyle
M\ddot{Q}(t)+\eta\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\dot{Q}(t)=f(t).$
(39)

Taking the Fourier transform of Eq. (39), rearranging and taking the inverse
Fourier transform, leads to the particular solution

$\displaystyle Q^{p}(t)$
$\displaystyle=\mathcal{F}^{-1}\left[\frac{f(\omega)}{-M\omega^{2}+\eta
i\omega|\omega|^{2\alpha}};t\right]$
$\displaystyle=\frac{1}{2\pi}f(t)\ast\mathcal{F}^{-1}\left(\frac{1}{-M\omega^{2}+\eta
i\omega|\omega|^{2\alpha}};t\right).$ (40)

The homogeneous solution of Eq. (39) is obtained by solving
$\omega\left(-M\omega+\eta i|\omega|^{2\alpha}\right)=0,$ thus yielding
$\omega=0\text{ or }\omega=i\left(\eta/M\right)^{1/(1-2\alpha)}.$ Hence, the
homogeneous solution reads

$\displaystyle Q^{h}(t)$
$\displaystyle=\mathcal{F}^{-1}\left\\{A\delta(\omega)+B\delta\left[\omega-i\left(\frac{\eta}{M}\right)^{\frac{1}{1-2\alpha}}\right]\right\\}$
$\displaystyle=C+D\text{exp}\left[-t\left(\frac{\eta}{M}\right)^{\frac{1}{1-2\alpha}}\right]$
$\displaystyle=C+D\text{exp}\left(-\frac{t}{t_{s}}\right),$ (41)

where $C$ and $D$ are constants that will be determined later and
$t_{s}=(M/\eta)^{1/(1-2\alpha)}$ is a characteristic time-scale of our system.
Combining Eqs. (40) and (41), we obtain

$\displaystyle Q(t)$ $\displaystyle=Q^{p}(t)+Q^{h}(t)$
$\displaystyle=\frac{1}{2\pi}f(t)\ast\mathcal{F}^{-1}\left(\frac{1}{-M\omega^{2}+\eta
i\omega|\omega|^{2\alpha}};t\right)$
$\displaystyle+C+D\text{exp}\left(-\frac{t}{t_{s}}\right).$ (42)

To determine the constants $C$ and $D$, we proceed in the same way as for the
fluctuating force. We consider that the particle is only coupled to the bath
for $t>0$, such that we can take $f(0)=0$. This implies that $C$ and $D$ are
determined solely by the homogeneous part in Eq. (42), which leads to

$\displaystyle\begin{split}\begin{cases}C=Q(0)+\dot{Q}(0)t_{s},\\\
D=-\dot{Q}(0)t_{s}.\end{cases}\end{split}$ (43)

Since our system is translation invariant, we can chose for simplicity
$Q(0)=0$. Eq. (42) then becomes

$\displaystyle Q(t)$
$\displaystyle=\frac{f(t)}{2\pi}\ast\mathcal{F}^{-1}\left(\frac{1}{-M\omega^{2}+\eta
i\omega|\omega|^{2\alpha}};t\right)$
$\displaystyle+\dot{Q}(0)t_{s}\left[1-\text{exp}\left(-\frac{t}{t_{s}}\right)\right].$
(44)

Next, we can use that $\langle f(t)\dot{Q}(0)\rangle=\langle f(t)\rangle=0$,
$\langle\dot{Q}(0)^{2}\rangle=k_{B}T/M$ and $\langle
f(t)f(t^{\prime})\rangle=2\eta
k_{B}T\prescript{W}{-}{\bm{D}_{t}^{\alpha}}\prescript{W}{-}{\bm{D}_{t^{\prime}}^{\alpha}}\delta(t-t^{\prime})$,
to calculate the MSD for the above solution. After partial integration, we
find

$\displaystyle\left\langle Q(t)^{2}\right\rangle$ $\displaystyle=\frac{\eta
k_{B}T}{2\pi^{2}}\int_{0}^{t}d\tau\left\\{\mathcal{F}^{-1}\left[\frac{\left(i\omega\right)^{\alpha}}{-M\omega^{2}+\eta
i\omega|\omega|^{2\alpha}};\tau\right]\right\\}^{2}$
$\displaystyle+\frac{k_{B}T}{M}t_{s}^{2}\left[1-\text{exp}\left(-\frac{t}{t_{s}}\right)\right]^{2}.$
(45)

This is the second main result of this paper. Next, we discuss its physical
implications.

Figure 1: The three different fundamental strips of a Mellin transform
involving the shown gamma functions in the legend. (a) the first fundamental
strip, (b) the second and (c) the third. Notice that green poles are on the
right of the fundamental strip in (a) and (b) but on the left and (c).

##  V Asymptotic expansions of the mean squared displacement

To determine the asymptotic behavior of Eq. (45), we need to perform the
asymptotic expansions of the inverse Fourier integral inside that equation. To
achieve this aim, we first simplify the integral by using the property of the
inverse Fourier of a product. Then, we take the Mellin transform of the
remaining integral. Next, we write the expression in a way that is compatible
with the definition of the Fox H-function. After that, we take the inverse
Mellin transform and identify the specific Fox H-function that solves the
problem. Finally, we use known asymptotic expansions of the Fox H-function to
obtain the asymptotic behavior of Eq. (45). This procedure was inspired by
Ref. [53]. Since the reader might not be familiar with the Mellin transforms
and the Fox H-function, we provide a short overview of their properties in
Apps. C and D, respectively. More specifically, we start by analysing

$\displaystyle
I(t)\equiv\mathcal{F}^{-1}\left[\frac{\left(i\omega\right)^{\alpha}}{-M\omega^{2}+\eta
i\omega|\omega|^{2\alpha}};t\right],$ (46)

which can be simplified to

$\displaystyle I(t)$
$\displaystyle=i\Theta(t)\ast\mathcal{F}^{-1}\left[\frac{\left(i\omega\right)^{\alpha}}{-M\omega+\eta
i|\omega|^{2\alpha}};t\right]$ $\displaystyle\equiv i\Theta(t)\ast
I^{\prime}(t).$ (47)

Now, we focus solely on $I^{\prime}(t)$. Because of convergence, we will take
$\alpha\neq 0,1/2$. The Mellin transform of $I^{\prime}(t)$ is given by

$\displaystyle F(s)$ $\displaystyle=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{s-\alpha}\times$
$\displaystyle\times\frac{\Gamma(s)\Gamma\left(1-\frac{\alpha-1}{2\alpha-1}-\frac{s}{2\alpha-1}\right)\Gamma\left(\frac{\alpha-1}{2\alpha-1}+\frac{s}{2\alpha-1}\right)}{\Gamma\left(\frac{\alpha^{2}}{2\alpha-1}+\frac{s}{2}\right)\Gamma\left(1-\frac{\alpha^{2}}{2\alpha-1}-\frac{s}{2}\right)}.$
(48)

Details of this derivation can be found in App. E. The pole structure of Eq.
(48) depends on $\alpha$, which has consequences when identifying a Fox
H-function. This results in a region of convergence of the Mellin transform,
known as the fundamental strip, that also depends on $\alpha$. The behavior of
the different gamma functions can be observed in Fig. 1. Therefore, we
separate the analysis in $3$ distinct cases: $\alpha>1$, $1/2<\alpha<1$ and
$\alpha<1/2$, which we will refer to as the first, second and third
fundamental strips, respectively. Here, we present the calculation for the
first fundamental strip. The calculations for the second and third fundamental
strips are similar and can be found in App. E.

In the first fundamental strip, for $\alpha>1$, we rewrite Eq. (48) in such a
way that after taking the inverse Mellin transform of that equation, we can
identify the Fox H-function given by

$I^{\prime}(t)=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{-\alpha}H^{21}_{23}\left[\frac{t}{t_{s}}\
\middle|\begin{array}[]{ccc}&\left(\frac{\alpha-1}{2\alpha-1},\frac{1}{2\alpha-1}\right),&\left(\frac{\alpha^{2}}{2\alpha-1},\frac{1}{2}\right)\\\
(0,1)&\left(\frac{\alpha-1}{2\alpha-1},\frac{1}{2\alpha-1}\right),&\left(\frac{\alpha^{2}}{2\alpha-1},\frac{1}{2}\right)\end{array}\right],$
(49)

where we again identified the characteristic timescale of the system, $t_{s}$.
Using the asymptotic expansion of the Fox H-function developed in App. D, we
find

$I^{\prime}(t)\sim\begin{cases}&1\qquad\>\text{ when }t\to 0,\\\
&t^{-\alpha}\quad\ \>\text{when }t\to\infty.\end{cases}$ (50)

Substituting Eq. (50) into Eq. (47) gives

$\displaystyle I(t)\sim\begin{cases}&t\qquad\>\ \ \ \text{ when }t\to 0,\\\
&t^{-\alpha+1}\quad\>\>\text{when }t\to\infty.\end{cases}$ (51)

For $t\to 0$, we Taylor expand the exponential inside Eq. (45) up to linear
order and discard it for $t\to\infty$. Substituting Eq. (51) into Eq. (45)
together with the approximations of the exponential function, as explained
before, results in

$\displaystyle\left\langle
Q(t)^{2}\right\rangle\sim\begin{cases}&t^{2}\quad\quad\ \ \>\>\text{ when
}t\to 0,\\\ &t^{-2\alpha+3}\quad\text{ when }t\to\infty,\end{cases}$ (52)

for $\alpha>1$. Combining this with the results from the second and third
fundamental strips derived in App. E, we find that the short-time behavior of
the MSD is given by

$\displaystyle\left\langle Q(t)^{2}\right\rangle\stackrel{{\scriptstyle t\to
0}}{{\sim}}\begin{cases}&t^{2}\quad\ \ \>\>\ \text{ for
}\alpha>\frac{-1-\sqrt{5}}{4},\alpha\notin\mathcal{P},\\\
&t^{2\alpha+3}\quad\>\text{for
}\alpha<\frac{-1-\sqrt{5}}{4},\alpha\notin\mathcal{P},\end{cases}$ (53)

with
$\mathcal{P}=\big{\\{}(-3-\sqrt{21})/4,(-1-\sqrt{3})/2,-1/2,0,(1+\sqrt{3})/2,(-3+\sqrt{21})/4,1/2,1\big{\\}}$.
On the other hand, the long-time limit yields

$\displaystyle\left\langle Q(t)^{2}\right\rangle\stackrel{{\scriptstyle
t\to\infty}}{{\sim}}\begin{cases}&t^{-2\alpha+3}\quad\>\text{ for
}\alpha>1,\alpha\notin\mathcal{P},\\\ &t^{2\alpha+1}\quad\ \>\ \text{ for
}\alpha<1,\alpha\notin\mathcal{P}.\end{cases}$ (54)

This is the third main result of this paper. Next, we discuss its implications
and compare it with results in the literature.

##  VI Discussion and Conclusion

We have introduced the Weyl fractional derivative which, due to its properties
with the boundary and Fourier transforms, can be directly used inside the
Lagrangian. We did this by imposing a fractional derivative in the coupling
term between the system and bath, which we called a fractional Caldeira-
Leggett model. Doing so required the introduction of the Weyl fractional
Euler-Lagrange equation to derive the equations of motion, which we solved for
the harmonic oscillators using Fourier transforms. We noted that the
fractional bath spectral function could be written in terms of a fractional
power prefactor and the original Caldeira-Leggett spectral function. Inserting
the commonly used Ohmic assumption for the Caldeira-Leggett spectral function,
we found an effective spectral function which has a power law of order
$2\alpha+1$. This led to a Weyl fractional Langevin equation with colored
noise. We provided the analytical solutions for several quantities in terms of
an inverse Fourier transform, which we could interpret as a Fox H-function
with the help of Mellin transforms. Finally, the asymptotic limits of these
solutions were calculated for a wide range of values. For short times, we
found both ballistic and non-ballistic regimes. This means that, in some
cases, friction can dominate the short-time dynamics. For long times we found
everything ranging from sub-, normal-, and super-diffusion, as well as
saturation, which could hint to a relation with glassy physics. Overall, we
found a very broad range of anomalous diffusion.

We are now in a position to compare our results, given by Eqs. (53) and (54),
with the results of Ref. [20], in which, to the best of our knowledge, a
fractional Langevin equation was investigated for the first time. First, we
note that we use the Weyl fractional derivative, while Ref. [20] used the
Caputo fractional derivative. It is known that the Caputo fractional
derivative and Weyl fractional derivatives lead to different results in a lot
of cases. One such example is when the fractional derivative of the
exponential function is calculated (see App. A). This makes it very
interesting to compare the results. The friction force in Eq. (39) is of order
$2\alpha+1$, which corresponds to the derivative order from Ref. [20] (which
we will call $\alpha_{L}$). Therefore, to compare the results, we must make
the substitution $\alpha=(\alpha_{L}-1)/2$ in Eq. (54). Second, we note that
Ref. [20] only considered the two cases, $0<\alpha_{L}<1$ and
$1<\alpha_{L}<2$. For those restrictions, our result from Eq. (54) becomes

$\displaystyle\left\langle Q(t)^{2}\right\rangle\stackrel{{\scriptstyle
t\to\infty}}{{\sim}}t^{\alpha_{L}},$ (55)

which is exactly what was found in Ref. [20]. It is remarkable that the same
result is obtained for a Caputo derivative instead of Weyl, and derived using
completely different techniques. Moreover, we generalized the result of the
Caputo derivative, which in Ref. [20] was obtained for $0<\alpha_{L}<2$
(corresponding to $-1/2<\alpha<1/2$), to the entire real line, except for
$\alpha=0$, $\alpha=\pm 1/2,\alpha=1$ (corresponding to $\alpha_{L}=0,1,2,3$,
which are well-known integer derivative cases) and
$\alpha=(-3\pm\sqrt{21})/4$, $\alpha=(-1\pm\sqrt{3})/2$, which are cases for
which a logarithm would appear because of integration. We can thus conclude
that our model, given by Eq. (7), describes sub-diffusion for $-1/2<\alpha<0$,
and super-diffusion for $0<\alpha<1$, while $1<\alpha<3/2$ yields sub-
diffusion. Remarkably, when $\alpha>3/2$ and when $\alpha<-1/2$ we find a
negative MSD exponent, which typically indicates saturation in the MSD [21],
possibly describing glassy states. The short-time behavior of the MSD exhibits
the usual ballistic ($\sim t^{2}$) behavior when
$\alpha>-(1+\sqrt{5})/4\approx-0.809$. However, when
$-3/2<\alpha<-(1+\sqrt{5})/4$ we find sub-ballistic behavior, meaning that the
MSD goes as $t^{p}$, where $0<p<2$, when $t\to 0$. In particular, we have that
$0<p<(5-\sqrt{5})/2\approx 1.382$. Finally, for $\alpha<-3/2$, we find super-
ballistic behavior of the MSD for $t\to 0$, which in our words means that the
short time behavior goes as a negative exponent. The different behavior of the
MSD for short and long times in terms of the parameter of our model, $\alpha$,
and the total order of the friction force, $\alpha_{L}$, is summarized in Fig.
2.

Figure 2: The behavior of the MSD for small and large times are shown in terms
of the used parameter in this paper, $\alpha$, but also in terms of the total
friction force order, $\alpha_{L}$, which are linked by
$\alpha_{L}=2\alpha+1$. Some specific values of $\alpha$, and in turn
$\alpha_{L}$, should be ignored in this plot, namely
$\alpha\in\big{\\{}(-3\pm\sqrt{21})/4,(-1\pm\sqrt{3})/2,\pm 1/2,0,1\big{\\}}$.
The light red crosses indicate the values excluded for simplicity while the
red crosses indicate the values for which our calculation did not apply.

When we calculated the spectral function, we were able to identify the regular
Caldeira-Leggett spectral function inside it, but there is some delicate
interpretation that needs to be done here. The spectral function describes the
susceptibility of the bath to an external system in the frequency domain.
However, as it turns out, this is not just a property of the bath, but also
depends on the way that you couple to it, which is almost always linear. This
is done because it is typically the lowest-order contribution, and it gives
excellent options to proceed mathematically, with e.g. completing the square.
Since baths are typically thought to be Ohmic (linear in frequency) in this
fashion, we believe that the most likely assumption should be to keep this
part identical, and thus retain an Ohmic $J_{CL}(\omega)\sim\omega$. In turn,
this means that the fractional coupling spectral function goes as
$J(\omega)\sim\omega^{2\alpha+1}$, which resulted in a Weyl fractional
Langevin equation. The results of the latter remarkably match with Ref. [20],
despite a different definition being used there, and even expand the domain of
available parameters. Previously, this equation has only been derived using
non-Ohmic baths, but this derivation shows that it can also be achieved for
Ohmic baths, as long as the coupling term is changed accordingly. Hence,
understanding the coupling mechanism is as important as understanding the bath
properties in order to describe (anomalous) quantum diffusion.

As we have demonstrated, the Weyl fractional Langevin equation can describe
anomalous diffusion for many different exponents. In experiments of quantum
transport on fractals [13], theory of diffusion on confining geometries [54],
and random walks [11, 55], several connections have been made to link either
the order of the derivative or the diffusion exponent itself to a fractal
dimension. Fractals can exhibit a dimension that is non-integer, such as the
Sierpiński gasket which has a dimension of $d_{f}\approx 1.58$. Ref. [13] has
shown that the MSD exhibits a regime where the exponent is equal to $d_{f}$.
Although the mechanism for this is still not fully understood, this could be
explained by an effect of the fractal dimension on the interaction between the
underlying material and the quantum particle. If we assume that
$\alpha=(d_{f}-1)/2$, then the Weyl fractional Langevin equation reproduces
this result.

Moving forward, an important aspect that requires further investigation is the
realization of a fractional derivative coupling as an effective theory. One
promising approach is to couple the bath to another bath and subsequently
integrate out the intermediate bath degrees of freedom to obtain an effective
bath with the desired fractional derivative coupling, or perhaps by some
interaction within the bath that could add frequency dependence. Another
potential method that could result in Eq. (7) is inspired by the projection
scheme used in Refs. [44, 45] and the fractal delta function used in Ref.
[56]. In this approach, one would project the initial theory down to a
fractal, which could result in a fractional coupling. Additionally, applying
fractal statistics to the bath could lead to the fractional derivative
coupling. Furthermore, it is also interesting to explore quantization for the
Weyl fractional derivative further than has already been done in the
literature [40, 57, 58]. To do this, one would first need to develop the
Hamiltonian formalism for these fractional derivatives, which can be inspired
by Ref. [59]. After that, one chooses either for canonical quantization [40],
path integral quantization [57] or the less known Weyl-Wigner quantization
[58]. The last is in our opinion the most natural one to use for the Weyl
fractional derivative and therefore the most promising method. Finally, we
believe that it would be interesting to see if previous works [49, 60, 61],
that assumes a non-Ohmic bath, could be reproduced by the approach used in
this paper, as it would introduce a new way to justify the models
microscopically.

###### Acknowledgements.

We thank Lars Fritz and Rodrigo Arouca for fruitful discussions. This work was
supported by the Netherlands Organization for Scientific Research (NWO, Grant
No. 680.92.18.05, C.M.S. and R.C.V.).

##  Appendix A Fractional derivative properties

As shortly mentioned in the introduction, there exist many fractional
derivatives. One that is often used in the literature [21, 20] is the Caputo
fractional derivative, defined as

$\displaystyle\prescript{C}{}{\bm{D}_{t}^{s}}f(t)=\frac{1}{\Gamma(n-s)}\int_{0}^{t}d\tau\left(t-\tau\right)^{n-s-1}f^{(n)}(\tau),$
(56)

where $n\in\mathbb{N}$ such that $n-1\leq s<n$. However, the Caputo derivative
is not compatible with Fourier transforms, and for this reason we use the Weyl
fractional derivative, which is defined in Eq. (1) in this paper. To
illustrate the difference between the Caputo and Weyl fractional derivative,
we shortly discuss what happens to the exponential function when applying both
of these fractional derivatives to it. One of the defining properties of the
exponential function is that it is an eigenfunction of the derivative, meaning
that the derivative of the exponential function is itself an exponential
function up to some scaling. When generalizing integer order derivatives, this
might be one of the properties that one wants to keep. The Weyl fractional
derivative is one of those generalizations that satisfies this property,
namely

$\displaystyle\prescript{W}{+}{\bm{D}_{t}^{\alpha}}e^{at}$
$\displaystyle=\mathcal{F}^{-1}\left\\{(i\omega)^{\alpha}\mathcal{F}\left[e^{at};\omega\right];t\right\\}$
$\displaystyle=\mathcal{F}^{-1}\left\\{(i\omega)^{\alpha}\delta(\omega+ia);t\right\\}$
$\displaystyle=a^{\alpha}e^{at}.$ (57)

This is in stark contrast with what happens when one applies the Caputo
fractional derivative to the exponential function. In that case, we obtain

$\displaystyle\prescript{C}{}{\bm{D}_{t}^{s}}e^{at}$
$\displaystyle=\sum_{i=0}^{\infty}\frac{\Gamma(n+1)}{\Gamma(n-s+1)}t^{n-\alpha}.$
(58)

Thus, the exponential function is clearly not an eigenfunction of the Caputo
fractional derivative. This is one of the many differences between the Caputo
and Weyl fractional derivatives.

In what follows, we present mathematical proofs of certain properties of the
Weyl fractional derivative. First, we note that the (pseudo) Weyl fractional
derivative is a linear operator, which directly follows from the definition
and linearity of integrals. Next, we show that these operators are indeed real
operators, meaning that

$\displaystyle\left[\prescript{W}{\pm}{\bm{D}_{t}^{\alpha}}f(t)\right]^{*}=\prescript{W}{\pm}{\bm{D}_{t}^{\alpha}}f(t).$
(59)

It indeed holds that

$\displaystyle\left[\prescript{W}{\pm}{\bm{D}_{t}^{\alpha}}f(t)\right]^{*}$
$\displaystyle=\left[\int_{-\infty}^{\infty}d\omega e^{i\omega t}(\pm
i\omega)^{\alpha}\int_{-\infty}^{\infty}\frac{d\tau}{2\pi}e^{-i\omega\tau}f(\tau)\right]^{*}$
$\displaystyle=\int_{-\infty}^{\infty}d\omega e^{-i\omega t}(\mp
i\omega)^{\alpha}\int_{-\infty}^{\infty}\frac{d\tau}{2\pi}e^{i\omega\tau}f(\tau)$
$\displaystyle=\int_{-\infty}^{\infty}d\omega e^{i\omega t}(\pm
i\omega)^{\alpha}\int_{-\infty}^{\infty}\frac{d\tau}{2\pi}e^{-i\omega\tau}f(\tau)$
$\displaystyle=\prescript{W}{\pm}{\bm{D}_{t}^{\alpha}}f(t).$ (60)

In addition,

$\displaystyle\prescript{W}{-}{\bm{D}_{t}^{\vphantom{\beta}\alpha}}\prescript{W}{+}{\bm{D}_{t}^{\vphantom{\alpha}\beta}}f(t)$
$\displaystyle=\mathcal{F}^{-1}\left\\{(-i\omega)^{\alpha}\mathcal{F}\left[\prescript{W}{+}{\bm{D}^{\alpha}_{t}}f(t);\omega\right];t\right\\}$
$\displaystyle=\mathcal{F}^{-1}\left\\{(-i\omega)^{\alpha}(i\omega)^{\beta}\mathcal{F}\left[f(t);\omega\right];t\right\\}$
$\displaystyle=\mathcal{F}^{-1}\left\\{(i\omega)^{\beta}\mathcal{F}\left[\prescript{W}{-}{\bm{D}^{\alpha}_{t}}f(t);\omega\right];t\right\\}$
$\displaystyle=\prescript{W}{+}{\bm{D}_{t}^{\vphantom{\alpha}\beta}}\prescript{W}{-}{\bm{D}_{t}^{\vphantom{\beta}\alpha}}f(t),$
(61)

which means that the pseudo and the Weyl fractional derivative commute. They
also satisfy the semi-group property, meaning that their orders add up. This
is proved by

$\displaystyle\prescript{W}{\pm}{\bm{D}_{t}^{\vphantom{\alpha}\beta}}\prescript{W}{\pm}{\bm{D}_{t}^{\vphantom{\beta}\alpha}}f(t)$
$\displaystyle=\mathcal{F}^{-1}\left\\{(\pm i\omega)^{\alpha}(\pm
i\omega)^{\beta}\mathcal{F}\left[f(t);\omega\right];t\right\\}$
$\displaystyle=\mathcal{F}^{-1}\left\\{(\pm
i\omega)^{\alpha+\beta}\mathcal{F}\left[f(t);\omega\right];t\right\\}$
$\displaystyle=\prescript{W}{\pm}{\bm{D}^{\alpha+\beta}_{t}}f(t).$ (62)

Next, we will proof Eq. (3). The Fourier transform of a convolution, where one
of the terms is a Weyl fractional derivative of a function $f$, can be
rewritten as

$\displaystyle\mathcal{F}\left[g(t)\ast\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}f(t);\omega\right]$
$\displaystyle=\mathcal{F}\left[g(t);\omega\right]\mathcal{F}\left[\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}f(t);\omega\right]$
$\displaystyle=\mathcal{F}\left[g(t);\omega\right](\pm
i\omega)^{\alpha}\mathcal{F}\left[f(t);\omega\right]$
$\displaystyle=\mathcal{F}\left[\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}g(t);\omega\right]\mathcal{F}\left[f(t);\omega\right]$
$\displaystyle=\mathcal{F}\left[\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}g(t)\ast
f(t);\omega\right].$ (63)

Furthermore, we know that

$\displaystyle\mathcal{F}\left\\{\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}\left[g(t)\ast
f(t)\right];\omega\right\\}$
$\displaystyle=\mathcal{F}\left[g(t);\omega\right](\pm
i\omega)^{\alpha}\mathcal{F}\left[f(t);\omega\right].$ (64)

Hence, we conclude, by uniqueness of the Fourier transform, that

$\displaystyle
g(t)\ast\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}f(t)=\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}g(t)\ast
f(t)=\prescript{W}{\pm}{\bm{D}^{\alpha}_{t}}\left[g(t)\ast f(t)\right].$ (65)

Next, we will proof the partial integration formula, namely Eq. (4). By
definition

$\displaystyle\int_{-\infty}^{\infty}dtg(t)\prescript{W}{+}{\bm{D}^{\alpha}_{t}}f(t)$
$\displaystyle=\int_{-\infty}^{\infty}dtg(t)\int_{-\infty}^{\infty}d\omega
e^{i\omega t}(i\omega)^{\alpha}\times$
$\displaystyle\times\frac{1}{2\pi}\int_{-\infty}^{\infty}d\tau
e^{-i\omega\tau}f(\tau)$ $\displaystyle=\int_{-\infty}^{\infty}d\tau
f(\tau)\int_{-\infty}^{\infty}d\omega
e^{i\omega\tau}(-i\omega)^{\alpha}\times$
$\displaystyle\times\frac{1}{2\pi}\int_{-\infty}^{\infty}dte^{-i\omega t}g(t)$
$\displaystyle=\int_{-\infty}^{\infty}dtf(t)\prescript{W}{-}{\bm{D}^{\alpha}_{t}}g(t),$
(66)

where, in the second line, we made the substitution $\omega\to-\omega$ and in
the last line we renamed $t\longleftrightarrow\tau$. Finally, we prove that
the equations of motion for Lagrangians of the form given in Eq. (5), are
given by Eq. (6). We start by taking a variation of the action, which results
in

$\displaystyle\delta S[y]=\int_{-\infty}^{\infty}dt\Bigg{(}$
$\displaystyle\frac{\partial\mathcal{L}}{\partial y}\delta
y+\sum_{k=1}^{n}\frac{\partial\mathcal{L}}{\partial\prescript{W}{+}{\bm{D}_{t}^{\alpha_{k}}}y}\delta\prescript{W}{+}{\bm{D}_{t}^{\alpha_{k}}}y$
$\displaystyle+\sum_{j=1}^{m}\frac{\partial\mathcal{L}}{\partial\prescript{W}{-}{\bm{D}_{t}^{\beta_{j}}}y}\delta\prescript{W}{-}{\bm{D}_{t}^{\beta_{j}}}y\Bigg{)}.$
(67)

Swapping the order of the variation and the fractional derivatives and
applying the partial integration formula yields

$\displaystyle\delta S[y]=\int_{-\infty}^{\infty}dt\Bigg{(}$
$\displaystyle\frac{\partial\mathcal{L}}{\partial
y}+\sum_{k=1}^{n}\prescript{W}{-}{\bm{D}_{t}^{\alpha_{k}}}\frac{\partial\mathcal{L}}{\partial\prescript{W}{+}{\bm{D}_{t}^{\alpha_{k}}}y}$
$\displaystyle+\sum_{j=1}^{m}\prescript{W}{+}{\bm{D}_{t}^{\beta_{j}}}\frac{\partial\mathcal{L}}{\partial\prescript{W}{-}{\bm{D}_{t}^{\beta_{j}}}y}\Bigg{)}\delta
y.$ (68)

Because Eq. (68) holds for all $\delta y$, setting that equation equal to zero
implies

$\displaystyle\frac{\partial\mathcal{L}}{\partial
y}+\sum_{k=1}^{n}\prescript{W}{-}{\bm{D}_{t}^{\alpha_{k}}}\frac{\partial\mathcal{L}}{\partial\prescript{W}{+}{\bm{D}_{t}^{\alpha_{k}}}y}+\sum_{j=1}^{m}\prescript{W}{+}{\bm{D}_{t}^{\beta_{j}}}\frac{\partial\mathcal{L}}{\partial\prescript{W}{-}{\bm{D}_{t}^{\beta_{j}}}y}=0.$
(69)

##  Appendix B Eliminating the bath

In this appendix, we show how to pass from Eq. (16) to Eq. (17), since it is
highly non-trivial. Taking the inverse Fourier transform of Eq. (16) yields

$\displaystyle q_{j}^{p}(t)$
$\displaystyle=-\mathcal{F}^{-1}\left[\frac{C_{j}}{m_{j}(\omega^{2}-\omega_{j}^{2})}(i\omega)^{\alpha}Q(\omega);t\right]$
$\displaystyle=\mathcal{F}^{-1}\left[\frac{(i\omega)^{\alpha}C_{j}}{m_{j}\omega_{j}^{2}}\left(1-\frac{\omega^{2}}{\omega^{2}-\omega_{j}^{2}}\right)Q(\omega);t\right]$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\Bigg{\\{}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)-\mathcal{F}^{-1}\left[(i\omega)^{\alpha}\frac{\omega^{2}}{\omega^{2}-\omega_{j}^{2}}Q(\omega);t\right]\Bigg{\\}}$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\Bigg{\\{}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)+\frac{d^{2}}{dt^{2}}\mathcal{F}^{-1}\left[(i\omega)^{\alpha}\frac{1}{\omega^{2}-\omega_{j}^{2}}Q(\omega);t\right]\Bigg{\\}}$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\Bigg{\\{}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)+\frac{1}{2\pi}\frac{d^{2}}{dt^{2}}\mathcal{F}^{-1}\left(\frac{1}{\omega^{2}-\omega_{j}^{2}};t\right)\ast\mathcal{F}^{-1}\left[(i\omega)^{\alpha}Q(\omega);t\right]\Bigg{\\}},$
(70)

where in the second line we separated a term that will cancel the term coming
from $\mathcal{L}_{CT}$ and in the last line we used the property of the
inverse Fourier of a product in our conventions. We further simplify the above
expression by making use of

$\displaystyle\mathcal{F}^{-1}\left(\frac{1}{\omega^{2}-\omega_{j}^{2}};t\right)$
$\displaystyle=-\frac{1}{2\omega_{j}}\Bigg{[}\mathcal{F}^{-1}\left(\frac{1}{\omega+\omega_{j}};t\right)-\mathcal{F}^{-1}\left(\frac{1}{\omega-\omega_{j}};t\right)\Bigg{]}=-\frac{2\pi}{\omega_{j}}\sin(\omega_{j}t)\Theta(t),$
(71)

where in the last line we interpreted the inverse Fourier integrals in the
positive $i\epsilon$-interpretation. By the positive
$i\epsilon$-interpretation, we mean that first we move the pole to the upper
half of the complex plane, see Fig. 3. After, we close the contour from above
(below) for positive (negative) times. This allows one to use the residue
theorem and Jordan’s lemma to yield a finite value for these integrals.
Finally, we put the added imaginary part to zero. This results in
$\mathcal{F}^{-1}\left[1/(\omega-b);t\right]=\Theta(t)2\pi ie^{ibt}$, where
$b$ is real. The positive $i\epsilon$-interpretation is justified, since in
the limit of $\alpha\to 0$, we would like to recover the usual Caldeira-
Leggett model. The physical justification is that the particle is only put in
the bath after time $t=0$. Therefore, the interaction term should not have an
impact before $t=0$. A Heaviside function, which automatically comes out of
the positive $i\epsilon$-interpretation, ensures that, only after $t=0$,
coupling between the bath and the particle occurs222The first term in Eq. (70)
still shows a coupling between the bath and the particle, even for $t<0$.
However, this term will exactly cancel the term coming from
$\mathcal{L}_{CT}$, such that it has no contribution in the equation of motion
of the particle.. Substituting Eq. (71) into Eq. (70) results in

$\displaystyle q_{j}^{p}(t)$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\left[\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)-\frac{1}{\omega_{j}}\frac{d^{2}}{dt^{2}}\sin(\omega_{j}t)\Theta(t)\ast\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)\right]$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\Bigg{\\{}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)-\frac{d}{dt}\left[\cos(\omega_{j}t)\Theta(t)+\frac{1}{\omega_{j}}\sin(\omega_{j}t)\delta(t)\right]\ast\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)\Bigg{\\}}$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\left[\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)-\frac{d}{dt}\cos(\omega_{j}t)\Theta(t)\ast\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)\right]$
$\displaystyle=\frac{C_{j}}{m_{j}\omega_{j}^{2}}\left\\{\prescript{W}{+}{\bm{D}_{t}^{\alpha}}Q(t)-\frac{d}{dt}\prescript{W}{+}{\bm{D}_{t}^{\alpha}}\left[\cos(\omega_{j}t)\Theta(t)\ast
Q(t)\right]\right\\},$ (72)

where in the first line we acted with a derivative on
$\sin(\omega_{j}t)\Theta(t)$, in the second line we used the fact that
$\frac{d}{dt}\Theta(t)=\delta(t)$ which results in a term of the form
$\delta(t)\sin(\omega_{j}t)\ast Q(t)\sim\sin(0)=0$, and in the last line we
applied Eq. 65 to switch the order of the convolution and fractional
derivative.

Figure 3: The pole structure of the function $1/(z-i\epsilon)$ is depicted in
the above figure. After applying complex analysis and taking $\epsilon\to 0$,
one can give a finite value to a divergent integral over $1/z$, which we call
the positive $i\epsilon$-interpretation. Note that the shown contour is only
chosen for $t>0$.

##  Appendix C The Mellin transform

The Mellin transform is defined by

$\displaystyle\mathcal{M}\left[f(t);s\right]=\int_{0}^{\infty}dtf(t)t^{s-1}.$
(73)

In general, the integral exists only for complex values of $s=a+ib$, such that
$a_{1}<a<a_{2}$, where $a_{1}$ and $a_{2}$ depend on the function $f(t)$ to
transform. This introduces what is called the strip of definition of the
Mellin transform, which will be denoted by $S(a_{1},a_{2})$. The Mellin
transform of the exponential function is, e.g., given by

$\displaystyle\mathcal{M}\left[e^{-at};s\right]$
$\displaystyle=a^{-s}\Gamma(s),$ (74)

which converges even for complex values of $a$ [62] and Re$(s)>0$.
Investigation of the Mellin transform of a more difficult function, namely
$1/\left(1+ct^{a}\right)$, is done by starting with the following equation

$\displaystyle\int_{0}^{\infty}du\frac{u^{\frac{s}{a}-1}}{1+u}=\Gamma\left(\frac{s}{a}\right)\Gamma\left(1-\frac{s}{a}\right),$
(75)

which can be derived from the beta function. Next, we apply the substitution
$u=ct^{a}$ on the left-hand side, which results in333Note that the modulus is
due to the reversal in the direction of integration for a negative power.

$\displaystyle\int_{0}^{\infty}dt\frac{t^{s-1}}{1+ct^{a}}=\frac{1}{|a|}c^{-\frac{s}{a}}\Gamma\left(\frac{s}{a}\right)\Gamma\left(1-\frac{s}{a}\right).$
(76)

Finally, we can identify the Mellin transform of the function $f(t)$ on the
left-hand side of equation (76), to conclude that

$\displaystyle\mathcal{M}\left[\frac{1}{1+ct^{a}};s\right]=\frac{1}{|a|}c^{-\frac{s}{a}}\Gamma\left(\frac{s}{a}\right)\Gamma\left(1-\frac{s}{a}\right),$
(77)

where the region of convergence is given by $0<\text{Re}(s)/a<1$, see Ref.
[62]. The inverse Mellin transform of a function, $F(s)$, defined on a strip
$S(a_{1},a_{2})$ (see Fig. 4), is given by

$\displaystyle f(t)=\mathcal{M}^{-1}\left[F(s);t\right]=\frac{1}{2\pi
i}\int_{a-i\infty}^{a+i\infty}dsF(s)t^{-s},$ (78)

where $a_{1}<a<a_{2}$.

Figure 4: A depiction of the fundamental strip, needed for the inverse Mellin
transform. Here, we chose $a_{1}=1.3$, $a=2.7$ and $a_{2}=3.9$. This means
that the fundamental strip in the figure is $S(1.3,3.9)$.

##  Appendix D The Fox H-function

In this appendix, we briefly introduce the Fox H-function, see Refs. [63, 64,
65]. The usefulness of the Fox H-function is that there are known asymptotic
expansions. Therefore, when we derive an equality between Eq. (46) and a Fox
H-function, we also find the asymptotic expansions of the original integral.

Let $m,n,p,q\in\mathbb{N}$ with $0\leq n\leq p$ and $1\leq m\leq q$,
$A_{i},B_{j}\in\mathbb{R}_{+}$, $a_{i},b_{j}\in\mathbb{C}$ for $i=1,...,p$ and
$j=1,...,q$. The Fox H-function with parameters
$m,n,p,q,a_{i},A_{i},b_{j},B_{j}$ and argument $z\neq 0$ is defined by

$\displaystyle
H^{m,n}_{p,q}\left[z\middle|\begin{array}[]{ccc}\left(a_{1},A_{1}\right),&...,&\left(a_{p},A_{p}\right)\\\
\left(b_{1},B_{1}\right),&...,&\left(b_{q},B_{q}\right)\end{array}\right]=\frac{1}{2\pi
i}\int_{L}ds\Upsilon(s)z^{-s},$ (81)

where $z^{-s}$ is not necessarily the principal value complex power function
and $\Upsilon(s)$ is given by

$\displaystyle\Upsilon(s)=\frac{\prod_{i=1}^{m}\Gamma(b_{i}+sB_{i})\prod_{l=1}^{n}\Gamma(1-a_{l}-sA_{l})}{\prod_{i^{\prime}=m+1}^{q}\Gamma(1-b_{i^{\prime}}-sB_{i^{\prime}})\prod_{l^{\prime}=n+1}^{p}\Gamma(a_{l^{\prime}}+sA_{l^{\prime}})}.$
(82)

Here, an empty product is always interpreted as unity and $L$ is a suitable
contour, which separates the poles of the two factors in the numerator of
$\Upsilon(s)$. For more information about the contour, $L$, see Ref. [63]. We
will also denote the Fox H-function as

$\displaystyle H^{m,n}_{p,q}(z)\equiv
H^{m,n}_{p,q}\left[z\middle|\begin{array}[]{ccc}\left(a_{1},A_{1}\right),&...,&\left(a_{p},A_{p}\right)\\\
\left(b_{1},B_{1}\right),&...,&\left(b_{q},B_{q}\right)\end{array}\right]$
(85)

when the specification of the parameters, $a_{l},A_{l},b_{i},B_{i}$, is not
necessary. Let us now consider a simple example, namely, how the exponential
function can be expressed in terms of the Fox H-function. We start with

$\displaystyle H^{1,0}_{0,1}\left[z\middle|\begin{array}[]{c}\\\
\left(0,1\right)\end{array}\right]$ $\displaystyle=\frac{1}{2\pi
i}\int_{L}ds\Gamma(s)z^{-s}.$ (88)

Now, we choose the contour L to be the vertical line, $\text{Re}(s)=\gamma$,
$\gamma>0$, which lies to the right of the poles of $\Gamma(z)$. Thus, we have

$\displaystyle H^{1,0}_{0,1}\left[z\middle|\begin{array}[]{c}\\\
\left(0,1\right)\end{array}\right]$ $\displaystyle=\frac{1}{2\pi
i}\int_{\gamma-i\infty}^{\gamma+i\infty}ds\Gamma(s)z^{-s}.$ (91)

Next, we assume that we can close the contour such that it picks up all the
poles of the gamma function444A more formal treatment of this example can be
found in Ref. [66]. In the same reference, more advanced examples are
presented.. Applying the residue theorem to the poles of the gamma function
yields

$\displaystyle H^{1,0}_{0,1}\left[z\middle|\begin{array}[]{c}\\\
\left(0,1\right)\end{array}\right]$
$\displaystyle=\sum_{\nu=0}^{\infty}z^{\nu}\frac{(-1)^{\nu}}{\nu!},$ (94)
$\displaystyle=e^{-z},$ (95)

which shows that

$\displaystyle e^{-z}=H^{1,0}_{0,1}\left[z\middle|\begin{array}[]{c}\\\
\left(0,1\right)\end{array}\right].$ (98)

Some quantities that will be important to determine the asymptotic behavior of
Fox-H functions are defined as

$\displaystyle\mu_{H}\equiv\sum_{j=1}^{q}B_{j}-\sum_{j=1}^{p}A_{j},$ (99)
$\displaystyle\delta_{H}\equiv\sum_{j=1}^{q}b_{j}-\sum_{j=1}^{p}a_{j}+\frac{p-q}{2},$
(100)
$\displaystyle\alpha_{H}\equiv\sum_{j=1}^{n}A_{j}-\sum_{j=n+1}^{p}A_{j}+\sum_{j=1}^{m}B_{j}-\sum_{j=m+1}^{q}B_{j}.$
(101)

Now, we can use Section 1.7 from Ref. [63], which shows the following
properties.

  1. (i)

If $\alpha_{H}>0$ and $\left|\text{arg}(z)\right|<\frac{1}{2}\pi\alpha_{H}$,
then the Fox H-function has the asymptotic expansions at zero given by

$\displaystyle H^{m,n}_{p,q}(z)$ $\displaystyle=\mathcal{O}(z^{c}),$
$\displaystyle|z|\to 0$ (102)

for $\mu_{H}\geq 0$ and

$\displaystyle H^{m,n}_{p,q}(z)$
$\displaystyle=\mathcal{O}\left(z^{c}\left|\ln(z)\right|^{N-1}\right),$
$\displaystyle|z|\to 0$ (103)

for $\mu_{H}<0$, with $c=\min_{1\leq i\leq
m}\left[\text{Re}(b_{i})/B_{i}\right]$ and $N$ is the order of the poles of
$\Gamma(b_{i}+sB_{i})$, to which some other poles of another gamma function
$\Gamma(b_{k}+B_{k}s)$ with $1\leq i\leq m$ and $1\leq k\leq m$ could
coincide.

  2. (ii)

If $\mu_{H}<0$, $\alpha_{H}=0$, then the asymptotic expansion near zero
becomes

$\displaystyle H^{m,n}_{p,q}(z)=\mathcal{O}(z^{\sigma}),\qquad|z|\to
0,\>|\text{arg}(z)|\leq\widetilde{\epsilon},$ (104)

where $\widetilde{\epsilon}$ is a constant with
$0<\widetilde{\epsilon}<(\pi/2)\min_{1\leq i\leq m;m+1\leq k\leq
q}\left(A_{i},B_{k}\right)$ and $\sigma=\min_{1\leq i\leq
m}\big{\\{}\text{Re}(b_{i})/B_{i},[\text{Re}(\delta_{H})+1/2]/\mu_{H}\big{\\}}$.

  3. (iii)

If $\alpha_{H}>0$, then the Fox H-function has the asymptotic expansion at
infinity given by

$\displaystyle H^{m,n}_{p,q}(z)$ $\displaystyle=\mathcal{O}(z^{d}),$
$\displaystyle|z|\to\infty$ (105)

for $\mu_{H}\leq 0$ and

$\displaystyle H^{m,n}_{p,q}(z)$
$\displaystyle=\mathcal{O}(z^{d}\left|\ln(z)\right|^{M-1}),$
$\displaystyle|z|\to\infty$ (106)

for $\mu_{H}>0$, with $d=\min_{1\leq l\leq
n}\left[(\text{Re}(a_{l})-1)/A_{l}\right]$ and $M$ is the order of the poles
$\Gamma\left(1-a_{l}-sA_{l}\right)$ to which some of the poles of
$\Gamma\left(1-a_{k}-sA_{k}\right)$ with $1\leq l\leq n$ and $1\leq k\leq n$
could coincide.

  4. (iv)

If $\mu_{H}>0$ and $\alpha_{H}=0$, then

$\displaystyle
H^{m,n}_{p,q}(z)=\mathcal{O}(z^{\rho}),\qquad|z|\to\infty,\>\left|\text{arg}(z)\right|\leq\epsilon,$
(107)

where $\epsilon$ is a constant with
$0<\epsilon<(\pi/2)\min_{n+1\leq j\leq p;1\leq k\leq
m}\left(A_{j},B_{k}\right)$ and
$\rho=\max_{1\leq l\leq
n}\big{\\{}[\text{Re}(a_{l})-1]/A_{l},[\text{Re}(\delta_{H})+1/2]/\mu_{H}\big{\\}}$.

##  Appendix E Asymptotic expansion of the MSD

Taking the Mellin transform of $I^{\prime}(t)$, see Eq. (47), we get

$\displaystyle F(s)$
$\displaystyle\equiv\mathcal{M}\left[I^{\prime}(t);s\right]$
$\displaystyle=\int_{0}^{\infty}dtt^{s-1}\int_{-\infty}^{\infty}d\omega
e^{i\omega t}\frac{(i\omega)^{\alpha}}{i\eta|\omega|^{2\alpha}-M\omega}$
$\displaystyle=\frac{1}{M}\Gamma(s)\int_{-\infty}^{\infty}d\omega\frac{\left(i\omega\right)^{\alpha}\left(-i\omega\right)^{-s}}{i\frac{\eta}{M}|\omega|^{2\alpha}-\omega},$
(108)

where in the last line we used Eq. 74. This already gives a restriction on
$s$, namely $\text{Re}(s)>0$. Next, we apply properties of the complex power
function, together with the fact that $\omega\in\mathbb{R}$, which results in

$\displaystyle F(s)$
$\displaystyle=\frac{1}{M}\Gamma(s)\int_{-\infty}^{\infty}d\omega\frac{|\omega|^{\alpha-s}e^{\text{arg}(i\omega)i\alpha-\text{arg}(-i\omega)is}}{i\frac{\eta}{M}|\omega|^{2\alpha}-\omega}$
$\displaystyle=\frac{1}{M}\Gamma(s)\int_{-\infty}^{\infty}d\omega\frac{|\omega|^{\alpha-s}e^{i\text{arg}(i\omega)\left(\alpha+s\right)}}{i\frac{\eta}{M}|\omega|^{2\alpha}-\omega}.$
(109)

We can further simplify the above equation by splitting the integral in its
positive and negative parts. Doing the substitution $\omega\to-\omega$ in the
latter yields

$\displaystyle F(s)$
$\displaystyle=\frac{1}{M}\Gamma(s)\left[e^{i\frac{\pi}{2}\left(\alpha+s\right)}\int_{0}^{\infty}d\omega\frac{\omega^{\alpha-s}}{i\frac{\eta}{M}\omega^{2\alpha}-\omega}+e^{-i\frac{\pi}{2}\left(\alpha+s\right)}\int_{0}^{\infty}d\omega\frac{\omega^{\alpha-s}}{i\frac{\eta}{M}\omega^{2\alpha}+\omega}\right]$
$\displaystyle=\frac{1}{M}\Gamma(s)\left[e^{-i\frac{\pi}{2}\left(\alpha+s\right)}\int_{0}^{\infty}d\omega\frac{\omega^{\alpha-s-1}}{1+i\frac{\eta}{M}\omega^{2\alpha-1}}-e^{i\frac{\pi}{2}\left(\alpha+s\right)}\int_{0}^{\infty}d\omega\frac{\omega^{\alpha-s-1}}{1-i\frac{\eta}{M}\omega^{2\alpha-1}}\right]$
$\displaystyle=\frac{1}{M|2\alpha-1|}\left[e^{-i\frac{\pi}{2}(\alpha+s)}\left(i\frac{\eta}{M}\right)^{\frac{s-\alpha}{2\alpha-1}}-e^{i\frac{\pi}{2}(\alpha+s)}\left(-i\frac{\eta}{M}\right)^{\frac{s-\alpha}{2\alpha-1}}\right]\Gamma(s)\Gamma\left(\frac{\alpha-s}{2\alpha-1}\right)\Gamma\left(1-\frac{\alpha-s}{2\alpha-1}\right)$
$\displaystyle=\frac{1}{M|2\alpha-1|}\left(\frac{\eta}{M}\right)^{\frac{s-\alpha}{2\alpha-1}}\left[e^{-i\frac{\pi}{2}\left(s+\frac{2\alpha^{2}}{2\alpha-1}\right)}-e^{i\frac{\pi}{2}\left(s+\frac{2\alpha^{2}}{2\alpha-1}\right)}\right]\Gamma(s)\Gamma\left(\frac{\alpha-s}{2\alpha-1}\right)\Gamma\left(1-\frac{\alpha-s}{2\alpha-1}\right)$
$\displaystyle=-\frac{2i}{M|2\alpha-1|}t_{s}^{s-\alpha}\sin\left[\frac{\pi}{2}\left(s+\frac{2\alpha^{2}}{2\alpha-1}\right)\right]\Gamma(s)\Gamma\left(\frac{\alpha-s}{2\alpha-1}\right)\Gamma\left(1-\frac{\alpha-s}{2\alpha-1}\right),$
(110)

where in the third line we applied Eq. (77) to each term separately and in the
last line we identified the characteristic time-scale of the system given by
$t_{s}=(M/\eta)^{1/(1-2\alpha)}$. However, using this example restricts the
values of $s$ even further. The additional constraint reads
$0<\text{Re}\left(\frac{\alpha-s}{2\alpha-1}\right)<1$. Together with the
previous constraint that we already mentioned, $\text{Re}(s)>0$, we obtain
three different regions of convergence that are important to define the
inverse Mellin transform, see Eq. (78). Just like in the definition, we call
them fundamental strips (FS) and they are given by

$\text{FS}=\begin{cases}&S(0,\alpha)\qquad\ \ \ \ \ \ \text{ if }\alpha>1,\\\
&S(1-\alpha,\alpha)\qquad\>\text{ if }\frac{1}{2}<\alpha<1,\\\
&S(0,1-\alpha)\qquad\ \text{ if }\alpha<\frac{1}{2}.\end{cases}$ (111)

Finally, we use Euler’s reflection formula to get

$\displaystyle F(s)=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{s-\alpha}\frac{\Gamma(s)\Gamma\left(\frac{\alpha-s}{2\alpha-1}\right)\Gamma\left(1-\frac{\alpha-s}{2\alpha-1}\right)}{\Gamma\left(\frac{\alpha^{2}}{2\alpha-1}+\frac{s}{2}\right)\Gamma\left(1-\frac{\alpha^{2}}{2\alpha-1}-\frac{s}{2}\right)},$
(112)

which completes the derivation of Eq. (48). In the case of the first and
second fundamental strip, where $2\alpha-1>0$, we rewrite Eq. (112) as

$\displaystyle F(s)=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{s-\alpha}\frac{\Gamma(s)\Gamma\left(1-\frac{\alpha-1}{2\alpha-1}-\frac{s}{2\alpha-1}\right)\Gamma\left(\frac{\alpha-1}{2\alpha-1}+\frac{s}{2\alpha-1}\right)}{\Gamma\left(\frac{\alpha^{2}}{2\alpha-1}+\frac{s}{2}\right)\Gamma\left(1-\frac{\alpha^{2}}{2\alpha-1}-\frac{s}{2}\right)}.$
(113)

Finally, we take the inverse Mellin transform, defined by Eq. (78), which
results in

$\displaystyle I^{\prime}(t)$
$\displaystyle=\mathcal{M}^{-1}\left[F(s);t\right]=\frac{1}{2\pi
i}\int_{a-i\infty}^{a+i\infty}dsF(s)t^{-s},$ (114)

where $F(s)$ is given by Eq. (113). Comparing the above equation with Eqs.
(81) and (82), we get

$\displaystyle I^{\prime}(t)=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{-\alpha}H^{21}_{23}\left[\frac{t}{t_{s}}\
\middle|\begin{array}[]{ccc}&\left(\frac{\alpha-1}{2\alpha-1},\frac{1}{2\alpha-1}\right),&\left(\frac{\alpha^{2}}{2\alpha-1},\frac{1}{2}\right)\\\
(0,1)&\left(\frac{\alpha-1}{2\alpha-1},\frac{1}{2\alpha-1}\right),&\left(\frac{\alpha^{2}}{2\alpha-1},\frac{1}{2}\right)\end{array}\right],$
(117)

which is exactly Eq. (49). Note that Eq. (117) holds for both the first and
second fundamental strips. Next, we focus on the third fundamental strip, and
rewrite Eq. (112) as

$\displaystyle F(s)=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{s-\alpha}\frac{\Gamma(s)\Gamma\left(\frac{\alpha}{2\alpha-1}+\frac{s}{1-2\alpha}\right)\Gamma\left(1-\frac{\alpha}{2\alpha-1}-\frac{s}{1-2\alpha}\right)}{\Gamma\left(\frac{\alpha^{2}}{2\alpha-1}+\frac{s}{2}\right)\Gamma\left(1-\frac{\alpha^{2}}{2\alpha-1}-\frac{s}{2}\right)},$
(118)

which, when taking the inverse Mellin transform and comparing with Eqs. (81)
and (82) results in

$\displaystyle I^{\prime}(t)$ $\displaystyle=-\frac{2\pi
i}{M|2\alpha-1|}t_{s}^{-\alpha}H^{21}_{23}\left[\frac{t}{t_{s}}\
\middle|\begin{array}[]{ccc}&\left(\frac{\alpha}{2\alpha-1},\frac{1}{1-2\alpha}\right),&\left(\frac{\alpha^{2}}{2\alpha-1},\frac{1}{2}\right)\\\
(0,1)&\left(\frac{\alpha}{2\alpha-1},\frac{1}{1-2\alpha}\right),&\left(\frac{\alpha^{2}}{2\alpha-1},\frac{1}{2}\right)\end{array}\right].$
(121)

What remains is to apply the known asymptotic expansions for the Fox
H-function to Eqs. (117) and (121). First, we note that in both cases we have
$\mu_{H},\alpha_{H}>0$. This means that the short and long time expansion are
given by Eqs. (102) and (106), respectively. The short-time behavior is given
by

$\displaystyle I^{\prime}(t)\stackrel{{\scriptstyle t\to
0}}{{\sim}}\begin{cases}&1\qquad\quad\ \ \ \text{ for }\alpha>1,\\\
&t^{\alpha-1}\qquad\ \>\text{ for }\frac{1}{2}<\alpha<1,\\\
&t^{\frac{2\alpha^{2}}{2\alpha-1}}\qquad\>\text{ for
}\frac{1}{4}<\alpha<\frac{1}{2},\\\ &t^{-\alpha}\qquad\ \ \ \text{ for
}0<\alpha\leq\frac{1}{4},\\\ &t^{\frac{2\alpha^{2}}{2\alpha-1}}\qquad\ \text{
for }\alpha<0,\\\ \end{cases}$ (122)

and the long-time behavior is

$\displaystyle I^{\prime}(t)\stackrel{{\scriptstyle
t\to\infty}}{{\sim}}\begin{cases}&t^{-\alpha}\quad\quad\ \>\text{ for
}\alpha>1,\\\ &t^{\alpha-1}\quad\>\quad\text{ for }\alpha<1,\alpha\neq
0.\end{cases}$ (123)

From now on we will discard the case where
$\alpha=-1/2,(-1\pm\sqrt{3})/2,(-3\pm\sqrt{21})/4$. This is because we want to
keep the calculation as simple as possible and therefore do not want any
logarithms while integrating the expansions. We thus have that
$\alpha\in\mathbb{R}\setminus\mathcal{P}$ with
$\mathcal{P}=\big{\\{}(-3-\sqrt{21})/4,(-1-\sqrt{3})/2,-1/2,0,(1+\sqrt{3})/2,(-3+\sqrt{21})/4,1/2,1\big{\\}}$.
Substituting the expansions given by Eqs. (122) and (123) into Eq. (47) adds a
one to the exponents of the asymptotic expansions due to the integration555To
calculate the integral over the small time asymptotic expansion of
$I^{\prime}(t)$ which results in the asymptotic expansions of $I(t)$, one has
to technically introduce a small time cut-off since otherwise the evaluation
at $t=0$ would diverge.. Thus, we obtain

$\displaystyle I(t)\stackrel{{\scriptstyle t\to
0}}{{\sim}}\begin{cases}&t\qquad\quad\ \ \>\>\quad\text{ for }\alpha>1,\\\
&t^{\alpha}\qquad\ \>\qquad\text{ for }\frac{1}{2}<\alpha<1,\\\
&t^{\frac{2\alpha^{2}}{2\alpha-1}+1}\qquad\>\text{ for
}\frac{1}{4}<\alpha<\frac{1}{2},\\\ &t^{-\alpha+1}\qquad\ \ \ \text{ for
}0<\alpha\leq\frac{1}{4},\\\ &t^{\frac{2\alpha^{2}}{2\alpha-1}+1}\qquad\
\text{ for }\alpha<0,\\\ \end{cases}$ (124)

where $\alpha\notin\mathcal{P}$. Finally, we can substitute these asymptotic
expansions into Eq. (45) and Taylor expand the exponential factor up to first
order, to get the small time behavior of the MSD, which is given by

$\displaystyle\left\langle Q(t)^{2}\right\rangle\stackrel{{\scriptstyle t\to
0}}{{\sim}}\begin{cases}&t^{2}\qquad\text{ for
}\alpha>\frac{-1-\sqrt{5}}{4},\alpha\notin\mathcal{P},\\\
&t^{2\alpha+3}\quad\text{for
}\alpha<\frac{-1-\sqrt{5}}{4},\alpha\notin\mathcal{P}.\end{cases}$ (125)

We can do exactly the same for the long-time behavior, except that now we
discard the exponential in Eq. (45) because it vanishes as $t\to\infty$. The
long-time behavior of the MSD is

$\displaystyle\left\langle Q(t)^{2}\right\rangle\stackrel{{\scriptstyle
t\to\infty}}{{\sim}}\begin{cases}&t^{-2\alpha+3}\quad\ \text{ for
}\alpha>1,\alpha\notin\mathcal{P},\\\ &t^{2\alpha+1}\quad\ \ \ \text{ for
}\alpha<1,\alpha\notin\mathcal{P}.\end{cases}$ (126)

## References

  * Bao  _et al._ [2021] C. Bao, M.-G. Suh, B. Shen, K. Şafak, A. Dai, H. Wang, L. Wu, Z. Yuan, Q.-F. Yang, A. B. Matsko, F. X. Kärtner, and K. J. Vahala, Quantum diffusion of microcavity solitons, Nat. Phys. 17, 462 (2021).
  * Kagan [1992] Y. Kagan, Quantum diffusion in solids, J. Low Temp. Phys. 87, 525 (1992).
  * Kleinert [2010] P. Kleinert, Theory of hot-electron quantum diffusion in semiconductors, Phys. Rep. 485, 1 (2010).
  * Storchak and Prokof’ev [1998] V. G. Storchak and N. V. Prokof’ev, Quantum diffusion of muons and muonium atoms in solids, Rev. Mod. Phys. 70, 929 (1998).
  * Guarneri [1989] I. Guarneri, Spectral properties of quantum diffusion on discrete lattices, EPL 10, 95 (1989).
  * Biagetti  _et al._ [2018] M. Biagetti, G. Franciolini, A. Kehagias, and A. Riotto, Primordial black holes from inflation and quantum diffusion, J. Cosmol. Astropart. Phys. 2018 (07), 032.
  * Sagi  _et al._ [2012] Y. Sagi, M. Brook, I. Almog, and N. Davidson, Observation of anomalous diffusion and fractional self-similarity in one dimension, Phys. Rev. Lett. 108, 093002 (2012).
  * Regner  _et al._ [2013] B. M. Regner, D. Vučinić, C. Domnisoru, T. M. Bartol, M. W. Hetzer, D. M. Tartakovsky, and T. J. Sejnowski, Anomalous diffusion of single particles in cytoplasm, Biophys. J. 104, 1652 (2013).
  * Banks and Fradin [2005] D. S. Banks and C. Fradin, Anomalous diffusion of proteins due to molecular crowding, Biophys. J. 89, 2960 (2005).
  * Einstein [1905] A. Einstein, Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen, Ann. Phys. (Leipzig) 322, 549 (1905).
  * Metzler and Klafter [2000] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339, 1 (2000).
  * Chen [2006] W. Chen, Time–space fabric underlying anomalous diffusion, Chaos Solitons Fractals 28, 923 (2006).
  * Xu  _et al._ [2021] X.-Y. Xu, X.-W. Wang, D.-Y. Chen, C. M. Smith, and X.-M. Jin, Quantum transport in fractal networks, Nat. Photonics 15, 703 (2021).
  * Caldeira [2014] A. O. Caldeira,  _An Introduction to Macroscopic Quantum Phenomena and Quantum Dissipation_ (Cambridge University Press, 2014).
  * Weiss [2012] U. Weiss,  _Quantum dissipative systems_ (World Scientific, 2012).
  * Caldeira and Leggett [1983] A. Caldeira and A. Leggett, Quantum tunnelling in a dissipative system, Ann. Phys. (NY) 149, 374 (1983).
  * Leggett  _et al._ [1987] A. J. Leggett, S. Chakravarty, A. T. Dorsey, M. P. Fisher, A. Garg, and W. Zwerger, Dynamics of the dissipative two-state system, Rev. Mod. Phys. 59, 1 (1987).
  * Caldeira and Leggett [1981] A. O. Caldeira and A. J. Leggett, Influence of dissipation on quantum tunneling in macroscopic systems, Phys. Rev. Lett 46, 211 (1981).
  * Ferialdi [2017] L. Ferialdi, Dissipation in the Caldeira-Leggett model, Phys. Rev. A 95, 052109 (2017).
  * Lutz [2001] E. Lutz, Fractional Langevin equation, Phys. Rev. E 64, 051106 (2001).
  * Verstraten  _et al._ [2021] R. C. Verstraten, R. F. Ozela, and C. Morais Smith, Time glass: A fractional calculus approach, Phys. Rev. B 103, L180301 (2021).
  * Burov and Barkai [2008] S. Burov and E. Barkai, Critical exponent of the fractional Langevin equation, Phys. Rev. Lett. 100, 070601 (2008).
  * Verstraten [2020] R. Verstraten,  _The fractional Langevin equation_ , Master’s thesis, Utrecht University, Heidelberglaan 8, 3584 CS Utrecht (2020).
  * Oldham and Spanier [1974] K. Oldham and J. Spanier,  _The fractional calculus theory and applications of differentiation and integration to arbitrary order_ (Elsevier, 1974).
  * Mainardi [2018] F. Mainardi, Fractional calculus: Theory and applications (2018).
  * de Oliveira and Machado [2014] E. C. de Oliveira and J. A. T. Machado, A review of definitions for fractional derivatives and integral, Math. Probl. Eng. https://doi.org/10.1155/2014/238459 (2014).
  * Bayın [2016] S. Ş. Bayın, Definition of the Riesz derivative and its application to space fractional quantum mechanics, J. Math. Phys. 57, 123501 (2016).
  * Ferrari [2018] F. Ferrari, Weyl and Marchaud derivatives: A forgotten history, Mathematics 6, 6 (2018).
  * Samko  _et al._ [1993] S. G. Samko, A. A. Kilbas, O. I. Marichev, _et al._ ,  _Fractional integrals and derivatives_ , Vol. 1 (Gordon and breach Science publishers, Yverdon Yverdon-les-Bains, Switzerland, 1993).
  * Miller [2006] K. S. Miller, The Weyl fractional calculus, in  _Fractional Calculus and Its Applications: Proceedings of the International Conference Held at the University of New Haven, June 1974_ (Springer, 2006) pp. 80–89.
  * Mulder  _et al._ [2021] R. A. Mulder, M. A. Caracanhas, and C. M. Smith, Quantizing Lévy flights, Phys. Rev. B 103, 174301 (2021).
  * Agrawal [2002] O. Agrawal, Formulation of Euler–Lagrange equations for fractional variational problems, J. Math. Anal. Appl. 272, 368 (2002).
  * Barkai [2001] E. Barkai, Fractional Fokker-Planck equation, solution, and application, Phys. Rev. E 63, 046118 (2001).
  * Magdziarz  _et al._ [2009] M. Magdziarz, A. Weron, K. Burnecki, and J. Klafter, Fractional Brownian motion versus the continuous-time random walk: A simple test for subdiffusive dynamics, Phys. Rev. Lett. 103, 180602 (2009).
  * Razminia  _et al._ [2014] K. Razminia, A. Razminia, and J. T. Machado, Analysis of diffusion process in fractured reservoirs using fractional derivative approach, Commun. Nonlinear Sci. Numer. Simul. 19, 3161 (2014).
  * Chen  _et al._ [2010] W. Chen, H. Sun, X. Zhang, and D. Korošak, Anomalous diffusion modeling by fractal and fractional derivatives, Comput. Math. with Appl. 59, 1754 (2010).
  * Tarasov [2013] V. E. Tarasov, Review of some promising fractional physical models, Int. J. Mod. Phys. B 27, 1330005 (2013).
  * Umarov [2015] S. Umarov,  _Introduction to Fractional and Pseudo-Differential Equations with Singular Symbols_ , Vol. 41 (2015).
  * Lim  _et al._ [2008] S. Lim, M. Li, and L. Teo, Langevin equation with two fractional orders, Phys. Lett. A 372, 6309 (2008).
  * Bagarello [2020] F. Bagarello, Fourier transforms, fractional derivatives, and a little bit of quantum mechanics, Rocky Mt. J. Math. 50, 415 (2020).
  * Haber [2011] H. Haber, The complex logarithm, exponential and power functions (2011), lecture Notes, University of California.
  * Musina and Nazarov [2016] R. Musina and A. I. Nazarov, On fractional Laplacians – 2, Ann. Inst. Henri Poincare (C) Anal. Non Lineaire 33, 1667 (2016).
  * Lischke  _et al._ [2020] A. Lischke, G. Pang, M. Gulian, F. Song, C. Glusa, X. Zheng, Z. Mao, W. Cai, M. M. Meerschaert, M. Ainsworth, and G. E. Karniadakis, What is the fractional Laplacian? a comparative review with new results, J. Comput. Phys. 404, 109009 (2020).
  * Marino [1993] E. Marino, Quantum electrodynamics of particles on a plane and the Chern-Simons theory, Nucl. Phys. B 408, 551 (1993).
  * Marino  _et al._ [2014] E. C. Marino, L. O. Nascimento, V. S. Alves, and C. Morais Smith, Unitarity of theories containing fractional powers of the d’Alembertian operator, Phys. Rev. D 90, 105003 (2014).
  * Marino  _et al._ [2015] E. Marino, L. O. Nascimento, V. S. Alves, and C. M. Smith, Interaction induced quantum valley Hall effect in graphene, Phys. Rev. X 5, 011040 (2015).
  * Marino  _et al._ [2018] E. C. Marino, L. O. Nascimento, V. S. Alves, N. Menezes, and C. M. Smith, Quantum-electrodynamical approach to the exciton spectrum in transition-metal dichalcogenides, 2D Mater. 5, 041006 (2018).
  * Kirichenko and Stephanovich [2021] E. Kirichenko and V. Stephanovich, The influence of Coulomb interaction screening on the excitons in disordered two-dimensional insulators, Sci. Rep. 11, 1 (2021).
  * Kehrein and Mielke [1996] S. K. Kehrein and A. Mielke, On the spin-boson model with a sub-ohmic bath, Phys. Lett. A 219, 313 (1996).
  * Anders  _et al._ [2007] F. B. Anders, R. Bulla, and M. Vojta, Equilibrium and nonequilibrium dynamics of the sub-ohmic spin-boson model, Phys. Rev. Lett. 98, 210402 (2007).
  * Vojta  _et al._ [2005] M. Vojta, N.-H. Tong, and R. Bulla, Quantum phase transitions in the sub-ohmic spin-boson model: Failure of the quantum-classical mapping, Phys. Rev. Lett. 94, 070604 (2005).
  * Zheng [2004] H. Zheng, Dynamics of a two-level system coupled to ohmic bath: a perturbation approach, Eur. Phys. J. B 38, 559 (2004).
  * Schneider [1986] W. R. Schneider, Stable distributions: Fox function representation and generalization, in  _Stochastic Processes in Classical and Quantum Systems_ , edited by S. Albeverio, G. Casati, and D. Merlini (Springer Berlin Heidelberg, Berlin, Heidelberg, 1986) pp. 497–511.
  * Jeon and Metzler [2010] J.-H. Jeon and R. Metzler, Fractional Brownian motion and motion governed by the fractional Langevin equation in confined geometries, Phys. Rev. E 81, 021103 (2010).
  * Kobelev and Romanov [2000] V. Kobelev and E. Romanov, Fractional Langevin equation to describe anomalous diffusion, Prog. Theor. Phys. 139, 470 (2000).
  * Khalili Golmankhaneh  _et al._ [2022] A. Khalili Golmankhaneh, K. Ali, R. Yilmazer, and M. Kaabar, Local fractal Fourier transform and applications, Comput. methods differ. equ. 10, 595 (2022).
  * Laskin [2000] N. Laskin, Fractional quantum mechanics, Phys. Rev. E 62, 3135 (2000).
  * Tarasov [2008] V. E. Tarasov, Weyl quantization of fractional derivatives, J. Math. Phys. 49, 10.1063/1.3009533 (2008).
  * Baleanu [2008] D. Baleanu, New applications of fractional variational principles, Rep. Math. Phys. 61, 199 (2008).
  * Verstraten  _et al._ [2022] R. C. Verstraten, T. Ludwig, R. A. Duine, and C. M. Smith, The fractional Landau-Lifshitz-Gilbert equation, (2022), arXiv:2211.12889 [cond-mat.mes-hall] .
  * Nalbach and Thorwart [2010] P. Nalbach and M. Thorwart, Ultraslow quantum dynamics in a sub-ohmic heat bath, Phys. Rev. B 81, 054308 (2010).
  * Bateman [1954] H. Bateman,  _Tables of Integral Transforms_ , Vol. 1 (McGraw-Hill Book Company, New York, 1954).
  * A.M. Mathai [2010] H. J. H. A.M. Mathai, Ram Kishore Saxena,  _The H-Function: Theory and Applications_ (2010).
  * Kilbas and Saigo [1999] A. A. Kilbas and M. Saigo, On the H-function, International Journal of Stochastic Analysis 12, 191 (1999).
  * Braaksma [1964] B. L. J. Braaksma, Asymptotic expansions and analytic continuations for a class of Barnes-integrals, Compos. Math. 15, 239 (1962-1964).
  * Koekoek [2021] R. Koekoek, Special functions – hypergeometric functions – Barnes’ integral representation (2021), lecture Notes, TU Delft.