# Anomalous and ultraslow diffusion of a particle driven by power-law-
correlated and distributed-order noises

Z. Tomovski  University of Ostrava, Faculty of Sciences, Department of
Mathematics, 30. Dubna 22701 03 Ostrava, Czech Republic
<EMAIL_ADDRESS>K. Górska<EMAIL_ADDRESS>Institute of
Nuclear Physics, Polish Academy of Sciences, ul.Eljasza-Radzikowskiego 152, PL
31342 Krakow, Poland  T. Pietrzak<EMAIL_ADDRESS>Institute of
Nuclear Physics, Polish Academy of Sciences, ul.Eljasza-Radzikowskiego 152, PL
31342 Krakow, Poland  R. Metzler  Institute of Physics and Astronomy,
University of Potsdam, D-14476 Potsdam-Golm, Germany<EMAIL_ADDRESS>Asia Pacific Center for Theoretical Physics, Pohang 37673, Republic of Korea
T. Sandev<EMAIL_ADDRESS>Research Center for Computer Science and
Information Technologies, Macedonian Academy of Sciences and Arts, Bul. Krste
Misirkov 2, 1000 Skopje, Macedonia  Institute of Physics, Faculty of Natural
Sciences and Mathematics, Ss. Cyril and Methodius University, Arhimedova 3,
1000 Skopje, Macedonia  Institute of Physics and Astronomy, University of
Potsdam, D-14476 Potsdam-Golm, Germany

###### Abstract

We study the generalized Langevin equation approach to anomalous diffusion for
a harmonic oscillator and a free particle driven by different forms of
internal noises, such as power-law-correlated and distributed-order noises
that fulfil generalized versions of the fluctuation-dissipation theorem. The
mean squared displacement and the normalized displacement correlation function
are derived for the different forms of the friction memory kernel. The
corresponding overdamped generalized Langevin equations for these cases are
also investigated. It is shown that such models can be used to describe
anomalous diffusion in complex media, giving rise to subdiffusion,
superdiffusion, ultraslow diffusion, strong anomaly, and other complex
diffusive behaviors.

##  I Introduction

The behavior of a test particle of mass $m$, that is coupled to a thermal bath
of temperature $T$, can be described by Newton’s second law for a particle in
the presence of a deterministic external potential $V(x)$ and a stochastically
varying force $\xi(t)$. When the friction acting on the test particle is given
by the constant value $\gamma_{0}$, the resulting dynamic is described by the
standard Langevin equation for a Brownian particle coffey ; langevin ; Zwanzig
,

$m\ddot{x}(t)+m\gamma_{0}\dot{x}(t)+\frac{dV(x)}{dx}=\xi(t),\quad\dot{x}(t)=v(t),$
(1)

where $x(t)$ is the particle displacement and $v(t)$ is its velocity. The
stochastic force, i.e., the ”noise”, $\xi(t)$ is Gaussian, has zero mean, and
is white such that its autocovariance is $\delta$-correlated,
$\langle\xi(t+t^{\prime})\xi(t)\rangle=mk_{B}T\gamma_{0}$ with thermal energy
$k_{B}T$. The mean-squared displacement (MSD) encoded in the stochastic
equation (1) scales quadratically (”ballistically”) at short times, at which
it is dominated by inertial effects, and then crosses over to a linear time
dependence beyond the characteristic time scale $1/\gamma_{0}$.

In the following we will consider generalizations of the stochastic equation
(1) that include ”memory”, by which we mean non-localities in time. A simple
example giving rise to memory is the Brownian harmonic oscillator (”Ornstein-
Uhlenbeck process”)

$\dot{x}(t)=v(t),\quad m\dot{v}=-m\omega^{2}x(t)-m\gamma_{0}v(t)+\xi_{v}(t),$
(2)

where $\xi_{v}(t)$ is zero-mean white Gaussian noise. The subscript $v$
denotes that this noise appears in the equation for $\dot{v}(t)$. On purpose,
we here use the phase space notation explicitly keeping the velocity $v(t)$ in
the second equation. Suppose that $v(0)=0$, so that we can integrate the
second equation from $t=0$, yielding

$v(t)=\int_{0}^{t}e^{-\gamma_{0}(t-t^{\prime})}\left(-\omega^{2}x(t^{\prime})+\xi_{v}(t^{\prime})/m\right)dt^{\prime}.$
(3)

Substituting this result back into the first equation in (2), we find

$\dot{x}(t)=-\int_{0}^{t}K(t^{\prime})x(t-t^{\prime})dt^{\prime}+\xi_{x}(t),$
(4)

where we defined the memory kernel $K(t)$ and the positional friction
$\xi_{x}(t)$,

$K(t)=\omega^{2}e^{-\gamma_{0}t},\quad\xi_{x}(t)=\frac{1}{m}\int_{0}^{t}e^{-\gamma_{0}t^{\prime}}\xi_{v}(t-t^{\prime})dt^{\prime}.$
(5)

At sufficiently long times $t\ll 1/\gamma_{0}$ (or when we start the process
at $t=-\infty$) the Ornstein-Uhlenbeck process reaches equilibium and
satisfies the fluctuation-dissipation theorem
$\langle\xi_{x}(t)\xi_{x}(t^{\prime})\rangle\sim\langle
x^{2}\rangle_{\mathrm{eq}}K(|t-t^{\prime}|)$ with the thermal value $\langle
x^{2}\rangle_{\mathrm{eq}}=k_{B}T/(m\omega^{2})$.

From our result (4), which is characterized by an exponential memory, one thus
cannot conclude that the underlying Ornstein-Uhlenbeck process is non-
Markovian. However, we showed a general property: namely, when integrating out
Markovian degrees of freedom, memory effects in the resulting equation for the
test particle of interest emerge Zwanzig . Indeed, much more pronounced
memories, such as power-law forms, can be effected by eliminating a quasi-
continuum of Markovian degrees of freedom, see, e.g., goychuk ; kupfer ;
lizana ; Zwanzig .

For a general friction memory kernel $\gamma(t)$ the dynamics of a stochastic
process is described in terms of the generalized Langevin equation (GLE)
Zwanzig

$m\ddot{x}(t)+m\int_{0}^{t}\gamma(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}+\frac{dV(x)}{dx}=\xi(t),\quad\dot{x}(t)=v(t).$
(6)

The second fluctuation-dissipation theorem (FDT) is then valid in a thermal
bath of temperature $T$, where fluctuations and dissipation come from the same
source. The FDT relates the friction memory kernel $\gamma(t)$ with the
correlation function $\xi(t)$ of the random force kubo ; Mainardi et al ;
Zwanzig . The FDT allows one to write

$\langle\xi(t+t^{\prime})\xi(t^{\prime})\rangle=C(t)=k_{B}T\gamma(t).$ (7)

Here the friction memory kernel is assumed to satisfy Desposito_Vinales

$\lim_{t\rightarrow\infty}\gamma(t)=\lim_{s\rightarrow 0}s\hat{\gamma}(s)=0,$
(8)

where the hat denotes the Laplace transform of $\gamma(t)$, i.e.,
$\hat{\gamma}(s)=\mathcal{L}\\{\gamma(t);s\\}=\int_{0}^{\infty}\gamma(t)e^{-st}dt$.
We note that when the noise is external in the sense of Klimontovich klimo ,
i.e., the noise is not provided by a heat bath in a non-equilibrium systems,
the relation (7) does not hold and
$\langle\xi(t)\xi(t^{\prime})\rangle=C(t-t^{\prime})$ is used instead. When we
again consider the harmonic oscillator and use the unit mass $m=1$ (i.e.,
$V(x)=\omega^{2}x^{2}/2$), then (6) can rewritten in the form
Vinales_Desposito

$x(t)=\langle
x(t)\rangle+\int_{0}^{t}G(t-t^{\prime})\xi(t^{\prime})\mathrm{d}t^{\prime},\quad
v(t)=\langle
v(t)\rangle+\int_{0}^{t}g(t-t^{\prime})\xi(t^{\prime})\mathrm{d}t^{\prime},$
(9)

where

$\langle x(t)\rangle=x_{0}\left[1-\omega^{2}I(t)\right]+v_{0}G(t),\quad\langle
v(t)\rangle=v_{0}g(t)-\omega^{2}x_{0}G(t)$ (10)

are the average particle displacement velocity, respectively, given the
initial conditions $x_{0}=x(0)$ and $v_{0}=v(0)$. Moreover, we introduced the
so-called relaxation functions $G(t)$, $I(t)=\int_{0}^{t}G(\xi)d\xi$ and
$g(t)=dG(t)/(dt)$, which in the Laplace space read

$\hat{G}(s)=\frac{1}{s^{2}+s\hat{\gamma}(s)+\omega^{2}},\quad\hat{g}(s)=s\hat{G}(s),\quad\mbox{and}\quad\hat{I}(s)=s^{-1}\hat{G}(s).$
(11)

These functions are used to calculate the following four fundamental
quantities (i) the MSD Desposito_Vinales

$\langle x^{2}(t)\rangle=2k_{B}TI(t)=2k_{B}T\int_{0}^{t}G(\xi)d\xi;$ (12)

(ii) the diffusion coefficient ${\cal D}(t)=(1/2)[d\langle
x^{2}(t)\rangle/(dt)]$ that, due to relation (12) can be expressed as

${\cal D}(t)=k_{B}TG(t)=k_{B}T\frac{d}{dt}I(t),$ (13)

compare the proof of relation (13) for the GLE in the free-force case in
Mainardi et al ; Pottier ;
(iii) the normalized displacement autocorrelation function (DACF), that can be
experimentally measured which under the initial conditions $\langle
x_{0}^{2}\rangle=k_{B}T/\omega^{2}$, $\langle x_{0}v_{0}\rangle=0$, and
$\left\langle\xi(t)x_{0}\right\rangle=0$, can be represented as
sandev_physmaced ; sandev_jmp

$C_{X}=\frac{\langle x(t)x_{0}\rangle}{\langle
x^{2}_{0}\rangle}=1-\omega^{2}I(t)=1-\omega^{2}\int_{0}^{t}G(\xi)d\xi.$ (14)

(Note that different DACFs for fractional GLE are studied in sandev_jmp .)
(iv) the normalized velocity autocorrelation function (VACF) Desposito_Vinales

$C_{V}(t)=\frac{\langle v(t)v_{0}\rangle}{\langle
v^{2}_{0}\rangle}=g(t)=\frac{d}{dt}G(t)=\frac{d^{2}}{dt^{2}}I(t).$ (15)

In the large friction limit, we neglect the inertial term $m\ddot{x}(t)$, and
the resulting overdamped GLE (6) has the form

$\displaystyle\int_{0}^{t}\gamma(t-t^{\prime})\dot{x}(t^{\prime})dt^{\prime}+\frac{dV(x)}{dx}=\xi(t),\quad\dot{x}(t)=v(t).$
(16)

The solution of this overdamped equation is of particular interest due to its
application in modeling the anomalous dynamics of colloidal (micron-sized)
test particles or the longer-time internal motion of proteins. Large friction,
which appears due to the liquid environment means that the acceleration
$\ddot{x}(t)$ is negligible in comparison to the effect of the friction term.
Single-particle tracking of colloidal particles in an optical tweezers trap
franosch ; jeon ; jeon1 ; norre ; selh  or the internal motion of proteins can
be considered as the effective motion in an harmonic potential Kou_Xie ; yang
, the motion is described by the GLE

$\displaystyle\int_{0}^{t}\gamma(t-t^{\prime})\dot{x}(t^{\prime})\mathrm{d}t^{\prime}+m\omega^{2}x(t)=\xi(t),\quad\dot{x}(t)=v(t).$
(17)

The relaxation functions in the overdamped limit read

$\displaystyle
G_{o}(t)=\mathcal{L}^{-1}\Big{[}\frac{1}{s\hat{\gamma}(s)+\omega^{2}};t\Big{]},\quad
g_{o}(t)=\mathcal{L}^{-1}\Big{[}\frac{s}{s\hat{\gamma}(s)+\omega^{2}};t\Big{]},\quad\text{and}\quad
I_{o}(t)=\mathcal{L}^{-1}\Big{[}\frac{s^{-1}}{s\hat{\gamma}(s)+\omega^{2}};t\Big{]}.$
(18)

For the case of a white Gaussian form for the noise $\xi(t)$ the GLE (6)
corresponds to the classical overdamped Ornstein-Uhlenbeck process equation
with friction coefficient $\gamma_{0}$. In the force-free case this further
reduces to the Langevin equation (1) for a force-free Brownian particle. At
times $t\ll 1/\gamma_{0}$ the MSD is then given by $\langle
x^{2}(t)\rangle\sim 2[k_{B}T/(\gamma_{0}m)]t$, and thus the diffusion
coefficient becomes $\mathcal{D}=\lim_{t\to\infty}\langle
x^{2}(t)\rangle/(2t)=k_{B}T/(m\gamma_{0})$, whose physical dimensions are
$[\mathcal{D}]=\mathrm{length}^{2}/\mathrm{time}$. The latter result for the
diffusion coefficient of a Brownian particle in fact represents the Einstein-
Smoluchowski-Sutherland relation einstein ; smolu ; spiecho ; suth .

