Localization and ballistic diffusion for the tempered fractional Brownian-Langevin motion

TL;DR

This paper investigates the properties of tempered fractional Brownian motion and Langevin equations, revealing transitions between localization, ballistic, and anomalous diffusion behaviors over different time scales.

Contribution

It introduces the tempered fractional Langevin equation (tfLe), analyzing its diffusion transitions and ergodic properties, extending understanding of tempered fractional stochastic processes.

Findings

Mean squared displacement transitions from ballistic to anomalous and back to ballistic diffusion.

Overdamped tfLe transitions from subdiffusive to ballistic behavior.

Harmonic potential influences the diffusion transition dynamics.

Abstract

This paper further discusses the tempered fractional Brownian motion, its ergodicity, and the derivation of the corresponding Fokker-Planck equation. Then we introduce the generalized Langevin equation with the tempered fractional Gaussian noise for a free particle, called tempered fractional Langevin equation (tfLe). While the tempered fractional Brownian motion displays localization diffusion for the long time limit and for the short time its mean squared displacement has the asymptotic form $t^{2H}$, we show that the asymptotic form of the mean squared displacement of the tfLe transits from $t^2$ (ballistic diffusion for short time) to $t^{2-2H}$, and then to $t^2$ (again ballistic diffusion for long time). On the other hand, the overdamped tfLe has the transition of the diffusion type from $t^{2-2H}$ to $t^2$ (ballistic diffusion). The tfLe with harmonic potential is also considered.

Equations

$$m\frac{\text{rmd}v}{\text{rmd}t}=-\xi v+F(t),$$

$$B_H(t)=\frac{1}{\Gamma (1-\alpha )}{\int }_{-\infty }^{+\infty }[(t-s{)}_{+}^{-\alpha }-(-s{)}_{+}^{-\alpha }]B(\text{rmd}s),$$

$$(x)_{+}=\cases{x&for $x>0$\\ 0&for $x\leq 0$\\ }$$

$$m\frac{{\text{rmd}}^{2}x(t)}{\text{rmd}{t}^2}=-\xi {\int }_0^t(t-\tau {)}^{2H-2}\frac{\text{rmd}x}{\text{rmd}\tau }\text{rmd}\tau +\eta \gamma (t),$$

$$B_{\alpha ,\lambda }(t)={\int }_{-\infty }^{+\infty }[{\text{rme}}^{-\lambda (t-x{)}_{+}}(t-x{)}_{+}^{-\alpha }-{\text{rme}}^{-\lambda (-x{)}_{+}}(-x{)}_{+}^{-\alpha }]B(\text{rmd}x),$$

$$\{B_{\alpha ,\lambda }(ct){\}}_{t\in \mathbb{R}}=\{c^H{B}_{\alpha ,c\lambda }(t){\}}_{t\in \mathbb{R}}$$

$$\text{Cov}[B_{\alpha ,\lambda }(t),B_{\alpha ,\lambda }(s)]=\frac{1}{2}\left[C_t^2|t{|}^{2H}+C_s^2|s{|}^{2H}-C_{t-s}^2|t-s{|}^{2H}\right]$$

$$C_t^2$$

$$P(x,t)=\frac{1}{\sqrt{2\pi C_t^2{t}^{2H}}}{\text{rme}}^{-\frac{x^2}{2C_t^2{t}^{2H}}}.$$

$$\frac{\partial P(x,t)}{\partial t}=-\frac{\Gamma (H+\frac{1}{2})}{\sqrt{\pi }(2\lambda {)}^H}\left[H{t}^{H-1}{K}_H(\lambda t)+t^H{\dot{K}}_H(\lambda t)\right]\frac{{\partial }^{2}P(x,t)}{\partial x^2}.$$

$$\frac{\partial P(x,t)}{\partial t}=\frac{\Gamma (H+\frac{1}{2}){\text{rme}}^{-\lambda t}}{(2\lambda {)}^{H+\frac{1}{2}}}\left[\lambda t^{H-\frac{1}{2}}-\left(H-\frac{1}{2}\right)t^{H-\frac{3}{2}}\right]\frac{{\partial }^{2}P(x,t)}{\partial x^2}.$$

$$E_B(B_{\alpha ,\lambda }(t))=\frac{\text{Var}[{\bar{\delta }}^2(B_{\alpha ,\lambda }(t))]}{⟨{\bar{\delta }}^2(B_{\alpha ,\lambda }(t)){⟩}^2}=\frac{⟨[{\bar{\delta }}^2(B_{\alpha ,\lambda }(t)){]}^2⟩-⟨{\bar{\delta }}^2(B_{\alpha ,\lambda }(t)){⟩}^2}{⟨{\bar{\delta }}^2(B_{\alpha ,\lambda }(t)){⟩}^2}.$$

$$⟨{\bar{\delta }}^2(B_{\alpha ,\lambda }(t))⟩=C_{\Delta }^2|\Delta {|}^{2H},$$

$$\text{Var}[{\bar{\delta }}^2(B_{\alpha ,\lambda }(t))]\simeq D\frac{4{\Gamma }^2(H+\frac{1}{2})}{(2\lambda {)}^{2H+1}}{\Delta }^{2H}{t}^{-1},$$

$$\gamma (t)=\frac{B_{\alpha ,\lambda }(t+h)-B_{\alpha ,\lambda }(t)}{h},$$

$$⟨\gamma (t_1)\gamma (t_2)⟩=\frac{1}{h^2}⟨(B_{\alpha ,\lambda }(t_1+h)-B_{\alpha ,\lambda }(t_1))(B_{\alpha ,\lambda }(t_2+h)-B_{\alpha ,\lambda }(t_2))⟩$$

$$=\frac{1}{2h^2}(C_{t_1-t_2+h}^2|t_1-t_2+h{|}^{2H}+C_{t_1-t_2-h}^2|t_1-t_2-h{|}^{2H}-2C_{t_1-t_2}^2|t_1-t_2{|}^{2H}),$$

$$m\frac{{\text{rmd}}^{2}x(t)}{\text{rmd}{t}^2}=-\xi {\int }_0^{t}K(t-\tau )\frac{\text{rmd}x}{\text{rmd}\tau }\text{rmd}\tau +F(t),$$

$$x(s)=\frac{F(s)+m{v}_0}{m{s}^2+\xi sK(s)}.$$

$$x(s)\simeq \frac{F(s)+m{v}_0}{m{s}^2+\xi s^{3}f(s)}$$

$$v(s)\simeq \frac{F(s)+m{v}_0}{ms+\xi s^{2}f(s)}.$$

$$x(s)\simeq \frac{F(s)+m{v}_0}{A{s}^2}$$

$$x(t)\simeq \frac{1}{A}{\int }_0^t\tau F(t-\tau )\text{rmd}\tau +\frac{m{v}_0}{A}t.$$

$$⟨x^2(t)⟩\simeq \frac{2}{A^2}{\int }_0^t{\int }_0^{t_1}\text{rmd}{t}_2\text{rmd}{t}_1{t}_1{t}_2⟨F(t-t_1)F(t-t_2)⟩+\frac{m^2{v}_0^2}{A^2}{t}^2,$$

$$⟨x^2(t)⟩$$

$$⟨(\Delta x{)}^2⟩=⟨[x(t)-⟨x(t)⟩{]}^2⟩\simeq \frac{2k_{B}T\xi \Gamma (2H)}{A^2(2\lambda {)}^{2H}}{t}^2,$$

$$⟨F(s_1)F(s_2)⟩=k_{B}T\xi \frac{s_1^{2}f(s_1)+s_2^{2}f(s_2)}{s_1+s_2}.$$

$$⟨x(s_1)x(s_2)⟩$$

$$⟨x(t_1)x(t_2)⟩\simeq \left(\frac{2k_{B}T\xi \Gamma (2H)}{A^2(2\lambda {)}^{2H}}+\frac{m^2{v}_0^2}{A^2}\right)t_1{t}_2,$$

$$v(t)$$

$$⟨v(t_1)v(t_2)⟩$$

Full text

Keywords: generalized Langevin equation, tempered fractional Brownian motion, ballistic diffusion, localization diffusion, harmonic potential

(Some figures may appear in colour only in the online journal)

1 Introduction

Revealing the mechanism of the motion of particles in complex disordered systems is a fundamental and challenging topic [1, 2]. Generally, the motion is no longer Brownian because of the complex environment and/or the properties of the particles themselves. The mean squared displacement (MSD) of Brownian motion goes like $\langle(\Delta x)^{2}\rangle=\langle[x(t)-\langle x(t)\rangle]^{2}\rangle\sim t^{\nu}$ with $\nu=1$; for a long time $t$ if $\nu\not=1$, it is called anomalous diffusion, being subdiffusion for $\nu<1$ and superdiffusion for $\nu>1$ [3]; in particular, it is termed as localization diffusion if $\nu=0$ and ballistic diffusion if $\nu=2$. There are two types of stochastic processes to model anomalous diffusion: Gaussian processes and non-Gaussian ones; the non-Gaussian processes [4, 5, 6, 7] include continuous time random walks (CTRWs), Lévy processes, subordinated Lévy processes, and the Gaussian ones [8, 9, 10, 11, 12] contain fractional Brownian motion (fBm), generalized Langevin equation, etc. In this paper, we mainly consider the Langevin type equation with correlated internal noise, being a Gaussian process.

The most basic Gaussian process, describing normal diffusion, is Brownian motion with its corresponding Langevin equation [13]

[TABLE]

where $v$ is the velocity of a Brownian particle with a mass $m$, $\xi$ is a frictional constant, and the random fluctuation force $F(t)$ is white noise. For large $t$, the mean squared displacement of a Brownian particle is $\langle(\Delta x(t))^{2}\rangle=\langle(\int_{0}^{t}v(s)\rmd s)^{2}\rangle\simeq 2Dt$ with $D=\frac{k_{B}T}{\xi}$ being the Einstein relation, where $k_{B}$ is the Boltzmann constant and $T$ the absolute temperature of the environment. As the extension of Brownian motion, fractional Brownian motion (fBm) is still a Gaussian process, which can be seen as the fractional derivative (or integral) of a Brownian motion. The fBm is defined as [14]

[TABLE]

where

[TABLE]