Different forms for the friction memory kernel, particularly power-law forms
Burov_Barkai ; Desposito_Vinales ; Lutz ; Mainardi et al ; sandev_tomovski_pla
; Wang  and Mittag-Leffler (ML) forms sandev physa ; sandev3 ;
sandev_physmaced ; sandev_jmp ; Vinales_Desposito  have been introduced to
model anomalous diffusion, for which the MSD scales non-linearly in time,

$\langle
x^{2}(t)\rangle=\frac{2\mathcal{D}_{\alpha}}{\Gamma(1+\alpha)}t^{\alpha},$
(19)

where $\mathcal{D}_{\alpha}$ is the generalized diffusion coefficient with
physical dimension
$\left[\mathcal{D}_{\alpha}\right]=\mathrm{length}^{2}/\mathrm{time}^{\alpha}$,
and where $\alpha$ is the anomalous diffusion exponent. We distinguish the
cases of subdiffusion ($0<\alpha<1$) and superdiffusion ($1<\alpha$) Metzler3
. Anomalous diffusion of this power-law form occurs in a multitude of systems
across many scales pt ; cas ; hoefling ; pt1 ; igor_wings ; vilk  and it is
non-universal in the sense that the MSD (19) for a given $\alpha$ may emerge
from a range of different anomalous stochastic processes stasrev ; pccp ; andi
; henrik ; igor_sm . We also mention that anomalous diffusion based on the GLE
with crossovers to a different $\alpha$ exponent or normal diffusion can be
modeled in terms of tempered power-law kernels as introduced in daniel  and
applied to the anomalous diffusion of lipids in bilayer membranes jeon_prl .
Similar crossovers can be achieved in terms of formulations with distributed-
order kernels as discussed below. We also note that non-Gaussian processes
with power-law correlated noise were shown to emerge from a superstatistical
approach based on the GLE jakub ; straeten  and interactions of GLE dynamics
with reflecting boundaries were analyzed in vojta .

In previous work sandev_tomovski_pla  we considered the GLE for a free
particle driven by a mixture of $N$ independent internal white Gausian noises

$\xi(t)=\sum_{i=1}^{N}\alpha_{i}\,\xi_{i}(t),$ (20)

such that each has zero mean, $\langle\xi_{i}(t)\xi_{j}(t^{\prime})\rangle=0$
and correlation

$\displaystyle\langle\xi_{i}(t)\xi_{i}(t^{\prime})\rangle=\delta_{ij}\,\zeta_{i}(t^{\prime}-t),$
(21)

where $\delta_{ij}$ is the Kronecker-$\delta$. The correlation function of the
additive noise $\xi(t)$ is then sandev_tomovski_pla

$\displaystyle\langle\xi(t)\xi(t^{\prime})\rangle=\left<\sum_{i=1}^{N}\alpha_{i}\,\xi_{i}(t)\sum_{j=1}^{N}\alpha_{j}\,\xi_{j}(t^{\prime})\right>=\sum_{i=1}^{N}\alpha_{i}^{2}\,\langle\xi_{i}(t)\xi_{i}(t^{\prime})\rangle.$
(22)

From the second FDT (7) we thus see that the noise fulfils

$\displaystyle\sum_{i=1}^{N}\alpha_{i}^{2}\,\zeta_{i}(t)=k_{B}T\gamma(t),$
(23)

where $\gamma(t)$ is the associated friction memory kernel. In
sandev_tomovski_pla  the GLE with internal noises of Dirac-$\delta$, power-
law, and ML types were analyzed and various different diffusive regimes
obtained. Moreover, it was shown that friction memory kernels of distributed
order can be used to describe ultraslow diffusion with a logarithmic time
dependence of the MSD. In what follows we consider a stochastic harmonic
oscillator driven by a mixture of internal noises, from which we recover the
results for a free particle in the limit of vanishing force constant.

Here we analyze the MSD and DACF for different forms of the friction memory
kernel. In Section II we specify the GLE approach to anomalous diffusion. We
analyze the relaxation functions, MSD, and DACF for Dirac-$\delta$, power-law,
and combinations friction kernels. The corresponding overdamped limits are
analyzed and the force-free limit recovered. Distributed-order friction memory
kernels are considered in Section III. It is shown that such kernels yield
ultraslow diffusion, strong anomaly, and other complex behaviors. The
overdamped motion of a harmonic oscillator driven by distributed-order noises
is investigated in detail. A summary is presented in Section V.

##  II GLE in presence of white and power-law noises

We study the GLE for both white and power-law noises, and their combinations.
At the end of this section we also derive the overdamped limit.

###  II.1 Additive white noises

The simplest case of the GLE is obtained for a test particle connected to a
thermal bath of temperature $T$ effecting $N$ additive internal white Gaussian
and zero-mean noises, equation (20), in which each component fulfils
$\zeta_{i}(t)=\delta(t)$ (i.e., $\hat{\zeta}_{i}(s)=1$). Physically such a
joint noise may stem from an environment with different components. As we
assume that the noise is internal, from the second FDT we conclude that the
friction memory kernel is given by

$\hat{\gamma}(s)=\sum_{i=1}^{N}\frac{\alpha_{i}^{2}}{k_{B}T}$ (24)

in Laplace space. From relation (11), we obtain the relaxation function $G(t)$
in the form

$\displaystyle G(t)=\mathcal{L}^{-1}\left\\{\frac{1}{s^{2}+\kappa
s+\omega^{2}};t\right\\}=\frac{2}{\sqrt{\kappa^{2}-4\omega^{2}}}\exp\left(-\frac{\kappa
t}{2}\right)\sinh\left(\frac{t}{2}\sqrt{\kappa^{2}-4\omega^{2}}\right),$ (25)

where $\kappa=\sum_{i=1}^{N}\alpha_{i}^{2}/(k_{B}T)>2\omega$. The limit
$\omega=0$ leads to the known result $G(t)=[1-\exp(-\kappa t)]/\kappa$
sandev_tomovski_pla , i.e., a constant diffusion coefficient in the linear
time dependence of the MSD on time in the long time limit
($t\rightarrow\infty$),

$\langle x^{2}(t)\rangle\sim
2\mathcal{D}t,\quad\text{where}\quad\mathcal{D}=\frac{\left(k_{B}T\right)^{2}}{\sum_{i=1}^{N}\alpha_{i}^{2}},$
(26)

Writing $\kappa=2\omega_{c}$ in terms of the compound frequency
$\omega_{c}=\sum_{i=1}^{N}\alpha_{i}^{2}/(2k_{B}T)$, then
$\mathcal{D}(t)=t\exp(-\omega_{c}t)$, and the MSD is given by

$\langle
x^{2}(t)\rangle=\omega_{c}^{-2}[1-(1+\omega_{c}t)\exp(-\omega_{c}t)].$ (27)

Therefore, the relaxation function $I(t)$ from (14) with
$I(t)=\omega^{-2}\left[1-C_{X}(t)\right]$ is defined in terms of

$C_{X}(t)=e^{-\kappa
t/2}\left[\cosh\left(\frac{t}{2}\sqrt{\kappa^{2}-4\omega^{2}}\right)+\frac{\kappa}{\sqrt{\kappa^{2}-4\omega^{2}}}\sinh\left(\frac{t}{2}\sqrt{\kappa^{2}-4\omega^{2}}\right)\right].$
(28)

With our notation $\kappa=2\omega_{c}$ it is given by
$C_{X}(t)=(1+\omega_{c}t)\exp(-\omega_{c}t)$. The case $\kappa<2\omega$, i.e.,
the case of underdamped motion, yields in the form

$C_{X}(t)=e^{-\kappa
t/2}\left[\cos\left(\frac{t}{2}\sqrt{4\omega^{2}-\kappa^{2}}\right)+\frac{\kappa}{\sqrt{4\omega^{2}-\kappa^{2}}}\sin\left(\frac{t}{2}\sqrt{4\omega^{2}-\kappa^{2}}\right)\right].$
(29)

Thus, the MSD in the long time limit approaches the equilibrium value $\langle
x^{2}(t)\rangle_{\mathrm{eq}}=2k_{B}T/\omega^{2}$. A graphical representation
of the NDCF $C_{X}(t)$ is given in figure 1. From panel (a) we conclude that
in the overdamped overdamped case the NDCF shows a monotonic decay to zero,
without zero crossings. The underdamped case shows oscillatory behavior of
$C_{X}(t)$ with zero crossings (panels (b) and (c)), where oscillations become
more pronounced for increasing oscillator frequency.

(a)  (b)

(c)

Figure 1: Normalized displacement correlation function for $\kappa=1$: (a)
overdamped motion $\kappa>2\omega$ (28), $\omega=\omega_{c}=1/2$ (solid line),
$\omega=3/8$ (dashed line); $\omega=1/4$ (dot-dashed line); $\omega=1/8$
(dotted line); (b) underdamped motion (29), $\omega=\omega_{c}=1/2$ (solid
line), $\omega=5/8$ (dashed line); $\omega=3/4$ (dot-dashed line);
$\omega=7/8$ (dotted line); (c) underdamped motion (29),
$\omega=\omega_{c}=9/8$ (solid line), $\omega=11/8$ (dashed line);
$\omega=13/8$ (dot-dashed line); $\omega=15/8$ (dotted line).

###### Remark 1

Let us analyze the diffusion coefficient in this case somewhat further. We see
that if we consider a single internal white noise the diffusion coefficient is
$\mathcal{D}_{i}=(k_{B}T)^{2}/\alpha_{i}^{2}$, for $i={1,2,\dots,N}$. From
relation (26), by using the relation between harmonic and arithmetic mean, we
obtain

$\mathcal{D}=\frac{1}{\sum_{i=1}^{N}\frac{1}{\mathcal{D}_{i}}}\leq\frac{\sum_{i=1}^{N}\mathcal{D}_{i}}{N^{2}}.$
(30)

So we see that if the particular diffusion coefficients $\mathcal{D}_{i}$ are
identical and equal to $\bar{\mathcal{D}}$, the diffusion coefficient for $N$
independent internal white Gaussian noise terms scales inversely to $N$,
$\mathcal{D}=\bar{\mathcal{D}}/N$.

If we write $\mathcal{D}_{i}=a^{2}/(2\tau_{i})$ in a random walk-like
notation, where $a^{2}$ represents the squared lattice constant or the
variance of the jump length PDF stasrev ; he ; pccp , we rewrite the compound
diffusion coeffcient as

$\mathcal{D}=\frac{1}{\sum_{i=1}^{N}\frac{2\tau_{i}}{a^{2}}}=\frac{a^{2}}{2\sum_{i=1}^{N}\tau_{i}}=\frac{a^{2}}{2N\langle\tau\rangle}=\frac{\overline{\mathcal{D}}}{N}.$
(31)

In this sense the equality in (30) always holds. We note that the mean time
$\langle\tau\rangle$ may be a misrepresentation of the largest time scale
$\tau_{m}\max_{i}\\{\tau_{i}\\}$, depending on the underlying set
$\\{\tau_{i}\\}$.

###  II.2 Different power-law noises

We now turn to generalizing the results of sandev_tomovski_pla , where the
authors studied the GLE for free particles and $N$ independent internal
noises, by using a harmonic potential, corresonding to the Hookean force
$F(x)=-m\omega^{2}x$, in the GLE (6). To this end we first recall the
considerations in Burov_Barkai , where the authors study the GLE for a
harmonic potential and one internal noise, and then generalize these results.
We mention that in Burov_Barkai  complementary polynomials are used to analyze
the overdamped, underdamped, and critical behaviors of the oscillator. Here,
our investigation is based on multinomial Prabhakar-type functions, ${\cal
E}_{(\vec{\mu}),\beta}(s;\vec{a})$ (see EBazhlekova21  and A). Despite the
somewhat complicated notation, their definition and use in analytic
calculations is in fact quite straightforward.

####  II.2.1 Power-law noises

We first consider the single-term internal power-law noise with correlation
function $\zeta_{1}(t)=t^{-\lambda_{1}}/\Gamma(1-\lambda_{1})$ (i.e.,
$\hat{\zeta}_{1}(s)=s^{\lambda_{1}-1}$) for $\lambda_{1}\in(0,1)$. Thus, the
friction memory kernel due to the second FDT (23) is given by
$\gamma_{1}(t)=\alpha_{1}^{2}\zeta_{1}(t)/(k_{B}T)$. Then the GLE (6) reads
(using unit mass)

$\ddot{x}(t)+\frac{\alpha_{1}^{2}}{k_{B}T}{{}^{C}\\!}D_{t}^{\lambda_{1}}x(t)+\omega^{2}x(t)=\xi(t),\qquad\dot{x}(t)=v(t).$
(32)

Due to the power-law friction memory, the friction terms can in fact be
compactly written in terms of the fractional Caputo derivative
$({{}^{C}\\!}D_{t}^{\lambda_{1}}$, and the GLE (32) is often referred to as
the fractional Langevin equation Lutz ; pccp . The Caputo operator for
$\lambda_{1}\in(0,1)$ is defined via
${{}^{C}\\!}D_{t}^{\lambda_{1}}x(t)=\int_{0}^{t}\dot{x}(t^{\prime})(t-t^{\prime})^{-\lambda_{1}}dt^{\prime}/\Gamma(1-\lambda_{1})$.
The relaxation functions $G_{1}(t)$, $I_{1}(t)$, and $g_{1}(t)$ are obtained
by inverting the Laplace transform (11). They are given in terms of the
multinomial Prabhakar-type functions ${\cal
E}_{(\vec{\mu}),\beta}(s;\vec{a})=t^{\beta-1}E_{(\vec{\mu}),\beta}(-a_{1}t^{\mu_{1}},\ldots,-a_{N}t^{\mu}_{N})$
where $\vec{\mu}=\mu_{1},\ldots\mu_{N}$ and $\vec{a}=a_{1},\ldots,a_{N}$ (see
EBazhlekova21  and A) as

$\displaystyle
G_{1}(t)=\mathcal{L}^{-1}\left\\{(s^{2}+A_{1}s^{\lambda_{1}}+\omega^{2})^{-1};t\right\\}=\mathcal{E}_{(2,2-\lambda_{1}),2}(t;\omega^{2},A_{1}),$
(33)

and

$I_{1}(t)={\cal E}_{(2,2-\lambda_{1}),3}(t;\omega^{2},A_{1}),\quad
g_{1}(t)=\mathcal{E}_{(2,2-\lambda_{1}),1}(t;\omega^{2},A_{1}),$ (34)

where $A_{1}=\alpha_{1}^{2}/(k_{B}T)$. For $t\gg 1$ these function are
estimated by the two-parameter ML function, see (110), whose asymptotic
behaviors follow from Eq. (112), such that

$\displaystyle
G_{1}(t)\sim_{t\rightarrow\infty}A_{1}^{-1}t^{\lambda_{1}-1}E_{\lambda_{1},\lambda_{1}}(-\omega^{2}t^{\lambda_{1}}/A_{1})\sim_{t\rightarrow\infty}-\frac{A_{1}}{\omega^{4}}\frac{t^{-\lambda_{1}-1}}{\Gamma(-\lambda_{1})}=\frac{A_{1}\lambda_{1}}{\omega^{4}}\frac{t^{-\lambda_{1}-1}}{\Gamma(1-\lambda_{1})},$
(35) $\displaystyle
I_{1}(t)\sim_{t\rightarrow\infty}A_{1}^{-1}t^{\lambda_{1}}E_{\lambda_{1},\lambda_{1}+1}(-\omega^{2}t^{\lambda_{1}}/A_{1})\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\Big{[}1-\frac{A_{1}}{\omega^{2}}\frac{t^{-\lambda_{1}}}{\Gamma(1-\lambda_{1})}\Big{]},$
(36)

and

$\displaystyle
g_{1}(t)\sim_{t\rightarrow\infty}A_{1}^{-1}t^{\lambda_{1}-2}E_{\lambda_{1},\lambda_{1}-1}(-\omega^{2}t^{\lambda_{1}}/A_{1})\sim_{t\rightarrow\infty}-\frac{A_{1}}{\omega^{4}}\frac{t^{-\lambda_{1}-2}}{\Gamma(-1-\lambda_{1})}.$
(37)

The same results can be obtained by calculating the appropriate overdamped
relaxation functions such that we have $\lim_{t\to\infty}G_{1}(t)=G_{o;1}(t)$,
$\lim_{t\to\infty}I_{1}(t)=I_{o;1}(t)$, and
$\lim_{t\to\infty}g_{1}(t)=g_{o;1}(t)$.

####  II.2.2 $N$ power-law noises

The compound fractional Langevin equation, i.e., the GLE in which we use $N$
independent internal noises, for $m=1$ reads

$\ddot{x}(t)+\sum_{r=1}^{N}A_{r}{{{}^{C}\\!}D_{t}^{\lambda_{r}}}x(t)+\omega^{2}x(t)=\xi(t).$
(38)

The friction memory kernel $\gamma_{N}(t)$ becomes
$\sum_{r=1}^{N}\gamma_{r}(t)$, where $\gamma_{r}(t)$ is defined analogously to
$\gamma_{1}(t)$ above, but instead of a single internal noise we now have $N$.
Here, we assume that $0<\lambda_{1}<\dots<\lambda_{N}<1$. Using short-hand
notation we denote the relaxation functions $g_{N}(t)$, $G_{N}(t)$, and
$I_{N}(t)$ by the function $F_{N;j}(t)$ indexed by $j=1,2,3$ such that
$F_{N;1}(t)=g_{N}(t)$ is obtained for $j=1$, $F_{N;2}(t)=G_{N}(t)$ for $j=2$,
and $F_{N;3}(t)=I_{N}(t)$ for $j=3$. The auxiliary functions $F_{N;j}(t)$ read

$F_{N;j}(t)={\cal
E}_{(2,2-\lambda_{1},\ldots,2-\lambda_{N}),j}(t;\omega^{2},A_{1},\ldots,A_{N})$
(39)

in terms of the Prabakhar-type function. Following EBazhlekova21  we check the
behavior of $F_{N;j}(t)$ in (39) at long $t$. In this case we obtain that
$F_{N;j}(t)$ tends to the relaxation function for $N=1$, i.e.,

$F_{N;j}(t)\sim_{t\rightarrow\infty}A_{1}^{-1}t^{\lambda_{1}+j-3}E_{\lambda_{1},\lambda_{1}+j-2}(-\omega^{2}t^{\lambda_{1}}/A_{1}).$
(40)

Their further asymptotics in the long-time limit $t\to\infty$ depend on the
value of $j$, such that for $j=1,2,3$ we get Eqs. (37), (36), and (35),
respectively. Making use of Eq. (40) we calculate the quantities (i) to (iv)
defined in the introductory section necessary to describe the time evolution
of stochatic systems coupled to a thermal bath.

Based on these results the MSD can be shown to be given by the expression

$\langle x^{2}(t)\rangle_{N}=2k_{B}T{\cal
E}_{(2,2-\lambda_{1},\ldots,2-\lambda_{N}),3}(t;\omega^{2},A_{1},\ldots,A_{N}).$
(41)

The force-free case $\omega=0$, rewritten in terms a series of three-parameter
ML functions (see A), yields the known result (sandev_tomovski_pla, , Eqs.
(17) and (27)) such that we have

$\langle x^{2}(t)\rangle_{N}\sim_{t\rightarrow\infty}\langle
x^{2}(t)\rangle_{\mathrm{eq}}\left(1-\frac{A_{1}}{\omega^{2}}\frac{t^{-\lambda_{1}}}{\Gamma(1-\lambda_{1})}\right).$
(42)

The MSD approaches the equilibrium (thermal) value $\langle
x^{2}(t)\rangle_{\mathrm{eq}}=2k_{B}T/\omega^{2}$ in power-law fashion,
instead of the exponentially fast relaxation for the normal Orstein-Uhlenbeck
process. Note that similar power-law relaxations are known from subdiffusive
continuous time random walks Metzler1 ; Metzler3  and from the time-averaged
MSD of the fractional Langevin equation jae_pre ; jeon1  in an external
harmonic potential. From Eq. (42) we see that the noise with the smaller
exponent $\lambda_{1}$ has the dominant contribution to the oscillator
behavior in the long time limit. We also note that the same result for the MSD
in the long time limit can be obtained by Tauberian theorems (see B) if we
analyze the behavior of $\hat{I}_{N}(s)$ in the limit $s\to 0$
Gorenflo_Mainardi . The short time limit $t\to 0$ yields $I_{N}(t)\simeq
t^{2}/2-A_{2}t^{4-\lambda_{2}}/\Gamma(5-\lambda_{2})$, so that we conclude
that the noise with the larger exponent $\lambda_{2}$ has dominates the
dynamic in the short time limit.

The DACF (14) in the case of thermal initial conditions $\langle
x_{0}^{2}\rangle=k_{B}T/\omega^{2}$, $\langle x_{0}v_{0}\rangle=0$, and
$\langle\xi(t)x_{0}\rangle=0$, from Eqs. (14) and (39), we obtain

$C_{X}(t)=1-\omega^{2}{\cal
E}_{(2,2-\lambda_{1},\ldots,2-\lambda_{N}),3}(t;\omega^{2},A_{1},\ldots,A_{N}).$
(43)

Its asymptotic behavior in the long time limit via Eq. (36) reads

$\displaystyle
C_{X}(t)\sim_{t\rightarrow\infty}1-\frac{\omega^{2}}{A_{1}}t^{\lambda_{1}}E_{\lambda_{1},\lambda_{1}+1}\Big{(}-\frac{\omega^{2}}{A_{1}}t^{\lambda_{1}}\Big{)}=E_{\lambda_{1}}\Big{(}-\frac{\omega^{2}}{A_{1}}t^{\lambda_{1}}\Big{)}\sim_{t\rightarrow\infty}\frac{A_{1}}{\omega^{2}}\frac{t^{-\lambda_{1}}}{\Gamma(1-\lambda_{1})}.$
(44)

Thus we obtain a power-law decay, approaching the zero line from positive
values for $\lambda_{1}\in(0,1)$, and from negative values for
$\lambda_{1}\in(1,2)$, if we consider a memory kernel defined in Laplace space
as $\hat{\gamma}(s)=\sum_{i=1}^{N}s^{\lambda_{i}-1}$, $1<\lambda_{i}<2$. Since
$E_{\alpha}(-x)$ is a completely monotone function for $\alpha\in(0,1)$
oliveira ; KGorska ; mainardi_arxiv ; Pollard , we conclude that the
normalized displacement correlation function is completely monotone in the
long time limit for $\lambda_{1}\in(0,1)$. A more detailed analysis of
$C_{X}(t)$ is provided in subsection (II.4) for the case of high damping,
which has more practical application in the theory of anomalous dynamics in
single particle tracking and protein dynamics. We mention that the equality in
Eq. (43) was proven in KGorska18  and (KGorska23, , Remark 3).

###  II.3 Combinations of white and power-law noises

In the same way as for the power-law noises we can analyze the relaxation
functions for a mixture of $\delta$\- and power-law distributed noises. The
approach given in sandev_tomovski_pla  can be applied in the case of a
harmonic oscillator driven by $P$ white noises and $Q$ power-law noises
($P+Q=N$) as well. Here we consider the special case
$\gamma(t)=B_{1}\delta(t)+B_{2}t^{-\lambda}/\Gamma(1-\lambda)$ where
$B_{1}=\alpha^{2}/(k_{B}T)$, $B_{2}=\beta^{2}/(k_{B}T)$, and $0<\lambda<1$
tateishi ; sandev_tomovski_pla .

With the help of Eq. (104) we see that

$G(t)={\cal E}_{(2,2-\lambda,1),2}(t;\omega^{2},B_{2},B_{1}).$ (45)

From Eqs. (105) and (106) we can calculate the associated integral and
derivative which, respectively, lead to $I(t)$ and $g(t)$. Then, the MSD
obtained from Eq. (12) reads

$\langle x^{2}(t)\rangle=2k_{B}T{\cal
E}_{(2,2-\lambda,1),3}(t;\omega^{2},B_{2},B_{1}).$ (46)

Notice that Eq. (46) in the force-free case $\omega=0$, after applying Eq.
(117), can be recognized as (sandev_tomovski_pla, , Eq. (27) for
$\lambda_{1}=\lambda$, $\lambda_{2}=1$), namely,

$\langle
x^{2}(t)\rangle=2k_{B}T\sum_{n=0}^{\infty}\left(-B_{2}\right)^{n}t^{(2-\lambda)n+2}E_{1,(2-\lambda)n+3}^{n+1}(-B_{1}t).$
(47)

The asymptotic of the MSD at short and long times are calculated from Eqs.
(107) and (108), yielding

$\langle x^{2}(t)\rangle\sim_{t\rightarrow
0}\frac{t^{2}}{2}-\frac{t^{3}}{3!}-\frac{t^{4}}{4!}-\frac{t^{4-\lambda}}{\Gamma(5-\lambda)}$
(48)

for $t\to 0$, and for $t\gg 1$ we have

$\displaystyle\langle
x^{2}(t)\rangle\sim_{t\rightarrow\infty}\frac{2k_{B}T}{B_{2}}t^{\lambda}E_{\lambda,1+\lambda}\Big{(}-\frac{\omega^{2}}{B_{2}}t^{\lambda}\Big{)}\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\left[1-\frac{B_{2}}{\omega^{2}}\frac{t^{-\lambda}}{\Gamma(1-\lambda)}\right],$
(49)

Which means that the MSD has a power-law decay to the equilibrium value
$\langle x^{2}(t)\rangle_{\mathrm{eq}}=2k_{B}T/\omega^{2}$. We conclude that
the power-law noise is dominant in the long time limit, and the white noise in
the short time limit, as naively expected.

For the DACF $C_{X}(t)$, we obtain from relation (14) that

$C_{X}(t)=1-\omega^{2}{\cal E}_{(2,2-\lambda,1),3}(t;\omega^{2},B_{2},B_{1}),$
(50)

from where the long time limit follows,

$\displaystyle
C_{X}(t)\sim_{t\rightarrow\infty}1-\frac{\omega^{2}}{B_{2}}t^{\lambda}E_{\lambda,1+\lambda}\Big{(}-\frac{\omega^{2}}{B_{2}}t^{\lambda}\Big{)}\sim_{t\rightarrow\infty}\frac{B_{2}}{\omega^{2}}\frac{t^{-\lambda}}{\Gamma(1-\lambda)}.$
(51)

We conclude that $C_{X}(t)$ in the long time limit is a completely monotone
function since $0<\lambda<1$, showing a power-law decay to zero.

Following the same procedure one may consider mixtures of white noises, power-
law noises, and ML type noises. The calculation of the relaxation function can
be represented in terms of multinomial Prabhakar functions sandev_tomovski_pla
; EBazhlekova21 .

###  II.4 Overdamped limit

Next, we analyze the high-damping limit, in which the inertial term
$m\ddot{x}(t)$ can be neglected. Then, the relaxation functions are defined by
Eqs. (18).

The case of a mixture of $N$ internal white noises considered in II.1, for the
MSD $\langle x^{2}(t)\rangle_{o}=2k_{B}TI_{o}(t)$ and the DACF $C_{X,o}(t)$
yields

$\displaystyle\langle x^{2}(t)\rangle_{o}$
$\displaystyle=2k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\kappa
s+\omega^{2}};t\right\\}=\frac{2k_{B}T}{\omega^{2}}\left[1-\exp\left(-\frac{\omega^{2}}{\kappa}t\right)\right],$
(52) $\displaystyle C_{X,o}(t)$
$\displaystyle=\exp\Big{(}-\frac{\omega^{2}}{\kappa}t\Big{)}=\exp\left(-\frac{\mathcal{D}\omega^{2}}{k_{B}T}t\right),$
(53)