with $-\frac{1}{2}<\alpha<\frac{1}{2}$, and the Hurst index $H=\frac{1}{2}-\alpha$. Note that Brownian motion is recovered when $H=\frac{1}{2}$. The variance of $B_{H}(t)$ is $2D_{H}t^{2H}$, where $D_{H}=[\Gamma(1-2H)\cos(H\pi)]/(2H\pi)$. Fractional Langevin equation [10, 12], still being Gaussian process, reads

[TABLE]

where $x(t)$ is the particle displacement, $\eta=[k_{B}T\xi/(2D_{H}H(2H-1))]^{1/2}$, $\gamma(t)=\rmd B_{H}(t)/\rmd t$ is the fractional Gaussian noise, and $1/2<H<1$ is the Hurst parameter. The mean squared displacement of the trajectory sample $x(t)$ for large $t$ is $\langle x^{2}(t)\rangle\simeq 2k_{B}T/(\xi\Gamma(2H-1)\Gamma(3-2H))t^{2-2H}$.

Along the direction of extension of Brownian motion, the new Gaussian process, called tempered fractional Brownian motion (tfBm) [15], was introduced by Meerschaert and Sabzikar with its definition

[TABLE]

where $\lambda>0$, $\alpha<\frac{1}{2}$, and the Hurst index $H=\frac{1}{2}-\alpha$; and its basic theory was developed with application to modeling wind speed. This paper naturally introduces the generalized Langevin equation, termed as tempered fractional Langevin equation (tfLe), with the tempered fractional noise as internal noise. We discuss the ergodicity of the tfBm, and derive the corresponding Fokker-Planck equation. The mean squared displacement (MSD) of the tfLe is carefully analyzed, displaying the transition from $t^{2}$ (ballistic diffusion for short time) to $t^{2-2H}$, and then to $t^{2}$ (again ballistic diffusion for long time). For the overdamped tfLe, its MSD transits from $t^{2-2H}$ to $t^{2}$ (ballistic diffusion). The properties of the correlation function of the tfLe with harmonic potential is also considered.

The outline of this paper is as follows. In Section 2, we review the tfBm, derive its Fokker-Planck equation, and discuss its ergodicity. The tfLe is introduced in Section 3, in which the underdamped case, overdamped case, and the tfLe with harmonic potential are discussed in detail. We conclude the paper with some discussions in the last section.

2 Tempered fractional Brownian motion and tempered fractional Gaussian noise

We introduce the definitions of the tfBm and the tempered fractional Gaussian noise (tfGn). The Fokker-Planck equation for tfBm is derived, and the ergodicity of the tfBm is discussed.

2.1 Tempered fractional Brownian motion

Tempered fractional Brownian motion is defined in (1.1), which modifies the power law kernel in the moving average representation of a fractional Brownian motion by adding an exponential tempering [15]; it is generalizedly self-similar in the sense that for any $c>0$

[TABLE]

in distribution and it has the covariance function

[TABLE]

for any $t,s\in\mathbb{R}$, where

[TABLE]

for $t\neq 0$ and $C_{0}^{2}=0$, where $K_{H}(x)$ is the modified Bessel function of the second kind [15]. It is obvious that the variance of tfBm is $\langle B_{\alpha,\lambda}^{2}(t)\rangle=C_{t}^{2}|t|^{2H}\simeq 2\Gamma(2H)(2\lambda)^{-2H}$ as $t\rightarrow\infty$ on account of $K_{H}(t)\simeq\sqrt{\pi}(2t)^{-\frac{1}{2}}\rme^{-t}$, which means that tfBm is localization. And $\langle B_{\alpha,\lambda}(t)\rangle=0$. Since $B_{\alpha,\lambda}(t)$ is a Gaussian process, when $t>0$, $B_{\alpha,\lambda}(t)\sim N(0,C_{t}^{2}t^{2H})$. Then it has the probability density function (PDF)

[TABLE]

The Fourier transform of (2.4) is $P(k,t)=\int_{-\infty}^{+\infty}\rme^{\rmi kx}P(x,t)\rmd x=\rme^{-C_{t}^{2}t^{2H}k^{2}/2}$; taking partial derivative w.r.t. $t$ and performing inverse Fourier transform on both sides of this equation lead to

[TABLE]

For large $t$, from (2.5), we have that the PDF is asymptotically governed by

[TABLE]

The non-Markovian property is implied by the time-dependent diffusion constant: $D_{\alpha,\lambda}(t)=(2\lambda)^{-H-1/2}\Gamma(H+1/2)\rme^{-\lambda t}[\lambda t^{H-1/2}-(H-1/2)t^{H-3/2}]$.

Next we consider the ergodicity of tfBm and the convergence speed of the variance of the time-average mean-squared displacement, $\bar{\delta}^{2}(B_{\alpha,\lambda}(t))=\int_{0}^{t-\Delta}[B_{\alpha,\lambda}(t^{\prime}+\Delta)-B_{\alpha,\lambda}(t^{\prime})]^{2}\rmd t^{\prime}/(t-\Delta)$, where $\Delta$ is the lag time. If the average of $\bar{\delta}^{2}(B_{\alpha,\lambda}(t))$ equals to the ensemble average of $B_{\alpha,\lambda}(t)$ and the variance of $\bar{\delta}^{2}(B_{\alpha,\lambda}(t))$ tends to zero when the measurement time is long, the process tfBm is ergodic [10]; it is indeed (see (2.7) and (2.8) and their derivations presented in A). In fact, the variance of $\bar{\delta}^{2}(B_{\alpha,\lambda}(t))$ is a measure of ergodicity breaking and the ergodicity breaking parameter is defined as