where $\kappa$ is given below Eq. (25) and ${\cal D}$ is defined in Eq. (26).
We find that $C_{X,o}(t)$ has a monotonic exponential decay, as expected.

For the case of $N$ power-law noises, it follows from relation (12) in the
overdamped case that

$\displaystyle\langle x^{2}(t)\rangle_{o}$
$\displaystyle=\frac{2k_{B}T}{A_{N}}{\cal
E}_{(\lambda_{N},\lambda_{N}-\lambda_{1},\ldots,\lambda_{N}-\lambda_{N-1});\lambda_{N}+1}\Big{(}t,\frac{\omega^{2}}{A_{N}},\frac{A_{1}}{A_{N}},\ldots,\frac{A_{N-1}}{A_{N}}\Big{)},$
(54) $\displaystyle C_{X,o}(t)$ $\displaystyle=1-\frac{\omega^{2}}{A_{N}}{\cal
E}_{(\lambda_{N},\lambda_{N}-\lambda_{1},\ldots,\lambda_{N}-\lambda_{N-1});\lambda_{N}+1}\Big{(}t,\frac{\omega^{2}}{A_{N}},\frac{A_{1}}{A_{N}},\ldots,\frac{A_{N-1}}{A_{N}}\Big{)}.$
(55)

The DACF shows an asymptotic power-law decay in the long time limit, i.e.,

$C_{X,o}(t)\sim_{t\rightarrow\infty}E_{\lambda_{1}}\left(-\frac{\omega^{2}}{A_{1}}t^{\lambda_{1}}\right)\sim_{t\rightarrow\infty}\frac{A_{1}}{\omega^{2}}\frac{t^{-\lambda_{1}}}{\Gamma\left(1-\lambda_{1}\right)}.$
(56)

Thus, we can use equation (16) in the overdamped limit case, which is simpler
instead of the GLE (6), to analyze the asymptotic behavior of the harmonic
oscillator in the long time limit.

A graphical representation of the DACF is presented in figure 2. From panel
(a) we see that by changing the values of the parameters $\lambda_{1}$ and
$\lambda_{2}$ for fixed frequency $\omega$ there appears a non-monotonic decay
of $C_{X,o}(t)$ without zero crossings, approaching the zero line at infinity
(see solid line), a monotonic decay without zero crossings to zero at infinity
(dashed line), a non-monotonic decay without zero crossings approaching zero
at a finite time instant and at infinity (dot-dashed line), or a non-monotonic
decay with zero crossings approaching zero at infinity (dotted line). All long
time decays of the DACF are of power-law form to zero, as it can be
anticipated from relation (56). The behavior of $C_{X}(t)$ for different
values of the frequency $\omega$ and fixed values of $\lambda_{1}$ and
$\lambda_{2}$ is shown in figure 2(b) and (c). We see that there are different
critical frequencies, i.e., the frequency at which $C_{X,o}(t)$ changes its
behavior, for instance, from non-monotonic to monotonic decay without zero
crossings, or the frequency at which $C_{X,o}(t)$ crosses the zero line. Such
different types of critical frequencies were discussed in Burov_Barkai .

Figure 3 depicts the MSD for unconfined and confined motion. For free motion
the MSD may have monotonic, non-monotonic, and oscillatory behavior, turning
into a power-law form in the long time limit. For the confined case the MSD
has different behaviors at intermediate times, and in the long time limit the
MSD has a power-law approach to the equilibrium value $\langle
x^{2}(t)\rangle_{\mathrm{eq}}=2k_{B}T/\omega^{2}$.

(a)  (b)

(c)

Figure 2: Normalized displacement correlation function (55) for the following
cases for $N=2$: (a) $\omega=1$, $\lambda_{2}=3/2$, $\lambda_{1}=1/8$, (solid
line), $\lambda_{1}=1/2$ (dashed line), $\lambda_{1}=7/8$ (dot-dashed line),
$\lambda_{1}=5/4$ (dotted line); (b) $\lambda_{1}=1/8$, $\lambda_{2}=3/4$,
$\omega=0.5$ (solid line), $\omega=0.75$ (dashed line), $\omega=1.8959706$
(dot-dashed line), $\omega=2.5$ (dotted line); (c) $\lambda_{1}=3/4$,
$\lambda_{2}=3/2$, $\omega=1$ (solid line), $\omega=1.617$ (dashed line),
$\omega=3$ (dot-dashed line).

(a)  (b)

(c)

Figure 3: MSD (54) for the following cases for $N=2$: (a) $\omega=0$,
$\lambda_{1}=1/2$, $\lambda_{2}=1$, (solid line), $\lambda_{2}=5/4$ (dashed
line), $\lambda_{2}=3/2$ (dot-dashed line), $\lambda_{2}=7/4$ (dotted line);
(b) $\omega=0$, $\lambda_{1}=0.1$, $\lambda_{2}=1$, (solid line),
$\lambda_{2}=5/4$ (dashed line), $\lambda_{2}=3/2$ (dot-dashed line),
$\lambda_{2}=7/4$ (dotted line); (c) $\lambda_{1}=1/2$, $\lambda_{2}=7/4$,
$\omega=0$ (upper solid line), $\omega=1$ (dashed line), $\omega=1.5$ (dot-
dashed line), $\omega=2$ (dotted line), $\omega=2.5$ (lower solid line).

In the same way as above we obtain the following DACF for the case of a
mixture of a white noise a and power-law noise (see results II.3),

$\displaystyle
C_{X,o}(t)=\sum_{n=0}^{\infty}\left(-\frac{\omega^{2}}{B_{1}}\right)^{n}t^{n}E_{1-\lambda,n+1}^{n}\left(-\frac{B_{2}}{B_{1}}t^{1-\lambda}\right),$
(57)

for $0<\lambda<1$, and

$\displaystyle
C_{X,o}(t)=\sum_{n=0}^{\infty}\left(-\frac{\omega^{2}}{B_{2}}\right)^{n}t^{\lambda
n}E_{\lambda-1,\lambda
n+1}^{n}\left(-\frac{B_{1}}{B_{2}}t^{\lambda-1}\right),$ (58)

for $1<\lambda<2$, where $B_{1}$ and $B_{2}$ are given at the beginning of
Sec. II.3. These results are equivalent to those obtained in subsection II.3
in the long time limit for the GLE when the inertial term is not neglected, as
it should be.

##  III Distributed-order Langevin equations

It has been shown that the distributed-order differential equations are
suitable tool for modeling ultraslow relaxation and diffusion processes
chechkin1 ; chechkin2 ; Eab2 ; kochubei ; mainardi_book ; Mainardi_distributed
; sandev_tomovski_pla . In case of distributed order differential equations
one uses the following memory kernel (see for example Mainardi_distributed )

$\gamma(t)=(k_{B}T)^{-1}\,\int_{0}^{1}p(\lambda)\frac{t^{-\lambda}}{\Gamma(1-\lambda)}d\lambda,$
(59)

where $p(\lambda)$ is a dimensionless, non-negative weight function with
$\int_{0}^{1}p(\lambda)\mathrm{d}\lambda=c$, where $c$ is a constant. When
$c=1$, $p(\lambda)$ is normalized. If we substitute the distributed-order
memory kernel in the GLE (6) we transform it to the following distributed-
order Langevin equation

$\displaystyle\ddot{x}(t)+\frac{1}{k_{B}T}\int_{0}^{1}p(\lambda)^{C\\!}D^{\lambda}_{t}x(t)d\lambda+\frac{dV(x)}{dx}=\xi(t),\quad\dot{x}(t)=v(t).$
(60)

Here we note that assumption (8) is satisfied for distributed-order Langevin
equations since $\lim_{s\rightarrow 0}s\hat{\gamma}(s)\simeq\lim_{s\rightarrow
0}\int_{0}^{1}p(\lambda)s^{\lambda}d\lambda=0$. Thus we can use the
representations of MSD, VACF, time-dependent diffusion coefficient and DACFs
in terms of the relaxation functions.

We now consider the distributed-order Langevin equation (60) in the presence
of the constant external force $F$. Then we have

$\displaystyle\ddot{x}(t)+\frac{1}{k_{B}T}\int_{0}^{1}p(\lambda)^{C\\!}D^{\lambda}_{t}x(t)\,d\lambda-F=\xi(t),\quad\dot{x}(t)=v(t),$
(61)

from which, by the Laplace transform method, we obtain

$x(t)=\langle
x(t)\rangle_{F}+\int_{0}^{t}G(t-t^{\prime})\xi(t^{\prime})dt^{\prime},$ (62)

where $\langle x(t)\rangle_{F}=FI(t)$ and $\hat{\gamma}(s)$ is given by Eq.
(59), and where

$G(t)=\mathcal{L}^{-1}\left\\{\frac{1}{s^{2}+s\hat{\gamma}(s)};t\right\\}\quad\text{and}\quad
I(t)=\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{s^{2}+s\hat{\gamma}(s)};t\right\\}.$
(63)

The latter expression corresponds to (11) for a free particle ($\omega=0$).
Thus, we conclude that the generalized Einstein relation $\langle
x^{2}(t)\rangle_{F}=[F/(2k_{B}T)]\langle x^{2}(t)\rangle_{F=0}$ is satisfied
for the distributed-order Langevin equation (60), where $\langle
x^{2}(t)\rangle_{F=0}=2k_{B}T\mathcal{L}^{-1}\\{s^{-1}[s^{2}+s\hat{\gamma}(s)]^{-1};t\\}$
is the MSD for the case of a free particle.

###  III.1 Force-free case

We first consider the distributed-order Langevin equation for a free particle
with $\omega=0$. Note that a weight function of the form
$p(\lambda)=\sum_{i=1}^{N}\alpha_{i}^{2}\delta(\lambda-\lambda_{i})$ yields
the fractional Langevin equation and we have the compound white Gaussian
noise, see Subsect. II.2.2. For $p(\lambda)=\alpha^{2}$ we obtain the
uniformly distributed noise kochubei

$k_{B}T\gamma(t)=\alpha^{2}\int_{0}^{1}\frac{t^{-\lambda}}{\Gamma\left(1-\lambda\right)}d\lambda,$
(64)

which was used by Kochubei in the theory of evolution equations. In
sandev_tomovski_pla  we analyzed the GLE for a free particle with friction
memory kernel of distributed order (64) and showed that the MSD in the long
time limit is given by

$\langle
x^{2}(t)\rangle\sim_{t\to\infty}\frac{2\left(k_{B}T\right)^{2}}{\alpha^{2}}\Big{[}C+\log
t+{\rm e}^{t}\mathrm{E}_{1}(t)\Big{]},$ (65)

where $C=0.577216$ is the Euler-Mascheroni (or Euler’s) constant,
$\mathrm{Ei}(-t)=-\int_{t}^{\infty}(e^{-x}/x)dx$ is the exponential integral
erdelyi , and $\mathrm{E}_{1}(t)=-\mathrm{Ei}(-t)$. From the asymptotic
expansion ${\rm E}_{1}(t)\sim_{t\to\infty}t^{-1}{\rm
e}^{-t}\sum_{k=0}^{n-1}(-1)^{k}k!t^{-k}$ mainardi_book , which has an error of
order ${\cal O}(n!t^{-n})$, we obtained that the particle shows ultraslow
diffusion, i.e., $\langle
x^{2}(t)\rangle\sim_{t\to\infty}[2(k_{B}T)^{2}/\alpha^{2}](C+\log t)$
sandev_tomovski_pla .

Consider now the power-law case $p(\lambda)=\alpha^{2}\beta\lambda^{\beta-1}$,
where $\beta>0$. Then the friction memory kernel becomes

$k_{B}T\gamma(t)=\alpha^{2}\int_{0}^{1}\beta\lambda^{\beta-1}\frac{t^{-\lambda}}{\Gamma(1-\lambda)}d\lambda.$
(66)

For the MSD in the long time limit, we obtain (see (63) based on Tauberian
theorems)

$\langle
x^{2}(t)\rangle\sim_{t\to\infty}\frac{2(k_{B}T)^{2}}{\alpha^{2}}\frac{\log^{\beta}t}{\Gamma(1+\beta)}.$
(67)

From this result we conclude that this power-law model leads to ultraslow
diffusion (for $\beta=1$) or a strong anomaly (in Ref. drager  a strong
anomaly means the behavior of form $\langle x^{2}(t)\rangle\sim\log^{\nu}t$,
which for $\nu=4$ has the same form as the MSD of the Sinai diffusion model in
a random force field aljaz ; sinai , compare also to the discussion in igorg
). This result for the MSD is equivalent to the one obtained by Eab and Lim
Eab2 , since the long time limit considered here corresponds to the case of
the overdamped limit (the inertial term $\ddot{x}(t)$ is neglected) studied in
Eab2 . The same result for the MSD can be obtained from the distributed-order
diffusion equation chechkin2  in which the weight function corresponds to the
one considered in the memory kernel (66).

Next we consider the weight function
$p(\lambda)=\alpha^{2}/(\lambda_{2}-\lambda_{1})$ for
$0\leq\lambda_{1}<\lambda<\lambda_{2}\leq 1$, and $p(\lambda)=0$ otherwise
chechkin1 . With the help of relation (63) the MSD yields in the form

$\displaystyle\langle
x^{2}(t)\rangle=2k_{B}T\left[\frac{t^{2}}{2}+\sum_{n=1}^{\infty}\Big{(}-\frac{\alpha^{2}}{k_{B}T}\Big{)}^{n}\frac{1}{(\lambda_{2}-\lambda_{1})^{n}}\sum_{k=0}^{n}\binom{n}{k}(-1)^{k}\mu(t,n-1,(2-\lambda_{2})n+(\lambda_{2}-\lambda_{1})k+2)\right],$
(68)

which in the long time limit becomes (use the Tauberian theorem in Eq. (63))

$\displaystyle\langle
x^{2}(t)\rangle\sim_{t\rightarrow\infty}\frac{2\left(k_{B}T\right)^{2}}{\alpha^{2}}\sum_{n=1}^{\infty}\frac{t^{\lambda_{2}-(\lambda_{2}-\lambda_{1})n}}{\Gamma(\lambda_{2}-(\lambda_{2}-\lambda_{1})n+1)}\left[\ln
t^{\lambda_{2}-\lambda_{1}}-(\lambda_{2}-\lambda_{1})\psi(\lambda_{2}-(\lambda_{2}-\lambda_{1})n+1)\right],$
(69)

where we used the Volterra function $\mu\left(t,\beta,\alpha\right)$ defined
in C, and $\psi=\frac{\Gamma^{\prime}}{\Gamma}$ is the digamma function. Thus,
we conclude that in the long time limit the MSD has the ultraslow form

$\langle
x^{2}(t)\rangle\sim_{t\rightarrow\infty}\frac{2(k_{B}T)^{2}}{\alpha^{2}}\frac{\lambda_{2}-\lambda_{1}}{\Gamma(1+\lambda_{1})}t^{\lambda_{1}}\ln
t.$ (70)

For $\lambda_{1}=0$ and $\lambda_{2}=1$ we arrive at result (64) obtained for
uniformly distributed noise with a logarithmic MSD dependence on time.

Such relaxation patterns of logarithmic form and ultraslow diffusion have been
observed in the analysis of distributed-order relaxation and diffusion
equations by Kochubei kochubei2 ; kochubei3 ; kochubei4 ; kochubei  as well as
by Mainardi et al. mainardi_book ; Mainardi_distributed .

###  III.2 Harmonic oscillator: overdamped limit

We now move on to consider the distributed-order GLE for an harmonic
oscillator with different weight functions in the overdamped limit. For the
uniformly distributed noise (64) we obtain the following form of the MSD

$\displaystyle\langle x^{2}(t)\rangle_{o}$
$\displaystyle=2k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\frac{\alpha^{2}}{k_{B}T}\frac{s-1}{\log
s}+\omega^{2}};t\right\\}=2k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}\log
s}{\omega^{2}\log s+As-A};t\right\\},$ (71)

where $A=\alpha^{2}/(k_{B}T)$. From Tauberian theorems (see B) we analyze the
asymptotic behavior in the long and short time limits. At long times
($t\rightarrow\infty$, equivalent to $s\rightarrow 0$) it follows that

$\displaystyle\langle
x^{2}(t)\rangle_{o}\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{1-\frac{A}{\omega^{2}\log
s}};t\right\\}\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{1}{s}+\frac{\frac{A}{\omega^{2}}}{s\log
s};t\right\\}=\frac{2k_{B}T}{\omega^{2}}\left[1+\frac{A}{\omega^{2}}\nu(t)\right],$
(72)

where $\nu(t)$ is the Volterra function, see C. Thus, the DACF in the long
time limit becomes

$\displaystyle
C_{X,o}(t)\sim_{t\rightarrow\infty}-\frac{A}{\omega^{2}}\nu(t).$ (73)

At short times ($t\to 0$, or $s\to\infty$), we find for the MSD and DACF that

$\langle x^{2}(t)\rangle_{o}\sim_{t\rightarrow
0}\frac{2k_{B}T}{A}\mathcal{L}^{-1}\left\\{\frac{\log
s}{s^{2}};t\right\\}=\frac{2k_{B}T}{A}t\left(\log\frac{1}{t}+1-\gamma\right),\quad$
(74)

and

$C_{X,o}(t)\sim_{t\rightarrow
0}1-\frac{\omega^{2}}{A}t\left(\log\frac{1}{t}+1-\gamma\right),$ (75)

respectively.

Now consider the distributed-order noise (66). From the relaxation function
$I_{o}(t)$ the MSD becomes

$\langle
x^{2}(t)\rangle_{o}=\frac{2k_{B}T}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\frac{A}{\omega^{2}}\nu(\log\frac{1}{s})^{-\nu}\gamma(\nu,-\log
s)+1};t\right\\},$ (76)

where $\gamma(a,\sigma)=\int_{0}^{\sigma}t^{a-1}\exp(-t)dt$ is the lower
incomplete Gamma function. Taking the long time limit yields

$\displaystyle\langle
x^{2}(t)\rangle_{o}\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\left\\{1-\frac{A}{\omega^{2}}\Gamma(1+\nu)\mathcal{L}^{-1}\left\\{\frac{1}{s}{\left(\log\frac{1}{s}\right)^{-\nu}};t\right]\right\\}\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\left[1-\frac{A}{\omega^{2}}\frac{\Gamma(1+\nu)}{\log^{\nu}t};t\right],$
(77)

from which we obtain the asymptotic behavior of the DACF,

$C_{X,o}(t)\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\frac{A}{\omega^{2}}\frac{\Gamma(1+\nu)}{\log^{\nu}t}.$
(78)

In the short time limit we find

$\langle x^{2}(t)\rangle_{o}\sim_{t\rightarrow
0}\frac{2k_{B}T}{A}\frac{\log^{\nu}\frac{1}{t}}{\Gamma(1+\nu)}$ (79)

and

$C_{X,o}(t)\sim_{t\rightarrow
0}1-\frac{\omega^{2}}{A}\frac{\log^{\nu}\frac{1}{t}}{\Gamma(1+\nu)}.$ (80)

In a similar way we obtain for the distributed-order noise (59) with weight
function $p(\lambda)=\alpha^{2}/(\lambda_{2}-\lambda_{1})$
($0\leq\lambda_{1}<\lambda<\lambda_{2}\leq 1$, and $p(\lambda)=0$ otherwise)
the asymptotic forms of the MSD:

$\langle
x^{2}(t)\rangle_{o}\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}\left\\{1-\frac{A/\omega^{2}}{\left(\lambda_{2}-\lambda_{1}\right)}\big{[}\nu(t,-\lambda_{2})-\nu(t,-\lambda_{1})\big{]}\right\\},$
(81)

in the long time limit, and

$\displaystyle\langle x^{2}(t)\rangle_{o}\sim_{t\rightarrow
0}\frac{2k_{B}T}{A}\sum_{n=0}^{\infty}\frac{t^{\left(\lambda_{2}-\lambda_{1}\right)n+\lambda_{2}}}{\Gamma[(\lambda_{2}-\lambda_{1})n+\lambda_{2}+1]}\left[\log
t^{-\left(\lambda_{2}-\lambda_{1}\right)}-\left(\lambda_{2}-\lambda_{1}\right)\psi\left(\left(\lambda_{2}-\lambda_{1}\right)n+\lambda_{2}+1\right)\right]$
(82)

in the short time limit. From here we easily can find the DACF.

From these results we conclude that the distributed-order GLE may be used to
model various anomalous diffusive behaviors, such as ultraslow diffusion,
strong anomaly, and other complex diffusive regimes.

##  IV Further generalizations

###  IV.1 LE with power-logarithmic distributed order noises

Let us consider the Langevin equation (6) with the logarithmically
distributed-order friction kernel

$\gamma(t)=(k_{B}T)^{-1}\int_{0}^{1}\Gamma(3/2-\lambda)\frac{\log^{\lambda-1}t}{\sqrt{t}}d\lambda\quad\mbox{with
}\left(\hat{\gamma}(s)=\frac{\pi}{k_{B}T}\frac{s-1}{\sqrt{s}\log s}\right),$
(83)

which satisfies the condition (8). The relaxation function $I(t)$ defined by
Eq. (11) assumes the form

$I(t)=\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{s^{2}+\frac{\pi}{k_{B}T}\frac{s(s-1)}{\sqrt{s}\log{s}}+\omega^{2}};t\right\\}.$
(84)

As it is challenging to calculate the exact form of Eq. (84) for general $t>0$
we concentrate on the asymptotic behavior for small $s\ll 1$, yielding

$\displaystyle\hat{I}(s)$ $\displaystyle\sim_{s\rightarrow
0}\frac{s^{-1}}{\omega^{2}}\frac{1}{1-\frac{\pi}{k_{B}T\omega^{2}}\frac{s(1-s)}{\sqrt{s}\log
s}}\simeq\frac{1}{\omega^{2}}\left[1+\frac{\pi}{k_{B}T\omega^{2}}\frac{1}{\sqrt{s}\log
s}-\frac{\pi}{k_{B}T\omega^{2}}\frac{\sqrt{s}}{\log s}\right],$

where we take the zeroth and first-order terms of the series expansion. If $s$
tends zero then the term $\sqrt{s}/\log s$ can be neglected. Thus, we have