[TABLE]

For tfBm, from A the average of $\bar{\delta}^{2}(B_{\alpha,\lambda}(t))$ is

[TABLE]

hence $\langle\bar{\delta}^{2}\rangle=\langle B_{\alpha,\lambda}^{2}\rangle$ for all times; for long times and moderate $\lambda$, the variance of $\bar{\delta}^{2}(B_{\alpha,\lambda}(t))$ is

[TABLE]

where $D=\int_{0}^{\infty}\rmd\tau[(1+\tau)^{H-\frac{1}{2}}\rme^{-\lambda\Delta(1+\tau)}+|\tau-1|^{H-\frac{1}{2}}\rme^{-\lambda|\Delta(\tau-1)|}-2\tau^{H-\frac{1}{2}}\rme^{-\lambda\Delta\tau}]^{2}$. Thus we have $E_{B}(B_{\alpha,\lambda}(t))\simeq 4D\Gamma^{2}(H+1/2)/[(2\lambda)^{2H+1}C_{\Delta}^{4}\Delta^{2H}]t^{-1}$, which tends to zero with the speed $t^{-1}$ for $H\in(0,1)$ as $t\rightarrow\infty$. The simulation result for ergodicity breaking parameter $E_{B}(B_{\alpha,\lambda}(t))$ is shown in Figure LABEL:fig:0. On the other hand, for the limit $\lambda t\rightarrow 0$, $C_{t}^{2}$ tends to a constant. In this case, (2.2) reduces to the covariance function of fractional Brownian motion (fBm). Then the evolution of the ergodicity breaking parameter $E_{B}(B_{\alpha,\lambda}(t))$ is recovered to that of fBm [10].

2.2 Tempered fractional Gaussian noise

Given a tfBm (1.1), we can define the tempered fractional Gaussian noise (tfGn)

[TABLE]

being similar to the definition of fractional Gaussian noise (fGn) [14], where $h$ is small and $h\ll t$.

Based on the covariance function (2.2) of $B_{\alpha,\lambda}(t)$ and its zero mean, one can easily obtain that the mean of tfGn is $\langle\gamma(t)\rangle=0$ and its covariance is

[TABLE]

which means that the tfGn is a stationary Gaussian process. For a fixed $\lambda$, the asymptotic behavior of the covariance function is $\langle\gamma(0)\gamma(t)\rangle\simeq-\Gamma(H+1/2)\lambda^{2}(2\lambda)^{-H-1/2}t^{H-1/2}\rme^{-\lambda t}$ for large $t$ and $\langle\gamma(0)\gamma(t)\rangle\simeq C_{t}^{2}H(2H-1)t^{2H-2}$ with a positive constant $C_{t}^{2}$ for small $t$ (for the details of derivation, see B).

3 Generalized Langevin equation with tempered fractional Gaussian noise

We discuss the dynamics of the generalized Langevin equation with tempered fractional Gaussian noise for free particles, including the underdamped and overdamped cases. The properties of the correlation function of the equation under external potential are also discussed.

3.1 Dynamical behaviors for free particles

3.1.1 Underdamped generalized Langevin equation

We know that for large $t$ the MSD of fBm is like $t^{2H}$ and the one of the corresponding Langevin equation is like $t^{2-2H}$. In the above section, it is shown that tfBm is localization, i.e., its MSD goes like $t^{0}$. Can we expect that the MSD of the tfLe goes like $t^{2}$ for large $t$? The answer is yes.

Based on the second fluctuation-dissipation theorem [17], which links the dissipation memory kernel $K(t)$ with the autocorrelation function of internal noise $F(t)$: $\langle F(t_{1})F(t_{2})\rangle=k_{B}T\xi K(t_{1}-t_{2})$, the tfLe can be written as

[TABLE]

with $\dot{x}(0)=v_{0}$, $x(0)=0$, where $v_{0}$ is the initial velocity, $F(t)=\sqrt{2k_{B}T\xi}\gamma(t)$ is the internal noise with tfGn $\gamma(t)$, and $K(t)=2\langle\gamma(0)\gamma(t)\rangle=\frac{1}{h^{2}}(C_{t+h}^{2}|t+h|^{2H}+C_{t-h}^{2}|t-h|^{2H}-2C_{t}^{2}|t|^{2H})$ with $0<H<1$. In this case, the fluctuation and dissipation stem from the same source and the system will finally reach the equilibrium state. Taking Laplace transform of (3.1) leads to

[TABLE]

For convenience, denoting $C_{t}^{2}t^{2H}$ by $f(t)$, then $K(s)=\mathcal{L}[K(t)]=\frac{1}{h^{2}}(\rme^{sh}f(s)+\rme^{-sh}f(s)-2f(s))$. No matter $s\rightarrow 0$ or $s\rightarrow\infty$, $sh$ is always small since $h\ll t$. So along with a Taylor’s series expansion, we have $K(s)\simeq s^{2}f(s),$ for $s\rightarrow 0$ or $s\rightarrow\infty$. That is to say, the asymptotic expression of $x(s)$ is

[TABLE]

and