$\displaystyle\hat{I}(s)\sim_{s\rightarrow
0}\frac{1}{\omega^{2}}\left[1+\frac{\pi}{k_{B}T\omega^{2}}\frac{1}{\sqrt{s}\log
s}\right]\quad\text{ and }\quad
I(t)\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\left[1+\frac{\pi}{k_{B}T\omega^{2}}\nu(t,-1/2)\right].$
(85)

This allows us to find the MSD and DACF for $t\rightarrow\infty$,

$\langle
x^{2}(t)\rangle\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}+\frac{2\pi}{\omega^{4}}\nu(t,-1/2)\quad\text{
and }\quad
C_{X}(t)\sim_{t\rightarrow\infty}-\frac{\pi}{k_{B}T\omega^{2}}\nu(t,-1/2).$

Note that for $\omega=0$ we get

$\displaystyle\langle x^{2}(t)\rangle$
$\displaystyle=2k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{s^{2}+\frac{\pi}{k_{B}T}\frac{s(s-1)}{\sqrt{s}\log
s}};t\right\\}\sim_{t\rightarrow\infty}2k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\frac{\pi}{k_{B}T}\frac{s(s-1)}{\sqrt{s}\log
s}};t\right\\}$
$\displaystyle\sim_{t\rightarrow\infty}\frac{2(k_{B}T)^{2}}{\pi}\mathcal{L}^{-1}\left\\{\frac{\log
s}{s^{3/2}(s-1)};t\right\\}\sim_{t\rightarrow\infty}\frac{2(k_{B}T)^{2}}{\pi}\mathcal{L}^{-1}\left\\{\frac{\log\frac{1}{s}}{s^{3/2}};t\right\\}$
$\displaystyle\sim_{t\rightarrow\infty}\frac{2(k_{B}T)^{2}}{\pi}\frac{1}{\Gamma(3/2)}\sqrt{t}\log
t=\frac{4(k_{B}T)^{2}}{\pi^{3/2}}\sqrt{t}\log t.$ (86)

###  IV.2 Langevin equation with distributed-order Mittag-Leffler noise

In the next example, we study the GLE for an harmonic oscillator with the ML
friction memory kernel

$\gamma(t)=\frac{1}{k_{B}T}\int_{0}^{1}E_{\lambda}\left(-t^{\lambda}\right)d\lambda\quad\left(\hat{\gamma}(s)=\frac{1}{k_{B}T}\frac{\log\frac{s+1}{2}}{s\log
s}\right),$ (87)

which satisfies the condition (8). The long time asymptotic of the relaxation
function $I(t)$ involved in the MSD and the DACF can be calculated by use of
the Tauberian theorem ( B) in which we consider the limit $s\to 0$ of
$\hat{I}(s)$. Thus, we have

$\displaystyle
I(t)=\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{s^{2}+\frac{1}{k_{B}T}\frac{\log\frac{s+1}{2}}{\log
s}+\omega^{2}};t\right\\}\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{1}{s}-\frac{1}{k_{B}T\omega^{2}}\frac{\log
2}{s\log s};t\right\\}=\frac{1}{\omega^{2}}\left[1-\frac{\log
2}{k_{B}T\omega^{2}}\nu(t)\right].$ (88)

In the force-free case with $\omega=0$ we observe ultraslow diffusion, since
the relaxation function $I(t)$ has a logarithmic time dependence,

$I(t)\sim_{t\rightarrow\infty}\frac{k_{B}T}{\log
2}\mathcal{L}^{-1}\left\\{\frac{\log\frac{1}{s}}{s};t\right\\}=\frac{{k_{B}T}}{\log
2}\log t=k_{B}T\log_{2}t.$ (89)

###  IV.3 Distributed-order Langevin equation with Caputo-Prabhakar derivative

In this part we now consider the GLE (6) for an harmonic oscillator with the
choice of the distributed-order Prabhakar friction memory kernel

$\gamma(t)=\frac{1}{k_{B}T}\int_{0}^{1}t^{-\lambda}E_{{\rho},1-\lambda}^{\delta}(-t^{{\rho}})d\lambda\quad\text{
with }\quad\hat{\gamma}(s)=\frac{1}{k_{B}T}\frac{s-1}{s\log
s}(1+{s^{-\rho}})^{-\delta}.$ (90)

We consider two cases, the first one for $0<\rho,\delta<1$ and the second one
for $-1<\rho,\delta<0$. In both cases the condition (8) is satisfied. Using
the Tauberian theorem we calculate the asymptotics of the relaxation function
$I(t)$.

(i) In the case $0<\rho,\delta<1$ and for $t\to\infty$ we have

$\displaystyle I(t)$
$\displaystyle\sim_{t\rightarrow\infty}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\omega^{2}+\frac{1}{k_{B}T}\frac{s-1}{\log
s}(1+s^{-\rho})^{-\delta}};t\right\\}\sim_{t\rightarrow\infty}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\omega^{2}+\frac{1}{k_{B}T}\frac{s-1}{\log
s}s^{\rho\delta}};t\right\\}$
$\displaystyle\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{1+\frac{1}{k_{B}T\omega^{2}}\frac{s^{\rho\delta}}{\log\frac{1}{s}}};t\right\\}\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{1}{s}-\frac{1}{k_{B}T\omega^{2}}\frac{s^{-1+\rho\delta}}{\log\frac{1}{s}};t\right\\}$
$\displaystyle\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\left[1-\frac{1}{k_{B}T\omega^{2}\Gamma(1-\rho\delta)}\frac{t^{-\rho\delta}}{\log{t}}\right].$
(91)

From the asymptotic form of $I(t)$ we obtain the long time behaviors of the
MSD and DACF,

$\langle
x^{2}(t)\rangle\sim_{t\rightarrow\infty}\frac{2k_{B}T}{\omega^{2}}-\frac{1}{\omega^{4}\Gamma(1-\rho\delta)}\frac{t^{-\rho\delta}}{\log
t}$ (92)

and

$C_{X}(t)\sim_{t\rightarrow\infty}\frac{1}{k_{B}T\omega^{2}\Gamma(1-\rho\delta)}\frac{t^{-\rho\delta}}{\log{t}}.$
(93)

In addition, note that in the free particle case ($\omega=0$) the long time
limit (or for the overdamped case when we neglect the term with $s^{2}$)
becomes

$\displaystyle I(t)$
$\displaystyle\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}(1+{s^{-\rho}})^{\delta}\log
s}{{s-1}};t\right\\}\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{\frac{\log
s}{s^{1+\rho\delta}(s-1)};t\right\\}$
$\displaystyle\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{s^{-1-\rho\delta}\log\frac{1}{s};t\right\\}\sim_{t\rightarrow\infty}k_{B}T\frac{t^{\rho\delta}}{\Gamma(1+\rho\delta)}\log{t}.$
(94)

(ii) Replacing $\rho$ by $-\rho$ and $\delta$ by $-\delta$ in Eq. (90) we find

$\gamma(t)=\frac{1}{k_{B}T}\int_{0}^{1}t^{-\lambda}E_{{-\rho},1-\lambda}^{-\delta}(-t^{{-\rho}})d\lambda\quad\text{and}\quad\left(\hat{\gamma}(s)=\frac{1}{k_{B}T}\frac{s-1}{s\log
s}(1+{s^{\rho}})^{\delta}\right),$ (95)

where $0<\rho,\delta<1$. The relaxation function $I(t)$ for $\omega\neq 0$
reads

$\displaystyle I(t)$
$\displaystyle\sim_{t\rightarrow\infty}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\frac{1}{k_{B}T}\frac{s-1}{\log
s}(1+{s^{\rho}})^{\delta}+\omega^{2}};t\right\\}\sim_{t\rightarrow\infty}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\frac{1}{k_{B}T}\frac{s-1}{\log
s}+\omega^{2}};t\right\\}$
$\displaystyle=\frac{1}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{1-\frac{1}{k_{B}T\omega^{2}}\frac{1-s}{\log{s}}};t\right\\}\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{1}{s}+\frac{1}{k_{B}T\omega^{2}}\frac{1}{s\log
s};t\right\\}$
$\displaystyle=\frac{1}{\omega^{2}}\left[1+\frac{\nu(t)}{k_{B}T\omega^{2}}\right].$
(96)

For $\omega=0$ the $t\to\infty$-asymptotic of $I(t)$ is equal to

$\displaystyle I(t)$
$\displaystyle=k_{B}T\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{\frac{s-1}{\log
s}(1+{s^{\rho}})^{\delta}};t\right\\}\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{\frac{\log
s}{s(s-1)};t\right\\}$ $\displaystyle=k_{B}T\big{[}C+\log
t+e^{t}E_{1}(t)\big{]}.$ (97)

Note that the relevant MSD is the same as obtained for the force-free case by
using the distributed-order fractional derivative, i.e., Eq. (65) for
$\alpha=1$.

###  IV.4 Langevin equation with distributed-order Volterra function

We consider the GLE for which the friction kernel is based on the Volterra
functions $\nu(t,-\alpha)$ and $\mu(t,-\beta)$ for $0<\alpha,\beta<1$, see C.
In the first case the friction memory kernel in $t$ and $s$ reads

$\gamma_{1}(t)=\frac{1}{k_{B}T}\int_{0}^{1}\nu(t,-\alpha)d\alpha\quad\text{and}\quad\hat{\gamma}_{1}(s)=\frac{1}{k_{B}T}\frac{s-1}{s\log^{2}s},$
(98)

thus satisfying condition (8). For the GLE for the stochastic harmonic
oscillator ($\omega\neq 0$) the long time limit ($t\to\infty$) of the
relaxation function $I_{1}(t)$ yields in the form

$\displaystyle I_{1}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{s^{2}+\frac{1}{k_{B}T}\frac{s-1}{\log^{2}s}+\omega^{2}};t\right\\}\sim_{t\rightarrow\infty}\frac{1}{\omega^{2}}\mathcal{L}^{-1}\left\\{\frac{1}{s}+\frac{1}{k_{B}T\omega^{2}}\frac{1}{s\log^{2}s};t\right\\}$
$\displaystyle=\frac{1}{\omega^{2}}\left[1+\frac{\mu(t,1)}{k_{B}T\omega^{2}}\right].$
(99)

In the force-free case ($\omega=0$) we obtain (KGorska, , Eq.(45))

$\displaystyle I_{1}(t)$
$\displaystyle=\mathcal{L}^{-1}\left\\{\frac{s^{-1}}{s^{2}+\frac{1}{k_{B}T}\frac{s-1}{\log^{2}s}};t\right\\}\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{\frac{\log^{2}s}{s(s-1)};t\right\\}$
$\displaystyle=k_{B}T\Big{[}-(C^{2}+\pi^{2}/6)e^{t}+2te^{t}\Phi_{1;1}^{\star,(1,1)}(-t,3,1)-2(C+\log
t)e^{t}E_{1}(t)$ $\displaystyle-(2C+\log t)(e^{t}+1)(\log
t)-C^{2}+\pi^{2}/6\Big{]}.$ (100)

Here $\Phi_{1;1}^{\star,(1,1)}$ is the Hurwitz-Lerch function, see LiSr

In the second case, the friction in $t$ and $s$ spaces turns out to be

$\gamma_{2}(t)=\frac{1}{k_{B}T}\int_{0}^{1}\mu(t,-\beta)d\beta\quad\text{and}\quad\hat{\gamma}_{2}(s)=\frac{1}{k_{B}T}\frac{(\log
s)-1}{s(\log s)(\log\log s)}.$ (101)

Here we consider only the force-free case, i.e., $\omega=0$. In that case we
find the asumptotic of $I_{2}(t)$ at long times. Using that $\log\log s$ is a
slowly varying function, by Tauberian theorems we find that

$\displaystyle I_{2}(t)$
$\displaystyle\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{\frac{(\log
s)(\log\log s)}{s(\log
s-1)};t\right\\}=k_{B}T\mathcal{L}^{-1}\left\\{\frac{\log\log
s}{s(1-\frac{1}{\log s})};t\right\\}$
$\displaystyle\sim_{t\rightarrow\infty}k_{B}T\mathcal{L}^{-1}\left\\{\frac{1}{s}{\log\log\frac{1}{s^{-1}}};t\right\\}=k_{B}T\log\log\frac{1}{t}.$
(102)

##  V Summary

We studied the GLE for an external harmonic potential and the free particle
limit for various cases of the friction memory kernel. In particular we
derived the associated MSD and DACF for the different chosen forms for the
friction memory kernel of Dirac delta, power-law, and distributed-order forms.
Various anomalous diffusive behaviors, such as subdiffusion, superdiffusion,
ultraslow diffusion, and strong anomaly, are observed. Special attention was
paid on distributed-order GLEs and distributed-order diffusion-type equation.
For different forms of the weight functions we obtained ultraslow diffusion,
strong anomaly, and other complex diffusive behaviors.

It will be interesting to compare the results obtained here for the GLE to the
case for external noise klimo , i.e., when the second FDT is not satisfied. In
both cases superstatistical and stochastic variations of the diffusion
coefficient and anomalous scaling exponent (Hurst exponent) have been analyzed
recently diego1 ; chaos ; elife ; jakub ; jakub1 ; wei . Studying such
concepts in the frameworks developed here will significantly enlarge our
current range of stochastic models for disordered systems. We also note
potential generalizations with respect to subordinated generalized Langevin
equation models such as those studied in fox ; yingjie , fractional GLE
sandev3 ; sandev_jmp , as well as presence of stochastic resetting maj  in the
system, including resetting in the memory kernel fractal and fractional ;
TSbook .

## Acknowledgements

RM and TS acknowledge financial support by the German Science Foundation (DFG,
Grant number ME 1535/12-1). TS is supported by the Alliance of International
Science Organizations (Project No. ANSO-CR-PP-2022-05). TS is also supported
by the Alexander von Humboldt Foundation. KG and TP acknowledge financial
support provided by the NCN Grant Preludium Bis 2 No. UMO-2020/39/O/ST2/01563.

##  Appendix A Multinomial Prabhakar type functions ${\cal
E}_{(\vec{\mu}),\beta}(t;\vec{a})$

The multinomial Prabhakar function ${\cal E}_{(\vec{\mu}),\beta}(t;\vec{a})$
with $\mu_{1}>\mu_{2}>\ldots>\mu_{m}>0$, see Ref. EBazhlekova21 , is related
to the multinomial ML function
$E_{(\vec{\mu}),\beta}(-a_{1}t^{\mu_{1}},\ldots,-a_{m}t^{\mu_{m}})$ as follows

${\cal
E}_{(\vec{\mu}),\beta}(t;\vec{a})=t^{\beta-1}E_{(\vec{\mu}),\beta}(-a_{1}t^{\mu_{1}},\ldots,-a_{m}t^{\mu_{m}}).$
(103)

Its Laplace transform reads

$\mathcal{L}\\{{\cal
E}_{(\vec{\mu}),\beta}(t;\vec{a});s\\}=\frac{s^{\mu_{1}-\beta}}{s^{\mu_{1}}+a_{m}s^{\mu_{1}-\mu_{m}}+\ldots+a_{2}s^{\mu_{1}-\mu_{2}}+a_{1}}.$
(104)

Theorem 2.3 from Ref. EBazhlekova21  implies that the following identities
hold true:

$\displaystyle\int_{0}^{t}\mathcal{E}_{(\vec{\mu}),\beta}(\xi;\vec{a})d\xi$
$\displaystyle=\mathcal{E}_{(\vec{\mu}),\beta+1}(t;\vec{a}),$ (105)
$\displaystyle\frac{d}{dt}\mathcal{E}_{(\vec{\mu}),\beta}(t;\vec{a})$
$\displaystyle=\mathcal{E}_{(\vec{\mu}),\beta-1}(t;\vec{a}),\quad\beta>1.$
(106)

The asymptotic behaviors of the multinomial Prabhakar function for $t\ll 1$
and $t\gg 1$ are equal to

$\displaystyle\mathcal{E}_{(\vec{\mu}),\beta}(t;\vec{a})$
$\displaystyle\sim_{t\rightarrow
0}\frac{t^{\beta-1}}{\Gamma(\beta)}-\sum_{j=1}^{m}\frac{t^{\beta-1+\mu_{j}}}{\Gamma(\beta+\mu_{j})},$
(107) $\displaystyle\mathcal{E}_{(\vec{\mu}),\beta}(t;\vec{a})$
$\displaystyle\sim_{t\to\infty}a_{2}^{-1}t^{\beta-\mu_{2}-1}E_{\mu_{1}-\mu_{2},\beta-\mu_{2}}(-a_{1}a_{2}^{-1}t^{\mu_{1}-\mu_{2}}),$
(108)

respectively. Eq. (103) implies the series representation of ${\cal
E}_{(\vec{\mu}),\beta}(s;\vec{a})$, that is

$\mathcal{E}_{(\vec{\mu}),\beta}(t;\vec{a})=\sum_{k=0}^{\infty}\sum_{k_{1}+\ldots+k_{m}=k\atop
k_{1},\ldots,k_{m}\geq 0}\frac{(-1)^{k}\,k!}{k_{1}!\ldots
k_{m}!}\,\frac{\prod_{j=1}^{m}a_{j}^{k_{j}}t^{\beta-1+\sum_{j=0}^{m}\mu_{j}k_{j}}}{\Gamma(\beta+\sum_{j=1}^{m}\mu_{j}k_{j})}$
(109)

from which it appears that we can represent the multinomial Prabhakar function
by the sums of the three-parameter ML function Prabhakar , i.e.,

$E_{\mu,\beta}^{\delta}(-\lambda
t^{\mu})=\frac{1}{\Gamma(\delta)}\,\sum_{r=0}^{\infty}\frac{\Gamma(\delta+r)\,(-\lambda
t^{\mu})^{r}}{r!\Gamma(\beta+\mu r)},$ (110)

whose the Laplace transform reads

$\mathcal{L}\\{t^{\beta-1}E_{\mu,\beta}^{\delta}(-\lambda
t^{\mu});s\\}=\frac{s^{\mu\delta-\beta}}{(\lambda+s^{\mu})^{\delta}}$ (111)

For $\delta=1$ the three-parameter ML function becomes a two-parameter ML
function, while for $\beta=\delta=1$ it becomes one parameter ML function. The
asymptotic expansion of the three-parameter ML function follows from the
expression RGarra18 ; AGusti20

$E_{\rho,\beta}^{\delta}(-z)=\frac{z^{-\delta}}{\Gamma(\delta)}\sum_{n=0}^{\infty}\frac{\Gamma(\delta+n)}{\Gamma(\beta-\rho(\delta+n))}\frac{(-z)^{-n}}{n!},$
(112)

with $z>1$, and $0<\rho<2$.

For example, for $m=2$, we have

$\displaystyle{\cal E}_{(\mu_{1},\mu_{2}),\beta}(t;a_{1},a_{2})$
$\displaystyle=\sum_{j=0}^{\infty}(-a_{2})^{j}t^{\beta+\mu_{2}j-1}E^{1+j}_{\mu_{1},\beta+\mu_{2}j}(-a_{1}t^{\mu_{1}})$
(113)
$\displaystyle=\sum_{j=0}^{\infty}(-a_{1})^{j}t^{\beta+\mu_{1}j-1}E^{1+j}_{\mu_{2},\beta+\mu_{1}j}(-a_{2}t^{\mu_{2}}).$
(114)

Moreover, Eq. (104) for $m=2$ under condition
$|a_{1}s^{-\mu_{1}}+a_{2}s^{-\mu_{2}}|>1$ gives

$\displaystyle\mathcal{E}_{(\mu_{1},\mu_{2}),\beta}(t;a_{1},a_{2})$
$\displaystyle=\sum_{r=0}^{\infty}\frac{(-1)^{r}}{a_{1}^{r+1}}t^{\beta-\mu_{2}(r+1)-1}E_{\mu_{1}-\mu_{2},\beta-\mu_{2}(r+1)}^{r+1}\Big{(}-\frac{a_{2}}{a_{1}}t^{\mu_{1}-\mu_{2}}\Big{)}$
(115)
$\displaystyle=\sum_{r=0}^{\infty}\frac{(-1)^{r}}{a_{2}^{r+1}}t^{\beta-\mu_{1}(r+1)-1}E_{\mu_{2}-\mu_{1},\beta-\mu_{1}(r+1)}^{r+1}\left(-\frac{a_{1}}{a_{2}}t^{\mu_{2}-\mu_{1}}\right).$
(116)

From Eq. (109) for $m=3$ after some laborious calculation can be derived,
e.g., following formula which is used in our manuscript

$\displaystyle{\cal
E}_{(\mu_{1},\mu_{2},\mu_{3}),\beta}(t;a_{1},a_{2},a_{3})=\sum_{j=0}^{\infty}\frac{(-a_{3})^{j}}{j!}\sum_{k=0}^{\infty}\frac{(-a_{1})^{k}}{k!}(k+j)!t^{\beta-1+\mu_{1}k+\mu_{3}j}E^{j+k+1}_{\mu_{2},\beta+\mu_{1}k+\mu_{3}j}(-a_{2}t^{\mu_{2}}).$
(117)

##  Appendix B Tauberian theorems Feller

If the Laplace transform pair $\hat{r}(s)$ of the function $r(t)$ behaves like

$\hat{r}(s)\simeq s^{-\rho}L(s^{-1}),\quad s\to 0,\quad\rho>0,$ (118)

where $L(t)$ is a slowly varying function at infinity, then $r(t)$ has the
following asymptotic behavior Feller

$r(t)\simeq\frac{1}{\Gamma(\rho)}t^{\rho-1}L(t),\quad t\to\infty.$ (119)

A slowly varying function at infinity means that

$\lim_{t\to\infty}\frac{L(at)}{L(t)}=1,\quad a>0.$ (120)

The Tauberian theorem works also for the opposite asymptotic, i.e., for $t\to
0$.

##  Appendix C The Volterra family functions

Volterra’s function is defined as follows AApelblat10 ; AApelblat13 ; Erdei

$\mu(t,\beta,\alpha)=\frac{1}{\Gamma(1+\beta)}\int_{0}^{\infty}\frac{t^{u+\alpha}\,u^{\beta}}{\Gamma(u+\alpha+1)}du,\qquad\Re(\beta)>-1\quad\text{and}\quad
t>0,$ (121)

whose particular cases are

$\displaystyle\alpha=\beta=0:\qquad$ $\displaystyle\nu(t)=\mu(t,0,0),$
$\displaystyle\alpha\neq 0,\,\beta=0:\qquad$
$\displaystyle\nu(t,\alpha)=\mu(t,0,\alpha),$
$\displaystyle\alpha=0,\,\beta\neq 0:\qquad$
$\displaystyle\mu(t,\beta)=\mu(t,\beta,0).$

The Laplace transform of the Volterra’s function $\mu(t,\beta,\alpha)$ is
given by Erdei

$\mathcal{L}\\{\mu(t,\beta,\alpha);s\\}=\frac{1}{s^{\alpha+1}\log^{\beta+1}s}.$
(122)

##  Appendix D Derivation of Eq. (83)

Since,

$\displaystyle\int_{0}^{\infty}e^{-st}t^{a}(\log^{\lambda-1}t)dt$
$\displaystyle=\frac{\partial^{\lambda-1}}{\partial
a^{\lambda-1}}\int_{0}^{\infty}e^{-st}t^{a}dt$
$\displaystyle=\frac{\partial^{\lambda-1}}{\partial
a^{\lambda-1}}\left\\{\frac{\Gamma(a+1)}{s^{a+1}}\right\\}=I^{1-\lambda}_{0}\left\\{\frac{\Gamma(a+1)}{s^{a+1}}\right\\}$
$\displaystyle=\frac{\Gamma(a+1)}{\Gamma(1-\lambda)}\int^{s}_{0}(s-t)^{-\lambda}t^{-a-1}dt=\frac{\Gamma(a+1)\Gamma(-a)}{s^{\lambda+a}\Gamma(1-\lambda-a)},$
(123)

for $a=-1/2$, we get

$\int_{0}^{\infty}e^{-st}t^{-1/2}(\log^{\lambda-1}t)dt=\frac{\pi}{s^{\lambda-1/2}\Gamma(3/2-\lambda)},\quad
0<\lambda<1.$ (124)

Then,

$\displaystyle k_{B}T\,\hat{\gamma}(s)$
$\displaystyle=\int_{0}^{1}\Gamma(3/2-\lambda)d\lambda\left[\int_{0}^{\infty}e^{-st}t^{-1/2}(\log^{\lambda-1}t)dt\right]$
$\displaystyle=\pi\sqrt{s}\int_{0}^{1}\frac{d\lambda}{s^{\lambda}}=\pi\frac{s-1}{\sqrt{s}\log
s}.$ (125)

## References

  * (1) A. Apelblat, Integral transforms and Volterra function (Nova Sci. Pub., Inc.: New York, 2010).
  * (2) A. Apelblat, Volterra functions (Nova Sci. Pub., Inc.: New York, 2013).
  * (3) H. Bateman, A. Erdélyi, Higher Transcendental Functions, vol 3, (Robert E. Krieger Pub. Comp.: Malabra, 1955).
  * (4) M. Balcerek, A. Wyłomańska, K. Burnecki, R. Metzler, and D. Krapf, E-print arXiv.2307.12919.
  * (5) E. Barkai, Y. Garini, and R. Metzler, Phys. Today 65(8), 29 (2012).
  * (6) E. Bazhlekova, Fract. Calc. Appl. Anal. 24, 88 (2021).
  * (7) S. Burov and E. Barkai, Phys. Rev. E 78, 031112 (2008).