[TABLE]

Based on (3.3) and (3.4), we do the dynamics analysis for different cases.

Case I ($t$ is large, i.e., $s\rightarrow 0$, and $H\in(0,1)$): For moderate $\lambda$, $\lambda t$ is large, say, $\lambda t>10$. Then $f(t)=C_{t}^{2}t^{2H}\simeq 2\Gamma(2H)(2\lambda)^{-2H}$. By final value theorem, we have $\lim\limits_{s\to 0}f(s)=2\Gamma(2H)(2\lambda)^{-2H}/s$; substituting it into (3.3) leads to

[TABLE]

with $A=m+2\xi\Gamma(2H)(2\lambda)^{-2H}$. Then

[TABLE]

Hence $\langle x(t)\rangle\simeq\frac{mv_{0}}{A}t$ for large $t$, and

[TABLE]

with $\langle F(t-t_{1})F(t-t_{2})\rangle=k_{B}T\xi h^{-2}(C_{t_{1}-t_{2}+h}^{2}|t_{1}-t_{2}+h|^{2H}+C_{t_{1}-t_{2}-h}^{2}|t_{1}-t_{2}-h|^{2H}-2C_{t_{1}-t_{2}}^{2}|t_{1}-t_{2}|^{2H})$. Let $g(t_{1})=\int_{0}^{t_{1}}t_{2}\langle F(t-t_{1})F(t-t_{2})\rangle\rmd t_{2}$. Using the property of convolution and Taylor’s series expansion, we have $g(s)=\mathcal{L}[g(t_{1})]\simeq k_{B}T\xi f(s)$. So $g(t_{1})\simeq k_{B}T\xi f(t_{1})$ for large $t_{1}$. Then

[TABLE]

for $t\rightarrow\infty$. Eq. (3.8) is confirmed by the numerical simulations; see Figure LABEL:fig:1. And the MSD of the system is

[TABLE]

displaying the ballistic diffusion.

Besides, let us further consider the autocorrelation function of the position $x(t)$. After making double Laplace transforms of the autocorrelation function of internal noise $\langle F(t_{1})F(t_{2})\rangle$, we have

[TABLE]

Then

[TABLE]

Taking inverse Laplace transforms of (3.11) and letting $t_{1},t_{2}$ tend to infinity, we obtain

[TABLE]

which results in (3.8) when $t_{1}=t_{2}$.

Furthermore, if $\lambda$ is moderate and $t\rightarrow\infty$, the velocity

[TABLE]

Then we have $\langle v(t)\rangle\simeq\frac{mv_{0}}{A}$ and the autocorrelation function of $v(t)$ is

[TABLE]

If $t_{1}=t_{2}=t$, then

[TABLE]

as $t\rightarrow\infty$. From (3.14), the second moment of position $x(t)$ is

[TABLE]

where $C_{t}^{2}t^{2H}\simeq 2\Gamma(2H)(2\lambda)^{-2H}$ for large $\lambda t$ is used; Eq. (3.8) is confirmed again. The correlation function of $x(t)$ and $v(t)$ is

[TABLE]

Remark: There is another way to obtain the asymptotic expression of $f(s)$ for small $s$. By using the formula [18]:

[TABLE]

where $\rho>0$, $\rm{Re}(s)>0$, $\rm{Re}(\alpha)>0$, $\rm{Re}(\mu/\rho)>|\rm{Re}(\nu)|$ and $K_{\nu}(\alpha t)$ is the modified Bessel function of the second kind, we obtain the Laplace transform of $f(t)\,(=C_{t}^{2}t^{2H})$ as

[TABLE]

So $f(s)\simeq\frac{2\Gamma(2H)(2\lambda)^{-2H}}{s}$ for small $s$.

Case II ($s$ is large, i.e., $t\rightarrow 0$, and $H\in(\frac{1}{2},1)$): For moderate or small $\lambda$, $\lambda t\rightarrow 0$. Then we have $C_{t}^{2}\simeq 2D_{H}\Gamma^{2}(H+1/2)$ and $f(t)\simeq 2D_{H}\Gamma^{2}(H+1/2)t^{2H}$ with its Laplace transform $f(s)\simeq 2D_{H}\Gamma^{2}(H+1/2)\Gamma(2H+1)s^{-1-2H}$. So (3.3) becomes

[TABLE]

for large $s$. Taking inverse Laplace transform on Eq. (3.18) leads to

[TABLE]

Then the mean $\langle x(t)\rangle\simeq v_{0}t$. Similarly to the case $t\rightarrow\infty$, we obtain

[TABLE]

for small $t$. Note that for short times we have $\langle x^{2}(t)\rangle\simeq\frac{k_{B}T}{m}t^{2}$ if the thermal initial condition $\langle v_{0}^{2}\rangle=\frac{k_{B}T}{m}$ is assumed. Figure LABEL:2 numerically confirms the theoretical result. Naturally, in this case, $v(t)\simeq\frac{1}{m}\int_{0}^{t}F(\tau)\rmd\tau+v_{0}$. So $\langle v(t)\rangle\simeq v_{0}$ and $\langle v^{2}(t)\rangle\simeq\frac{4k_{B}T\xi D_{H}\Gamma^{2}(H+\frac{1}{2})}{m^{2}}t^{2H}+v_{0}^{2}\simeq v_{0}^{2}$ for small $t$.