S. Burov and E. Barkai, Phys. Rev. Lett. 100, 070601 (2008).

  * (8) S. Burov, J. -H. Jeon, R. Metzler, and E. Barkai, Phys. Chem. Chem. Phys. 13, 1800 (2011).
  * (9) P. Castiglione, A. Mazzino, P. Muratore-Ginanneschi, and A. Vulpiani, Physica D 134, 75 (1999).
  * (10) E. Capelas de Oliveira, F. Mainardi, and J. Vaz Jr., Eur. Phys. J. 193, 161 (2011).
  * (11) A. V. Chechkin, R. Gorenflo, and I. M. Sokolov, Phys. Rev. E 66, 046129 (2002).
  * (12) A. V. Chechkin, J. Klafter, and I. M. Sokolov, Europhys. Lett. 63, 326 (2003).
  * (13) M. Balcerek, K. Burnecki, S. Thapa, A. Wyłomańska, and A. Chechkin, Chaos 32. 093114 (2022).
  * (14) W. T. Coffey and Y. P. Kalmykov, The Langevin Equation (World Scientific: Singapore, 2012).
  * (15) M. A. Despósito and A. D. Viñales, Phys. Rev. E 77, 031123 (2008).
M. A. Despósito and A. D. Viñales, Phys. Rev. E 80, 021111 (2009).

  * (16) J. Dräger and J. Klafter, Phys. Rev. Lett. 84, 5998 (2000).
  * (17) C. H. Eab and S. C. Lim, Phys. Rev. E 83, 031136 (2011).
  * (18) A. Einstein, Ann. Phys. (Leipzig) 322, 549 (1905).
  * (19) A. Erdelyi, W. Magnus, F. Oberhettinger, and F. G. Tricomi, Higher Transcedential Functions, Vol. 3 (McGraw-Hill: New York, 1955).
  * (20) M. R. Evans and S. N. Majumdar, Phys. Rev. Lett. 106, 160601 (2011).
M. R. Evans, S. N. Majumdar, and G. Schehr, J. Phys. A: Math. Theor. 53,
193001 (2020).

  * (21) W. Feller, An Introduction to Probability Theory and Its Applications, Vol. II (Wiley: New York, 1968).
  * (22) Z. R. Fox, E. Barkai, and D. Krapf, Nat. Commun. 12, 6162 (2021).
  * (23) T. Franosch, M. Grimm, M. Belushkin, F. M. Mor, G. Foffi, L. Forro, and S. Jeney, Nature 478, 7367 (2011).
  * (24) R. Garra and R. Garrappa, Commun. Nonlinear Sci. Numer. Sim. 56, 314 (2018).
  * (25) A. Godec, A. V. Chechkin, E. Barkai, H. Kantz, and R. Metzler, J. Phys. A: Math. Theor. 47, 492002 (2014).
  * (26) R. Gorenflo and F. Mainardi, J. Phys.: Conf. Ser. 7, 1 (2005).
  * (27) D. Han, N. Korabel, R. Chen, M. Johnston, A. Gavrilova, V. J. Allan, S. Fedotov, and T. A. Waigh, eLife 9, e52224 (2020).
  * (28) K. Górska, T. Pietrzak, T. Sandev, Ž. Tomovski, J. Comput. Appl. Math. 433, 115306 (2023).
  * (29) K. Górska, A. Horzela, L. Bratek, G. Dattoli, and K. A. Penson, J. Phys. A: Math. Theor. 51, 135202 (2018).
  * (30) K. Górska, A. Horzela, and K. A. Penson, J. Phys. A: Math. Theor. 56, 313001 (2023).
  * (31) I. Goychuk, V. I. Kharchenko, and R. Metzler, Phys. Rev. E 96, 052134 (2017).
  * (32) A. Giusti, I. Colombaro, R. Garra, R. Garrappa, F. Polito, M. Popolizio, and F. Mainardi, Frac. Calc. Appl. Anal. 23, 9 (2020).
  * (33) I. Goychuk, Adv. Chem. Phys. 150, 187 (2012).
  * (34) Y. He, S. Burov, R. Metzler, and E. Barkai, Phys. Rev. Lett. 101, 058101 (2008).
  * (35) F. Höfling and T. Franosch, Rep. Progr. Phys. 76, 046602 (2013).
  * (36) J. -H. Jeon, V. Tejedor, S. Burov, E. Barkai, C. Selhuber-Unkel, K. Berg-Sørensen, L. Oddershede, and R. Metzler, Phys. Rev. Lett. 106, 048103 (2011).
  * (37) J. -H. Jeon and R. Metzler, Phys. Rev. E 85, 021147 (2012).
  * (38) J. -H. Jeon, N. Leijnse, L. B. Oddershede, and R. Metzler, New J. Phys. 15, 045011 (2013).
  * (39) J. -H. Jeon, H. Martinez-Seara Monne, M. Javanainen, and R. Metzler, Phys. Rev. Lett. 109, 188103 (2012).
  * (40) Y. L. Klimontovich, Turbulent Motion and the Structure of Chaos (Kluwer: Dordrecht, 1991).
  * (41) A. N. Kochubei, J. Math Anal. Appl. 340, 252 (2008).
  * (42) A. N. Kochubei, Ukranian Math. J. 60, 551 (2008).
  * (43) A. N. Kochubei, J. Phys. A: Math. Theor 42, 315203 (2009).
  * (44) A. N. Kochubei, Integr. Equ. Oper. Theory 71, 583 (2011).
  * (45) S. C. Kou and X. S. Xie, Phys. Rev. Lett. 93, 180603 (2004).
  * (46) D. Krapf and R. Metzler, Phys. Today 72(9), 48 (2019).
  * (47) R. Kubo, Rep. Prog. Phys. 29, 255 (1966).
  * (48) R. Kupferman, J. Stat. Phys., 114, 291 (2004).
  * (49) P. Langevin, Comptes Rendus 146, 530 (1908).
  * (50) Y. Liang, W. Wang, and R. Metzler, Phys. Rev. E, at press; E-print arXiv:2302.04872.
  * (51) S. -D. Lin and H. M. Srivastava, Appl. Math. Comput. 154, 725 (2004).
  * (52) L. Lizana, T. Ambjörnsson, A. Taloni, E. Barkai, and M. L. Lomholt, Phys. Rev. E 81, 051118 (2010).
  * (53) E. Lutz, Phys. Rev. E 64, 051106 (2001).
  * (54) F. Mainardi, Discrete Contin. Dyn. Syst. - B 19, 2267 (2014).
  * (55) F. Mainardi and P. Pironi Extr. Math. 10, 140 (1996).
  * (56) F. Mainardi, Fractional Calculus and Waves in Linear Viscoelesticity: An Introduction to Mathematical Models (Imperial College Press: London, 2010).
  * (57) F. Mainardi, A. Mura, R. Gorenflo, and M. Stojanovic J. Vib. Control 13, 1249 (2007).
  * (58) R. Metzler, E. Barkai, and J. Klafter, Phys. Rev. Lett. 82, 3563 (1999).
  * (59) R. Metzler and J. Klafter, Phys. Rep. 339, 1 (2000).
R. Metzler and J. Klafter, J. Phys. A: Math. Gen. 37, R161 (2004).

  * (60) R. Metzler, J. -H. Jeon, A. G. Cherstvy, and E. Barkai, Phys. Chem. Chem. Phys. 16, 24128 (2014).
  * (61) D. Molina-Garcia, T. Sandev, H. Safdari, G. Pagnini, A. Chechkin, and R. Metzler, New J. Phys. 20, 103027 (2018).
  * (62) G. Muñoz-Gil, G. Volpe, M. A. Garcia-March, E. Aghion, A. Argun, C. B. Hong, T. Bland, S. Bo, J. A. Conejero, N. Firbas, Ò. Garibo i Orts, A. Gentili, Z. Huang, J. -H. Jeon, H. Kabbech, Y. Kim, P. Kowalek, D. Krapf, H. Loch-Olszewska, M. A. Lomholt, J. -B. Masson, P. G. Meyer, S. Park, B. Requena, I. Smal, T. Song, J. Szwabiński, S. Thapa, H. Verdier, G. Volpe, A. Widera, M. Lewenstein, R. Metzler, and C. Manzo, Nature Comm. 12, 6253 (2021).
  * (63) K. Nørregaard, R. Metzler, C. Ritter, K. Berg-Sørensen, and L. Oddershede, Chem. Rev. 117, 4342 (2017).
  * (64) I. Petreska, L. Pejov, T. Sandev, L. Kocarev, and R. Metzler, Fractal Fract. 6, 88 (2022).
  * (65) H. Pollard, Bull. Amer. Math. Soc. 54, 1115 (1948).
  * (66) N. Pottier, Physica A 317, 371 (2003).
N. Pottier and A. Mauger, Physica A 282, 77 (2000).

  * (67) T. R. Prabhakar, Yokohama Math. J. 19, 7 (1971).
  * (68) T. Sandev, Physica Macedonica 61, 59 (2012).
  * (69) T. Sandev and A. Iomin, Special Functions of Fractional Calculus: Applications to Diffusion and Random Search Processes (World Scientific: Singapore, 2022).
  * (70) T. Sandev, R. Metzler, and Ž. Tomovski, Fract. Calc. Appl. Anal. 15, 426–450 (2012).
  * (71) T. Sandev, R. Metzler, and Ž. Tomovski, J. Math. Phys. 55, 023301 (2014).
  * (72) T. Sandev, Ž. Tomovski, and J. L. A. Dubbeldam, Physica A 390, 3627 (2011).
  * (73) T. Sandev and Ž. Tomovski, Phys. Lett. A 378, 1 (2014).
  * (74) H. Seckler and R. Metzler, Nature Comm. 13, 6717 (2022).
  * (75) C. Selhuber-Unkel, M. Weiss, S. Scott, Y. F. Barooji, A. Sabri, J. T. Erler, R. Metzler, and L. B. Oddershede, Phys. Chem. Chem. Phys. 25, 1513 (2023).
  * (76) Y. G. Sinai, Theor. Probab. Appl. 27, 256 (1982).
  * (77) J. Ślȩzak, R. Metzler, and M. Magdziarz, New J. Phys. 20, 023026 (2018).
  * (78) J. Ślȩzak and R. Metzler, J. Phys. A 56, 35LT01 (2023).
  * (79) M. von Smoluchowski, Ann. Phys. (Leipzig) 21, 756 (1906).
  * (80) I. M. Sokolov, Soft Matter 8, 9043 (2012).
  * (81) J. Klafter and I. M. Sokolov, Physics World 18, 29 (2005).
  * (82) J. Spiechowicz, I. G. Marchenko, P. Hänggi, and J. Łczka, Entropy 25, 42 (2022).
  * (83) W. Sutherland, Philos. Mag. 9, 781 (1905).
  * (84) A. A. Tateishi, E. K. Lenzi, L. R. da Silva, H. V. Ribeiro, S. Picoli Jr., and R. S. Mendes, Phys. Rev. E 85, 011147 (2012).
  * (85) E. van der Straeten and C. Beck, Physica A 390, 951 (2011).
  * (86) O. Vilk, E. Aghion, T. Avgar, C. Beta, O. Nagel, A. Sabri, R. Sarfati, D. K. Schwartz, M. Weiss, D. Krapf, R. Nathan, R. Metzler, and M. Assaf, Phys. Rev. Res. 4, 033055 (2022).
  * (87) A. D. Viñales and M. A. Despósito, Phys. Rev. E 75, 042102 (2007).
A. D. Viñales, K. G. Wang and M. A. Despósito, Phys. Rev. E 80, 011101 (2009).
A. D. Viñales and M. A. Despósito, Phys. Rev. E 73, 016111 (2006).

  * (88) T. Vojta, S. Skinner, and R. Metzler, Phys. Rev. E 100, 042142 (2019).
  * (89) K. G. Wang, Phys. Rev. A 45, 833 (1992).
K. G. Wang and M. Tokuyama, Physica A 265, 341 (1999).

  * (90) W. Wang, M. Balcerek, K. Burnecki, A. V. Chechkin, S. Janušonis, J. Ślȩzak, T. Vojta, A. Wyłomańska, and R. Metzler, Phys. Rev. Res. 5, L032025 (2023).
  * (91) H. Yang, G. Luo, P. Karnchanaphanurach, T.-M. Louie, I. Reich, S. Cova, L. Xun, and X. S. Xie, Science 302, 262 (2003).
  * (92) R. Zwanzig, Nonequilibrium Statistical Mechanics (Oxford Univ. Press: New York, 2001).