Case III ($\lambda$ is small, $t$ is large, and $H\in(\frac{1}{2},1)$): If $\lambda t\rightarrow\infty$, it reduces to Case I. Now we only consider the case that $\lambda t$ is small. As the second case, when $\lambda t$ is small, $C_{t}^{2}\simeq 2D_{H}\Gamma^{2}(H+1/2)$. Then (3.3) becomes

[TABLE]

for small $s$, which has the inverse Laplace transform

[TABLE]

with $P=\frac{1}{2\xi D_{H}\Gamma^{2}(H+\frac{1}{2})\Gamma(2H+1)\Gamma(2-2H)}$. So finally we obtain $\langle x(t)\rangle\simeq mv_{0}Pt^{1-2H}$ and the MSD is

[TABLE]

Similarly, $v(t)\simeq\frac{1}{m}\int_{0}^{t}F(\tau)E_{2H}(-B(t-\tau)^{2H})\rmd\tau+v_{0}E_{2H}(-Bt^{2H})$ by making the inverse Laplace transform of $v(s)\simeq\frac{F(s)+mv_{0}}{ms+2\xi D_{H}\Gamma^{2}(H+\frac{1}{2})\Gamma(2H+1)s^{1-2H}}$, where $B=2\xi D_{H}\Gamma^{2}(H+\frac{1}{2})\Gamma(2H+1)m^{-1}$ and $E_{\alpha}(z)$ is the particular case of the generalized Mittag-Leffler function [19]

[TABLE]

with $E_{\alpha}(z)=E_{\alpha,1}(z)$. The asymptotic expression of the generalized Mittag-Leffler function for $z\rightarrow\infty$ is $E_{\alpha,\beta}(-z)\simeq\frac{z^{-1}}{\Gamma(\beta-\alpha)}$, where $0<\alpha<1$ or $1<\alpha<2$.

It can be noted that $\langle F(t_{1})F(t_{2})\rangle\simeq 4k_{B}T\xi D_{H}\Gamma^{2}(H+1/2)H(2H-1)|t_{1}-t_{2}|^{2H-2}$ when $\lambda t$ is small, being similar to the autocorrelation function of fGn [10]; and the expression of $v(t)$ is similar to the one of fLe. So in this case, the asymptotic behavior of $\langle v(t)\rangle$ and $\langle v^{2}(t)\rangle$ is consistent with that of fLe. Thus $\langle v(t)\rangle\simeq v_{0}E_{2H}(-Bt^{2H})\simeq t^{-2H}$ and $\langle v^{2}(t)\rangle$ decays like $t^{-4H}$ for large $t$ [12].

From the above discussions, we know that the MSD of tfLe strongly depends on the value of $\lambda$ and $t$. When $\lambda t$ is small, it displays the same asymptotic behavior as the one of fLe [10]. As time goes on, for fixed $\lambda$ the MSD of tfLe transmits as $t^{2}\to t^{2-2H}\to t^{2}$; the smaller $\lambda$ is, the longer the time of the MSD behaves as $t^{2-2H}$. See Figure LABEL:3 for the demonstration of the results.

3.1.2 Overdamped generalized Langevin equation

The overdamped generalized Langevin equation without Newton’s acceleration term reads as

[TABLE]

where $K(t)$ and $F(t)$ are the same as the ones in Eq. (3.1). For moderate $\lambda$, taking the same procedure as above subsection, we can easily obtain that $\langle x(t)\rangle=0$ and

[TABLE]

for large $t$, being the same as the underdamped case (3.1). For short times,

[TABLE]

then we have $\langle x(t)\rangle=0$ and the MSD

[TABLE]

It can be noted that for short times, the MSD of the overdamped tfLe transits from $t^{2-2H}$ to $t^{2}$, and it behaves the same as the one of overdamped fLe [10] for small $\lambda t$. See Figure LABEL:3-1 for the numerical simulations, which verify (3.25) and (3.27).

3.2 Harmonic potential

Now we further consider the tfLe (3.1) with external potential $U(x)$, i.e.,

[TABLE]

where $-U^{\prime}(x)$ is an external force. If the external potential is a harmonic potential $U(x)=\frac{1}{2}m\omega^{2}x^{2}(t)$, where $\omega$ is the frequency of the oscillator, then we have

[TABLE]

In what follows we analyze the normalized displacement correlation function, which is defined by

[TABLE]

under the thermal initial conditions $\langle F(t)x(0)\rangle=0$, $\langle x^{2}(0)\rangle=\frac{k_{B}T}{m\omega^{2}}$, and $\langle x(0)v(0)\rangle=0$. Making Laplace transform of (3.29) leads to

[TABLE]

Then

[TABLE]

which results in

[TABLE]

with $I(s)=\frac{s^{-1}}{ms^{2}+\xi sK(s)+m\omega^{2}}$.

$K(s)$ and $f(t)$ are given below Eq. (3.2). Here we carefully take the asymptotic expression of $f(t)$ as

[TABLE]

Then

[TABLE]

where $a=m+2\xi\Gamma(2H)(2\lambda)^{-2H}$ and $b=2\xi\Gamma(H+\frac{1}{2})^{2}(2\lambda)^{-H-\frac{1}{2}}$. The approximation (3.34) is valid just for small $\omega$; for large $\omega$, the convergence region of the approximation of $I(s)$ may be different from the one of $I(s)$. One can note that the approximation in （(3.34) is not sensitive to the value of $H$, which is also illustrated by simulations (see Figs. LABEL:you and LABEL:0-8). We simply take $H=\frac{1}{2}$; for other values of $H$, one can use the techniques in [20, 21] to make inverse Laplace transform of $I(s)$, but the calculations are very complicated. Then

[TABLE]

Rewrite (3.35) as $a\lambda I(s)=\frac{1}{P(s)}+\frac{s^{-1}\lambda}{P(s)}$ with $P(s)=s^{2}+\frac{m\omega^{2}}{a\lambda}s+\frac{m\omega^{2}}{a}$. Then the solutions of $P(s)=0$ are $z_{1}=\alpha+\rmi\beta$ and $z_{2}=\alpha-\rmi\beta$, where $\alpha=-\frac{m\omega^{2}}{2a\lambda}$ and $\beta=\frac{\omega\sqrt{4ma\lambda^{2}-m^{2}\omega^{2}}}{2a\lambda}$. The roots are complex on account of the small $\omega$ satisfing $m\omega^{2}<4a\lambda^{2}$. Denoting

[TABLE]

we have

[TABLE]

and

[TABLE]

where $A_{1}=\frac{1}{z_{1}-z_{2}}$ and $A_{2}=\frac{1}{z_{2}-z_{1}}$.

Plugging the concrete expressions of $A_{k}$ into (3.37) leads to

[TABLE]

and the normalized displacement correlation function is

[TABLE]

where $\theta=\arctan(-\frac{\sqrt{4ma\lambda^{2}-m^{2}\omega^{2}}}{m\omega})$. It shows that for small $\omega$, the phase, amplitude, and period are $\theta$, $2\lambda\sqrt{\frac{a}{4a\lambda^{2}-m\omega^{2}}}\rme^{-\frac{m\omega^{2}}{2a\lambda}t}$, and $\frac{4\pi a\lambda}{\omega\sqrt{4ma\lambda^{2}-m^{2}\omega^{2}}}$, respectively, being verified by the simulations given in Figure LABEL:you (a) and Figure LABEL:you (b). The simulation results for $H=0.8$ are presented in Figure LABEL:0-8, in which (a) shows that after changing the value of $H$, besides a little bit of the difference of the amplitude, the simulation results are still consistent with Eq. (3.39). We numerically detect that $C_{x}(t)$ is always zero crossing. For $\lambda\rightarrow 0$, the results for $C_{x}(t)$ of the fLe [20] are recovered; see Figure LABEL:wu.

4 Conclusion

Tempered fractional Brownian motion is a recently introduced stochastic process, displaying localization diffusion. We further discuss its ergodicity and derive its Fokker-Planck equation. Then we introduce the tempered fractional Langevin equation with the tempered fractional Gaussian noise as internal noise. Both the undamped and overdamped cases are considered, and the evolution of the MSD is carefully investigated, revealing the procedure of transition and finally turning to the ballistic diffusion for long time. The normalized displacement correlation function $C_{x}(t)$ of the tempered fractional Langevin equation with harmonic potential is explicitly derived for small frequency. By the algorithm provided in C, almost all of the theoretical results are verified by numerical simulations.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under Grant No. 11671182.

Appendix A Derivation of (2.7) and (2.8)

For tfBm, denoting $g(z)=C_{z}^{2}|z|^{2H}$ with $C_{z}^{2}$ defined by (2.1), the average of $\bar{\delta}^{2}(x(t))$ is

[TABLE]

and the variance of $\bar{\delta}^{2}(x(t))$ is

[TABLE]

for long times and moderate $\lambda$, where

[TABLE]

[TABLE]

[TABLE]

Appendix B Asymptotic behavior of the covariance function of tfGn

Fixing the value of $\lambda$ and considering sufficiently small $t$, we have

[TABLE]

where $C_{t}^{2}$ is a constant. Figure LABEL:5 is plotting the asymptotic second integration of $K(t)$, i.e., $C_{t}^{2}|t|^{2H}$. The convexity (concavity) of $C_{t}^{2}|t|^{2H}$ implies the positivity (negativity) of $K(t)$. It is clear that $K(t)>0$ (see Figure LABEL:5(a)) for $\frac{1}{2}<H<1$, and $K(t)<0$ (see Figure LABEL:5(b)) for $0<H<\frac{1}{2}$.

Fixing $\lambda$ and letting $t\rightarrow+\infty$, i.e., $\lambda t$ is large, we have

[TABLE]

Using Taylor’s series expansion, along with $C_{t}^{2}=\frac{2\Gamma(2H)}{(2\lambda|t|)^{2H}}-\frac{2\Gamma(H+\frac{1}{2})K_{H}(\lambda|t|)}{\sqrt{\pi}(2\lambda|t|)^{H}}=\bar{A}|t|^{-2H}-\bar{B}|t|^{-H}K_{H}(\lambda|t|)$ and $K_{H}(\lambda t)\simeq\sqrt{\pi}(2\lambda t)^{-\frac{1}{2}}\rme^{-\lambda t}$ as $\lambda t$ is large, the first term in (B) is $2C_{t}^{2}H(2H-1)t^{2H-2}\simeq 2\bar{A}H(2H-1)t^{-2}-2\bar{B}H(2H-1)\sqrt{\pi}(2\lambda)^{-\frac{1}{2}}t^{H-\frac{5}{2}}\rme^{-\lambda t},$ and the second term in (B) is

[TABLE]

where $g(x)=x^{H-\frac{1}{2}}\rme^{-\lambda x}$. Combining the above two estimations leads to

[TABLE]

for large $t$; see Figure LABEL:5. The $K(t)$ can also be obtained by making second derivative on the asymptotic expression of $C_{t}^{2}t^{2H}$ for large $t$. So the asymptotic behavior of tfGn’s covariance function is $\langle\gamma(0)\gamma(t)\rangle\simeq-\frac{\Gamma(H+\frac{1}{2})\lambda^{2}}{(2\lambda)^{H+\frac{1}{2}}}t^{H-\frac{1}{2}}\rme^{-\lambda t}$ for large $t$ and $\langle\gamma(0)\gamma(t)\rangle\simeq C_{t}^{2}H(2H-1)t^{2H-2}$ with a constant $C_{t}^{2}$ for small $t$.

Appendix C Algorithm for numerical simulations

To generate tempered fractional Gaussian noise $X_{0},X_{1},\cdots$, we adopt the Hosking method [22], which works for the general stationary Gaussian process. The key observation of this algorithm is to generate $X_{n+1}$ by the conditional distribution of $X_{n+1}$ given $X_{n},\cdots,X_{0}$ recursively. The covariance function of tfGn is

[TABLE]

for $n,k=0,1,2,\cdots$. Note that $\langle B_{\alpha,\lambda}^{2}(t)\rangle=C_{t}^{2}|t|^{2H}$ from (2.2) and $B_{\alpha,\lambda}(t)$ has the stationary increments, $X_{n}:=B_{\alpha,\lambda}(n+1)-B_{\alpha,\lambda}(n)\sim N(0,C_{1}^{2})$, which means $\rho(0)=\mathbb{E}(X_{n}^{2})=C_{1}^{2}$. Furthermore, let $D(n)=\left(\rho(i-j)\right)_{i,j=0,\cdots,n}$ be the covariance matrix and $c(n)$ be the $(n+1)-$column vector with elements $c(n)_{k}=\rho(k+1),k=0,\cdots,n$. Define the $(n+1)\times(n+1)$ flipping matrix $F(n)=(\mathbf{1}(i=n-j))_{i,j=0,\cdots,n}$, where $\mathbf{1}$ denotes the indicator function.

We claim that the conditional distribution of $X_{n+1}$ is Gaussian with expectation $\mu_{n}$ and variance $\sigma_{n}^{2}$ given by

[TABLE]

To avoid matrix inversion in (3.1) in each step, define $d(n):=D(n)^{-1}c(n)$, and $\tau_{n}:=d(n)^{\prime}F(n)c(n)=c(n)^{\prime}F(n)d(n)$. Split the matrix $D(n+1)$ as follows:

[TABLE]

With some simple calculations, one gets

[TABLE]

From (3.6), for each $x\in\mathbb{R}^{n+1}$ and $y\in\mathbb{R}$, we have

[TABLE]

This implies that the conditional distribution of $X_{n+1}$ is indeed Gaussian with expectation $\mu_{n}$ and variance $\sigma_{n}^{2}$. On the other hand, by (3.9), some recursions are as follows:

[TABLE]

and

[TABLE]

with

[TABLE]

We start the recursion with $\mu_{0}=\rho(1)X_{0},~{}\sigma_{0}^{2}=C_{1}^{2}-\rho(1)^{2}$ and $\tau_{0}=\frac{1}{C_{1}^{2}}\rho(1)^{2}$. Taking cumulative sums on the generated tfGn samples $X_{0},\cdots,X_{n}$, one obtains the tempered fractional Brownian motion sample $B_{\alpha,\lambda}(k),k=0,\cdots,n$.

Next, we consider to generate a series tfGn samples of number $N$ for obtaining $B_{\alpha,\lambda}(T)$. The scaling property of tfBm shows that $B_{\alpha,\lambda}(T)=\left(\frac{T}{N}\right)^{H}B_{\alpha,\lambda T/N}(N)$ [15]. Denoting $\gamma(T):=\gamma_{\alpha,\lambda}(T)=\left(\frac{T}{N}\right)^{H-1}\gamma_{\alpha,\lambda T/N}(N)$, i.e., $h=\frac{T}{N}$ in (2.9), for simulating the second moment of position $\langle x^{2}(t)\rangle$, we solve (3.1) with the scheme:

[TABLE]

where $h$ is the step length, and $K(t)=\frac{1}{h^{2}}(C_{t+h}^{2}|t+h|^{2H}+C_{t-h}^{2}|t-h|^{2H}-2C_{t}^{2}|t|^{2H})$. In simulating the normalized displacement correlation function $C_{x}(t)$, we take the initial conditions as $v_{0}\sim N(0,\frac{k_{B}T}{m})$, and $x_{0}\sim U[-r,r]$, where $r=\sqrt{\frac{3k_{B}T}{m}}\frac{1}{\omega}$.