Subdiffusion in an external force field

TL;DR

This paper investigates how external forces influence subdiffusive processes modeled by subordinated Langevin equations, analyzing statistical properties and deriving Fokker-Planck equations for non-stationary, non-ergodic, aging systems.

Contribution

It introduces a comprehensive analysis of external force effects on subdiffusion using subordinated Langevin equations, including statistical measures and Fokker-Planck equations, highlighting different force interactions.

Findings

External forces affect subdiffusion differently depending on their timing and nature.

Statistical quantities reveal non-stationarity, non-ergodicity, and aging in the processes.

Derived Fokker-Planck equations describe the dynamics under various external forces.

Abstract

The phenomena of subdiffusion are widely observed in physical and biological systems. To investigate the effects of external potentials, say, harmonic potential, linear potential, and time dependent force, we study the subdiffusion described by subordinated Langevin equation with white Gaussian noise, or equivalently, by the single Langevin equation with compound noise. If the force acts on the subordinated process, it keeps working all the time; otherwise, the force just exerts an influence on the system at the moments of jump. Some common statistical quantities, such as, the ensemble and time averaged mean squared displacement, position autocorrelation function, correlation coefficient, generalized Einstein relation, are discussed to distinguish the effects of various forces and different patterns of acting. The corresponding Fokker-Planck equations are also presented. All the stochastic processes discussed here are non-stationary, non-ergodicity, and aging.

Equations

$$⟨(\Delta y(t){)}^2⟩=⟨[y(t)-⟨y(t)⟩{]}^2⟩\simeq t^{\beta }(\beta \ne 1),$$

$$\dot{x}(s)=\sqrt{2\sigma }\xi (s),\dot{t}(s)=\eta (s),$$

$$\dot{y}(t)=\sqrt{2\sigma }\overline{\xi }(t).$$

$$s(t)=\underset{s>0}{\inf }\{s:t(s)>t\}.$$

$${\mathcal{L}}_{t\to \lambda }[h(s,t)]={\int }_0^{\infty }{e}^{-\lambda t}h(s,t)dt={\lambda }^{\alpha -1}{e}^{-s{\lambda }^{\alpha }}.$$

$$p(y,t)={\int }_0^{\infty }{p}_0(y,s)h(s,t)ds,$$

$${\mathcal{L}}_{t\to \lambda }⟨y^n(t)⟩={\lambda }^{\alpha -1}{\mathcal{L}}_{s\to {\lambda }^{\alpha }}⟨x^n(s)⟩$$

$$\begin{array}{rl}& p(y_2,t_2;y_1,t_1)\\ & ={\int }_0^{\infty }{\int }_0^{\infty }{p}_0(y_2,s_2;y_1,s_1)h(s_2,t_2;s_1,t_1)d{s}_{1}d{s}_2,\end{array}$$

$$\begin{array}{rl}& ⟨y({\lambda }_1)y({\lambda }_2)⟩\\ & ={\int }_0^{\infty }{\int }_0^{\infty }⟨x(s_1)x(s_2)⟩h(s_2,{\lambda }_2;s_1,{\lambda }_1)d{s}_{1}d{s}_2,\end{array}$$

$$\begin{array}{rl}& h(s_2,{\lambda }_2;s_1,{\lambda }_1)\\ & =\delta (s_2-s_1)\frac{{\lambda }_1^{\alpha }-({\lambda }_1+{\lambda }_2{)}^{\alpha }+{\lambda }_2^{\alpha }}{{\lambda }_1{\lambda }_2}{e}^{-s_1({\lambda }_1+{\lambda }_2{)}^{\alpha }}\\ & +\Theta (s_2-s_1)\frac{{\lambda }_2^{\alpha }[({\lambda }_1+{\lambda }_2{)}^{\alpha }-{\lambda }_2^{\alpha }]}{{\lambda }_1{\lambda }_2}\\ & \cdot e^{-({\lambda }_1+{\lambda }_2{)}^{\alpha }{s}_1}{e}^{-{\lambda }_2^{\alpha }(s_2-s_1)}+\Theta (s_1-s_2)\\ & \cdot \frac{{\lambda }_1^{\alpha }[({\lambda }_1+{\lambda }_2{)}^{\alpha }-{\lambda }_1^{\alpha }]}{{\lambda }_1{\lambda }_2}{e}^{-({\lambda }_1+{\lambda }_2{)}^{\alpha }{s}_2}{e}^{-{\lambda }_1^{\alpha }(s_1-s_2)}.\end{array}$$

$$⟨y^2(t)⟩=\frac{2\sigma }{\Gamma (1+\alpha )}{t}^{\alpha }$$

$$\overline{{\delta }^2(\Delta )}\simeq \frac{2\sigma }{\Gamma (1+\alpha )}\Delta T^{\alpha -1}$$

$$\frac{\partial p(y,t)}{\partial t}=\sigma \frac{{\partial }^2}{\partial y^2}{D}_t^{1-\alpha }p(y,t).$$

$$D_t^{1-\alpha }p(y,t)=\frac{1}{\Gamma (\alpha )}\frac{\partial }{\partial t}{\int }_0^t(t-t^{\prime }{)}^{\alpha -1}p(y,t^{\prime })d{t}^{\prime },$$

$$\dot{x}(s)=-\gamma x(s)+\sqrt{2\sigma }\xi (s),\dot{t}(s)=\eta (s),$$

$$\dot{y}(t)=-\gamma y(t)\dot{s}(t)+\sqrt{2\sigma }\overline{\xi }(t),$$

$$\begin{array}{rl}y(t)& =x(s(t))\\ & =-\gamma {\int }_0^{s(t)}x(s^{\prime })d{s}^{\prime }+\sqrt{2\sigma }{\int }_0^{s(t)}\xi (s^{\prime })d{s}^{\prime }\\ & =-\gamma {\int }_0^{t}x(s(\tau ))ds(\tau )+\sqrt{2\sigma }{\int }_0^t\xi (s(\tau ))ds(\tau )\\ & =-\gamma {\int }_0^{t}y(\tau )ds(\tau )+\sqrt{2\sigma }{\int }_0^t\xi (s(\tau ))ds(\tau ).\end{array}$$

$$⟨y^2(t)⟩=\frac{\sigma }{\gamma }-\frac{\sigma }{\gamma }{E}_{\alpha }(-2\gamma t^{\alpha }),$$

$$⟨y^2(t)⟩\simeq \frac{2\sigma }{\Gamma (1+\alpha )}{t}^{\alpha },$$

$$⟨y^2(t)⟩\simeq \frac{\sigma }{\gamma }-\frac{\sigma }{2{\gamma }^2\Gamma (1-\alpha )}{t}^{-\alpha }$$

$$⟨y(t_1)y(t_2)⟩\simeq \frac{\sigma }{\gamma }\frac{B(t_1/t_2,\alpha ,1-\alpha )}{\Gamma (\alpha )\Gamma (1-\alpha )}$$

$$\begin{array}{rl}r[y(t_1),y(t_2)]& =\frac{⟨(y(t_1)-⟨y(t_1)⟩)(y(t_2)-⟨y(t_2)⟩)⟩}{\sqrt{⟨(y(t_1)-⟨y(t_1)⟩{)}^2⟩⟨(y(t_2)-⟨y(t_2)⟩{)}^2⟩}}\\ & \simeq \frac{1}{\Gamma (1-\alpha )\Gamma (1+\alpha )}{\left(\frac{t_1}{t_2}\right)}^{\alpha },\end{array}$$

$$r_0[y(t_1),y(t_2)]={\left(\frac{t_1}{t_2}\right)}^{\frac{\alpha }{2}}.$$

$$\frac{\partial p(y,t)}{\partial t}={\mathcal{L}}_{\text{FP}}{D}_t^{1-\alpha }p(y,t)$$

$$\begin{array}{rl}\overline{{\delta }^2(\Delta )}& =\frac{1}{T-\Delta }{\int }_0^{T-\Delta }[(y(t+\Delta )-y(t))\\ & -⟨y(t+\Delta )-y(t)⟩{]}^{2}dt,\end{array}$$

$$⟨\overline{{\delta }^2(\Delta )}⟩\simeq \frac{2\sigma }{\gamma }\frac{\text{sin}(\alpha \pi )}{\alpha (1-\alpha )\pi }{\left(\frac{\Delta }{T}\right)}^{1-\alpha }$$

$$\mathcal{E}\mathcal{B}\simeq {\left(\frac{T}{\Delta }\right)}^{\alpha -1},$$

$$\mathcal{E}\mathcal{B}\simeq \frac{2\sin (\alpha \pi )}{\alpha (1-\alpha )\pi }{\left(\frac{T}{\Delta }\right)}^{\alpha -1}.$$

$$\dot{y}(t)=-\gamma y(t)+\sqrt{2\sigma }\overline{\xi }(t),$$

$$\begin{array}{rl}⟨\overline{\xi }(t_1)\overline{\xi }(t_2)⟩& ={\mathcal{L}}_{{\lambda }_1\to t_1,{\lambda }_2\to t_2}^{-1}[({\lambda }_1+{\lambda }_2{)}^{-\alpha }]\\ & =t_1^{\alpha -1}\delta (t_1-t_2)/\Gamma (\alpha ).\end{array}$$

Full text

I Introduction

In the natural world, it is hard to find the real free particles; actually almost all the time, they are in some kinds of external potentials. The motion of particles in complex disordered systems generally is no longer Brownian, exhibiting anomalous diffusion behavior Jeon and Metzler (2012); Eule and Friedrich (2009); Cairoli and Baule (2015a); Magdziarz et al. (2008); Chen et al. (2017); Fedotov and Korabel (2015), which is characterized by the nonlinear evolution in time of the mean squared displacement (MSD) of particles; i.e.,

[TABLE]

which represents subdiffusion for $0<\beta<1$ and superdiffusion for $\beta>1$; for the case $\beta=2$, it is called ballistic diffusion and $\beta=0$ the localization diffusion Metzler and Klafter (2000); Bouchaud and Georges (1990).

One of the most powerful and popular models to describe anomalous diffusion is continuous-time random walk (CTRW), which was originally introduced by Montroll and Weiss in 1965 Montroll and Weiss (1965), extending regular random walks on lattices to a continuous-time variable. It has been successfully applied in various fields, such as, the charge carrier transport in amorphous semiconductors Scher and Montroll (1975), electron transfer Nelson (1999), dispersion in turbulent systems Solomon et al. (1993), and so on.

Another special model to describe the complex dynamics is Langevin equation; its classical version should be the differential equation form of Newton’s second law. Compared with CTRW, the Langevin picture has a striking advantage in characterizing external fields. Of course, it also has a close connection with CTRW model. In 1994, Fogedby Fogedby (1994) used the stochastic time changed method to introduce an equivalent form of the continuum limit of the subdiffusive CTRW — a Langevin equation coupled with a subordinator, i.e.,

[TABLE]

where $x(s)$ is named as original process with respect to internal time $s$, $\xi(s)$ is a white Gaussian noise with null mean value and autocorrelation function $\langle\xi(s)\xi(s^{\prime})\rangle=\delta(s-s^{\prime})$, and $\eta(s)$ is a fully skewed $\alpha$-stable Lévy noise with $0<\alpha<1$ Schertzer et al. (2001) and usually regarded as the formal derivative of the $\alpha$-stable subordinator $t(s)$ Applebaum (2009). The time changed process $y(t):=x(s(t))$, where $s(t)$ is the inverse $\alpha$-stable subordinator Kumar and Vellaisamy (2015); Alrawashdeh et al. (2017), is an equivalent stochastic trajectory to the continuum limit of the CTRW with power-law distributed waiting times. Since then, the subordination Applebaum (2009), which was put forward by Bochner Bochner (1949) in 1949, has become a useful method to describe the time-changed stochastic processes exhibiting anomalous diffusion. Especially in recent years, the coupled Langevin equations have been widely investigated Meerschaert and Scheffler (2004); Gajda and Wyłomańska (2015); Wyłomańska et al. (2016); Chen et al. (2018), and the time-changed stochastic processes are important models in many fields, such as, biology Golding and Cox (2006), physics Nezhadhaghighi et al. (2011), ecology Scher et al. (2002), etc.

In fact, the coupled Langevin equation (2) describing the subdiffusion dynamics can also be rewritten into a single Langevin equation in physical time $t$ with an additive compound noise $\overline{\xi}(t)=\int_{0}^{+\infty}\xi(\tau)\delta(t-t(\tau))d\tau$ Cairoli and Baule (2015b):

[TABLE]

Besides the discussions on the models (2) or (3), there are some research works for them with external forces Fedotov and Korabel (2015); Eule et al. (2007); Burov et al. (2010); Cairoli and Baule (2015b, a); Dieterich et al. (2015); Magdziarz et al. (2008); Weron and Magdziarz (2008); Chen et al. (2018), which mainly presented the asymptotic expression of the MSD of the stochastic process for long times, depending on the one-point probability density function (PDF) of the stochastic process. The studies on more general statistical quantities, such as, the correlation coefficient, which reflects the correlation of positions at two different times, as well as the time averaged MSD, are sometimes ignored. These statistical quantities are significantly important to distinguish the processes with the same diffusion behavior.

In this paper, we investigate the influence of three kinds of common external forces — position-dependent force, constant force, and time-dependent force. These forces may act on the original process $x(s)$ in (2) or on the subordinated process $y(t)$ in (3) for different physical realities. These two acting patterns are, respectively, for the cases, where the external force only modifies the dynamical behavior at the moments of jump or exerts effect for the whole time. The comparisons are made for various effects exerted by different acting patterns with different external forces through some common statistical quantities, such as, ensemble and time averaged MSD, correlation coefficient, and ergodicity breaking parameter. These quantities mainly depend on the two-point joint PDF of the observed processes, except the ensemble averaged MSD. For different patterns of the force acting on the Langevin equation (2) or (3), the methods of obtaining the position autocorrelation function are different, which are fully demonstrated in this paper.

One interesting finding is that the position-independent force acting on the subordinated process in (3) does not change the diffusion behavior, ergodic property, and the correlation coefficient, while the position-dependent external force does. But if acting on the original process in (2), the external forces (position-dependent or position-independent) produce different results for almost all the statistical quantities, compared with the ones of free particles. In addition, the exponent of ergodicity breaking parameter does not depend on the forces and the acting patterns.

The Fokker-Planck equations govern the PDF $p(y,t)$; generally, they vary with the change of the processes described by the Langevin equations (2) and (3) with forces. For the equations, it is found that the Riemann-Liouville fractional derivative with respect to time is included when the external force affects the process only at the moments of jump, while the fractional substantial derivative and another kind of novel fractional derivative are needed when the external force acts on the system for the whole time.

The structure of this paper is as follows. In Sec. II, we review the subordinator as well as the inverse subordinator, and briefly present the method of subordination we mainly use. Then we consider the effects of position-dependent force (harmonic potential), constant force (linear potential), and time-dependent force, respectively, in Sec. III–V. The potential properties of the diffusion behaviors are revealed through various statistical quantities. Finally, we make the summaries in Sec. VI and the detailed derivations of some of the results of the paper are presented in Appendix.

II Subordinator

Subordinator is a non-decreasing Lévy process with stationary and independent increments Applebaum (2009) and it can be regarded as a stochastic model of time evolution. The subordinator $t(s)$ in this paper is taken to be $\alpha$-stable one with $0<\alpha<1$ Applebaum (2009), which has the characteristic function $\langle\textrm{e}^{-\lambda t(s)}\rangle=\textrm{e}^{-s\lambda^{\alpha}}$. The brackets $\langle\cdots\rangle$ denote the statistical average over stochastic realizations. The corresponding inverse process, called inverse $\alpha$-stable subordinator $s(t)$ Kumar and Vellaisamy (2015); Alrawashdeh et al. (2017), is the first-passage time of the subordinator $\{t(s),\,s\geq 0\}$ , defined as

[TABLE]

In addition, we denote the PDF of the inverse $\alpha$-stable subordinator $s(t)$ as $h(s,t)$ and its Laplace transform ($t\rightarrow\lambda$) $h(s,\lambda)$ is Baule and Friedrich (2005)

[TABLE]

The PDF $p(y,t)$ of the subordinated process $y(t):=x(s(t))$ can be written as Baule and Friedrich (2005); Barkai (2001); Chen et al. (2018)

[TABLE]

where $p_{0}(x,s)$ is the PDF of the original process $x(s)$. The moments of the subordinated process $y(t)$ can be obtained through the relation

[TABLE]

in Laplace space. Similarly, the two-point joint PDF $p(y_{2},t_{2};y_{1},t_{1})$ of $y(t)$ can be obtained through the two-point joint PDF $p_{0}(x_{2},s_{2};x_{1},s_{1})$ of the original stochastic process $x(s)$,

[TABLE]

where $h(s_{2},t_{2};s_{1},t_{1})$ is the two-point joint PDF of the inverse subordinator $s(t)$. The correlation function of $y(t)$ in Laplace space ($t_{1}\rightarrow\lambda_{1}$, $t_{2}\rightarrow\lambda_{2}$) is

[TABLE]

with Baule and Friedrich (2005)

[TABLE]

Based on the formulae above, we have the MSD of the stochastic process $y(t)$ in (2) Metzler and Klafter (2000)

[TABLE]

and the autocorrelation function Baule and Friedrich (2005) $\langle y(t_{1})y(t_{2})\rangle=\frac{2\sigma}{\Gamma(1+\alpha)}t_{1}^{\alpha}$ for $t_{1}\leq t_{2}$. In addition, the time averaged MSD is He et al. (2008); Lubelski et al. (2008)

[TABLE]

for $\Delta\ll T$. The corresponding Fokker-Planck equation, governing the PDF $p(y,t)$ of finding the particle at position $y$ at time $t$, is Metzler and Klafter (2000, 2004)

[TABLE]

The symbol $D_{t}^{1-\alpha}$ is the Riemann-Liouville fractional derivative Podlubny (1999), defined as

[TABLE]

which is a nonlocal time derivative and indicates the non-Markovian property of the process $y(t)$.

III Subdiffusive dynamics in harmonic potential

In the following two subsections, we respectively discuss two cases: acting on the Langevin equations (2) and (3) by the harmonic potential. By comparing some statistical quantities, including ensemble and time averaged MSD, correlation coefficient, and ergodicity breaking parameter, we find some significant differences and interesting phenomena, especially in the latter case where the position-dependent external force acts on the system all the time.

III.1 Force acting on original process $x(s)$

Consider the Langevin system with a harmonic potential on the original process $x(s)$ Eule et al. (2007); Burov et al. (2010)

[TABLE]

where $\gamma$ is a positive constant, $\xi(s)$ and $\eta(s)$ are two independent noises defined in (2). The harmonic potential $V(x)=\gamma x^{2}/2$ leads to a friction-like force $F(x)=-dV(x)/dx=-\gamma x$ in the first equation of (15). Based on (15), a new single Langevin equation in physical time $t$ of the subordinated process $y(t)=x(s(t))$ can be obtained as

[TABLE]

with $\overline{\xi}(t)=\int_{0}^{+\infty}\xi(\tau)\delta[t-t(\tau)]d\tau$, or equivalently, $\overline{\xi}(t)=\xi(s(t))\dot{s}(t)$, since

[TABLE]

The noise $\overline{\xi}(t)$ here can be regarded as the formal derivative of the time-changed Brownian motion $B(s(t))$. The external force in (15) only changes the motion of the particles at the instant of jumps; in fact, this mechanism can be easily found from the equivalent Langevin equation in physical time (16), i.e., when a particle suffers a trapping event before next jump, the internal time process $s(t)$ remains a constant and the external force becomes zero due to $\dot{s}(t)=0$ in (16).

Using formula (7), it can be got that the first moment of the stochastic process $y(t)$ is zero due to symmetry and the MSD is

[TABLE]

by utilizing $\langle x^{2}(s)\rangle=\frac{\sigma}{\gamma}(1-e^{-2\gamma s})$ with the initial position $x_{0}=0$. Considering the asymptotic expression of the Mittag-Leffler function Erdelyi (1981) for small $t$: $E_{\alpha}(-2\gamma t^{\alpha})\simeq 1-\frac{2\gamma t^{\alpha}}{\Gamma(1+\alpha)}$, the asymptotic form of the MSD for short times $t\ll(2\gamma)^{-\frac{1}{\alpha}}$ is

[TABLE]

which coincides with the MSD of a free particle in (11) and implies that the harmonic potential does not affect the diffusion dynamics in short times. But for long times $t\gg(2\gamma)^{-\frac{1}{\alpha}}$, one has

[TABLE]

because of the asymptotic expression $E_{\alpha}(-2\gamma t^{\alpha})\simeq\frac{t^{-\alpha}}{2\gamma\Gamma(1-\alpha)}$ for large $t$. The saturation plateau value, denoted as $\langle y^{2}\rangle_{\textrm{th}}=\frac{\sigma}{\gamma}$, is approached at the power-law rate. The simulation results for different $\alpha$ are shown in Fig. 1. It can be seen that the MSD with a smaller $\alpha$ tends to the saturation plateau value more slowly, being an expected dynamical behavior within a confined harmonic potential due to smaller $\alpha$ corresponding to longer waiting time. This process behaves as a localization diffusion for long times. Compared with the original process $x(s)$, the MSD of which relaxes to the value $\frac{\sigma}{\gamma}$ exponentially, the subordinator $s(t)$ in this model only changes the convergence rate but keeps the same saturation plateau value.

In addition, the position autocorrelation function is Burov et al. (2010)

[TABLE]

for large $t_{1}$, $t_{2}$ ($t_{2}\geqslant t_{1}$), where $B(z,a,b)$ is the incomplete Beta function Abramowitz and Stegun (1972). Therefore, this process is non-stationary since the position autocorrelation function cannot be expressed as a function of time difference $|t_{1}-t_{2}|$. For fixed $t_{1}$ and $t_{2}\rightarrow\infty$, the correlation coefficient $r[y(t_{1}),y(t_{2})]$ of the stochastic process $y(t)$, which characterizes the correlation of position at two different times, can be obtained as

[TABLE]

where we have used the asymptotic expression of the incomplete Beta function for small $z$, i.e., $B(z,a,b)\simeq z^{a}/a$. It can be seen that the process $y(t)$ described by model (15) is long-range dependent due to the power-law decay with $0<\alpha<1$ of the correlation coefficient in (22). On the other hand, as a result of the harmonic potential, the correlation becomes weaker than that of a free particle, the correlation coefficient of which is

[TABLE]

The Fokker-Planck equation corresponding to the Langevin equation (15) is Metzler et al. (1999); Metzler and Klafter (2000); Gajda and Magdziarz (2010),

[TABLE]

with the Fokker-Planck operator $\mathcal{L}_{\textrm{{FP}}}=-\frac{\partial}{\partial y}F(y)+\sigma\frac{\partial^{2}}{\partial y^{2}}$ (here $F(y)=-\gamma y$), which can be derived by three different methods. The first one is based on the relation (6) between the PDF of subordinated process and original process Chen et al. (2018). The second one is to take the parameter $p=0$ in the Feynman-Kac equation Cairoli and Baule (2017); Wang et al. (2018). As for the last method, (24) can be got from the master equation in CTRW model within a harmonic potential Metzler and Klafter (2000); Carmi and Barkai (2011).

Let us turn to the time averaged MSD, defined as Metzler et al. (2014); Burov et al. (2011); Akimoto et al. (2018); Hou et al. (2018)

[TABLE]

where $\Delta$ is the lag time, and $T$ is the measurement time. We emphasize that the lag time $\Delta$ separating the displacement between trajectory points is much shorter than the measurement time $T$. Sometimes, the time averaged MSD of some stochastic systems confined in a harmonic potential approaches twice the ensemble averaged MSD $\langle y^{2}\rangle_{\textrm{th}}$ for long times, such as, overdamped Brownian motion, fractional Brownian motion, and fractional Langevin dynamics Jeon and Metzler (2012); Jeon et al. (2013). Different from it, the time averaged MSD of the confined model (15) is sublinear in lag time $\Delta$ Burov et al. (2010, 2011); Neusius et al. (2009)

[TABLE]

with $T\gg\Delta\gg(1/\gamma_{1})^{1/\alpha}$. Here $\gamma_{1}$ is the smallest nonzero eigenvalue of the Fokker-Planck operator $\mathcal{L}_{\textrm{{FP}}}$. The disagreement between the ensemble and time averaged MSD, the former is constant $\Delta^{0}$ and the latter scales as $\Delta^{1-\alpha}$, indicates non-ergodicity of the stochastic process $y(t)$. At short lag times $\Delta$, the linear scaling in lag times is observed Burov et al. (2010) $\langle\overline{\delta^{2}(\Delta)}\rangle\simeq\frac{2\sigma}{\Gamma(1+\alpha)}\Delta T^{\alpha-1}$, which is the same as the one of a free particle in (12) since the potential has not begun to affect the stochastic process.

Now, we consider the influence of the harmonic potential in terms of the ergodicity breaking parameter, defined as $\mathcal{EB}=\langle\overline{\delta^{2}}(\Delta)\rangle/\langle x^{2}(\Delta)\rangle.$ For the process of free particle, the ergodicity breaking parameter is

[TABLE]

while for the case with harmonic potential, the ergodicity breaking parameter becomes

[TABLE]

They have the same exponents of $T/\Delta$, but the coefficient of the latter is larger.

III.2 Force acting on subordinated process $y(t)$

The external force in (15) only makes an influence on the dynamical behavior at the moments of jump; contrary to it, the force may keep acting on the system all the time, even when the particle is trapped. Recently, such a model has been proposed in Cairoli and Baule (2015b), where the harmonic potential is assumed to directly act on the subordinated process $y(t)$ in physical times,

[TABLE]

where $\overline{\xi}(t)$ is the same noise as the one in (16). The two point correlation function of this compound noise $\overline{\xi}(t)$ Cairoli and Baule (2015b) could be gotten through the inverse Laplace transform

[TABLE]

The harmonic potential acts as a friction-like force $-\gamma y(t)$; even in the constant period of inverse subordinator $s(t)$, it still influences the stochastic dynamics. Actually, the Langevin equation (29) can be rewritten as a coupled Langevin system with a subordinator as

[TABLE]

More precisely, the solution of (29) is

[TABLE]

with initial condition $y_{0}=0$, which is equivalent to

[TABLE]

by replacing $s$ with $s(t)$. Compared with (15), the friction term $-\gamma x(s)$ is multiplied by the Lévy noise $\eta(s)$, which acts as a multiplicative noise in the first equation in (31).

From (32), it can be calculated that the mean of $y(t)$ is zero and the MSD is

[TABLE]

with the confluent hypergeometric function Abramowitz and Stegun (1972) ${}_{1}F_{1}(a,b;z)=\frac{\Gamma(b)}{\Gamma(a)\Gamma(b-a)}\int_{0}^{1}e^{zu}u^{a-1}(1-u)^{b-a-1}du$. The asymptotic expansion of MSD for short times $t\ll(2\gamma)^{-1}$ is

[TABLE]

which coincides to that of the free particle (11). For long times $t\gg(2\gamma)^{-1}$, using the asymptotic expansion ${}_{1}F_{1}(a,b;z)\simeq\Gamma(b)\left(e^{z}z^{a-b}/\Gamma(a)+(-z)^{-a}/\Gamma(b-a)\right)$ for large $z$ Abramowitz and Stegun (1972), we get

[TABLE]

which tends to zero at the power-law rate. The consistency between simulation and the theoretical results about the MSD of model (29) can be found in Fig. 2.

Different from the model (15), the subordinator in this model changes not only the convergence rate but also the stationary value of MSD for long times. The external force in this model damps the oscillation of the particle in harmonic potential and drags it towards zero for all times; while the subordinated process (15) does not get dragged to zero position during waiting times since the external force is zero during these time periods.

Let us pay attention to the critical time distinguishing two different scales in two models (15) and (29). It is $t=(2\gamma)^{-\frac{1}{\alpha}}$ in the first model, depending on the parameter $\alpha$ and influenced by the inverse subodinator $s(t)$. On the contrary, the critical time is $t=(2\gamma)^{-1}$ in the second model, which is independent of $\alpha$ and as same as that of original process $x(s)$ in (15). It means that the critical time in the second model is independent of the subordinator $s(t)$ and fully determined by the harmonic potential itself. On the other hand, the size relation between these two critical time is uncertain, depending on $\gamma$. If $\gamma<\frac{1}{2}$, the time during which harmonic potential does not work is longer in (15) than in model (29).

Taking Laplace transform in (29) yields $\langle y(\lambda_{1})y(\lambda_{2})\rangle=\frac{2\sigma}{(\lambda_{1}+\gamma)(\lambda_{2}+\gamma)(\lambda_{1}+\lambda_{2})^{\alpha}}$, from which, one arrives at the position autocorrelation function

[TABLE]

for $t_{2}\geq t_{1}$. This position autocorrelation function shows the process described by model (29) is non-stationary, as well as the process in model (16). Then using the MSD (34) and autocorrelation function (37) of process $y(t)$, one obtains the correlation coefficient of the stochastic process $y(t)$ for fixed $t_{1}$ and $t_{2}\rightarrow\infty$

[TABLE]

with $G_{1}(t_{1})=2\gamma t_{1}^{\alpha}/\alpha\cdot{}_{1}F_{1}(\alpha,1+\alpha;2\gamma t_{1})$ being a constant for fixed $t_{1}$. Different from (22), the correlation coefficient here exponentially decays to zero. It means that the positions at two different times are no longer long-range dependent as a result of the continuous influence of harmonic potential in model (29).

The ensemble averaged integrand in (25) could be obtained as

[TABLE]

This result not only depends on the lag time $\Delta$ but also the time $t$, implying the aging phenomenon by regarding $t$ as the aging time $t_{a}$ in a system. It means that the observation time impacts the statistical quantities of a system, which was initially prepared. After some calculations, the ensemble-time averaged MSD is obtained,

[TABLE]

Here, $M(\Delta)=\frac{\Delta^{\alpha+1}}{\alpha+1}~{}_{2}F_{2}(\alpha+1,1;\alpha+2,\alpha+1;-2\gamma\Delta)$ and ${}_{2}F_{2}(a,b;c,d;z)$ is the hypergeometric function Abramowitz and Stegun (1972).

For short $\Delta$, i.e., $\Delta\ll\gamma^{-1}$, the time averaged MSD is the same as the one for a free particle in (12), growing linearly with the lag time. For large $\Delta$, i.e., $\gamma^{-1}\ll\Delta\ll T$, by using the asymptotic expression of H-function Mathai et al. (2009), we have

[TABLE]

and find that the time averaged MSD approaches to a constant

[TABLE]

which is different from the time averaged MSD (26) in model (16). The disagreement between the ensemble and time averaged MSD, scaling as $\Delta^{\alpha-1}$ and $\Delta^{0}$ respectively, means the non-ergodicity of this system. Fig. 3 shows the consistency of the simulation results and analytical ones of the ensemble-time averaged MSD for different $\alpha$. It also can be found that the turning point is almost $\gamma^{-1}$ in Fig. 3, beyond which the plateau value (42) is approached. The ergodicity breaking parameter of this model is

[TABLE]

which is also similar to the one of free particle in (27) but with a larger coefficient.

IV Subdiffusive dynamics in linear potential

This section focuses on the influence of the linear potential acting on the original process in (2) and on the subordinated process in (3).

IV.1 Force acting on original process $x(s)$

We first consider the Langevin system with linear potential, i.e., a constant external force Cairoli and Baule (2015a); Dieterich et al. (2015); Chen et al. (2018)

[TABLE]

The corresponding single Langevin equation of $y(t)=x(s(t))$ in physical time is

[TABLE]

which evidently shows that the external force only acts at the moments of jump and does not affect the particle during waiting times (trapping events). In addition, the Fokker-Planck equation with respect to this Langevin system is also (24) by replacing $F(y)$ with $F$.

With the procedures similar to the case of harmonic potential in Section III, the first moment and MSD of the subordinated process $y(t)=x(s(t))$ are obtained as Compte et al. (1997); Metzler and Klafter (2000); Dieterich et al. (2015); Chen et al. (2018)

[TABLE]

with $0<\alpha<1$. The simulation results are presented in Fig. 4. This subordinated Langevin system shows subdiffusion when $0<\alpha<\frac{1}{2}$ and superdiffusion when $\frac{1}{2}<\alpha<1$. It is more or less interesting that the waiting time with infinite mean value produces superdiffusion. Compte et al. Compte et al. (1997) explained this phenomenon that some stagnated particles are not continuously dragged by the stream and thus slow down the advancement of the center of mass of the particles, hence the main dispersion mechanism should be convection. The deviation of the MSD of Langevin systems with constant force from the one of free particle implies this external force is a biasing force. In addition, the generalised Einstein relation Bouchaud and Georges (1990); Metzler et al. (1999); Barkai et al. (2000); Metzler and Klafter (2000); Blickle et al. (2007) connects the first moment of the particle displacements under a constant force to the second moment of a free particle, $\langle y(t)\rangle_{F}=\frac{F}{2k_{B}\mathcal{T}}\langle y^{2}(t)\rangle_{0}$. Here $k_{B}$ is the Boltzman constant and $\mathcal{T}$ is absolute temperature. Here we emphasize that the generalised Einstein relation holds for a subordinated process if it is valid for the original process. More precisely, the subordinator here affects the moments, simultaneously, as

[TABLE]

Let us now see the ergodicity of model (44). By using the technique of Laplace transform, we get the position autocorrelation function of $y(t)$ for $t_{2}\geq t_{1}$ as

[TABLE]

with a hypergeometric function ${}_{2}F_{1}(a,b;c;z)$ Abramowitz and Stegun (1972), which means that the process is non-stationary. Using the MSD and autocorrelation function of process $y(t)$, one can obtain the ensemble averaged integrand in (25)

[TABLE]

which shows the aging phenomenon because of the explicit dependence on $t$. Then we obtain the ensemble-time averaged MSD

[TABLE]

for $\Delta\ll T$. See the simulation results in Fig. 5. For short lag time, the time averaged MSD is linearly dependent on lag time $\Delta$ as the free subdiffusion case in (12), which can be explained by the fact that the particles are not affected by external forces in short times. With the increase of lag time, the time averaged MSD becomes super-linear in lag time, and is proportional to the square of the external force. The disagreement between time and ensemble averaged MSD shows the non-ergodicity of the stochastic process in (44), although it could exhibit “superdiffusion” phenomenon. Similar to the case of harmonic potential, the constant force in this model also only increases the coefficient of ergodicity breaking parameter,

[TABLE]

For model (44), there is the generalized Einstein relation, being similar to the CTRW model for the subdiffusive process He et al. (2008); Froemberg and Barkai (2013),

[TABLE]

where $\langle\overline{\delta^{1}(\Delta)}\rangle_{F}=\frac{1}{T-\Delta}\int_{0}^{T-\Delta}\langle y(t+\Delta)-y(t)\rangle_{F}dt$. From the discussions above, we recognize that if the generalized Einstein relation is satisfied by the original process $x(s)$ Bouchaud and Georges (1990); Barkai and Fleurov (1998), it will still be valid for a subordinated process $y(t)$ in (44).

As for the correlation of process $y(t)$ in (44), by using the asymptotic expression ${}_{2}F_{1}(a,b;c;z)\simeq 1+\frac{ab}{c}z$ for small $z$, one could obtain that

[TABLE]

for fixed $t_{1}$ and $t_{2}\rightarrow\infty$. Here $G_{2}(t_{1})$ could be regarded as a constant for fixed $t_{1}$. It indicates the long-range dependence of process $y(t)$, although the correlation is weaker than the one of free particle in (23).

IV.2 Force acting on subordinated process $y(t)$

What about the differences if the external constant force $F$ acts directly on the subordinated process $y(t)$ and continues to affect the stochastic process all the time. In this case, the Langevin equation is Cairoli and Baule (2015b); Eule and Friedrich (2009)

[TABLE]

with the equivalent coupled Langevin equation

[TABLE]

Using the solution of the exact trajectory, $y(t)=Ft+\sqrt{2\sigma}\int_{0}^{t}\overline{\xi}(\tau)d\tau$, the MSD of stochastic process $y(t)$ is

[TABLE]

which coincides with the MSD of free particle in (11); see Fig. 6 for the simulation results. It implies that the external force does not change the subdiffusion behavior and behaves as a decoupled force. Hence, the subdiffusion model (54) is Galilei invariant Metzler and Klafter (2000). In addition, the generalised Einstein relation in this case is not fulfilled since $\langle y(t)\rangle_{F}$ grows linearly with time $t$ while $\langle y^{2}(t)\rangle_{0}$ scales as $t^{\alpha}$.

The corresponding Fokker-Planck equation of model (54) is Cairoli et al. (2018)

[TABLE]

where

[TABLE]

is fractional substantial derivative Friedrich et al. (2006) with the Fourier-Laplace transform

[TABLE]

When $F=0$, $\mathcal{D}_{t}^{1-\alpha}$ recovers the Riemann-Liouville fractional derivative $D_{t}^{1-\alpha}$ and the Fokker-Planck equation goes back to the free subdiffusion case (13).

After some simple calculations, we obtain the position autocorrelation function of the Langevin system (54), $\langle(y(t_{1})-\langle y(t_{1})\rangle)(y(t_{2})-\langle y(t_{2})\rangle)\rangle=\frac{2\sigma}{\Gamma(1+\alpha)}t_{1}^{\alpha}$ for $t_{2}\geq t_{1}$, which reveals the non-stationary of this process. Then the correlation coefficient is obtained as

[TABLE]

for fixed $t_{1}$ and large $t_{2}$, which is identical with the one of free particle in (23), implying that the decoupled force does not affect the correlation of positions at two different times. In addition, the ensemble-time averaged MSD for $\Delta\ll T$ is

[TABLE]

which is the same as the free particle case in (12); see the simulation results in Fig. 7. Obviously, the ergodicity breaking parameter also remains unchanged. But the generalized Einstein relation is not satisfied when one compares $\langle\overline{\delta^{1}(\Delta)}\rangle_{F}=F\Delta$ with (58).

The above results show that the decoupled constant force, which affects a stochastic process all the time, will change the first moment of the stochastic process $y(t)$, but maintain the ensemble and time averaged MSD, as well as the correlation coefficient. The ergodic or non-ergodic behavior will not be changed in a decoupled force field. However, the Einstein relation is not valid any more due to the changes of the first moment.

Interestingly, the results for the case of the biased constant force in model (44) are completely opposite. The ensemble and time averaged MSD are changed while the exponent of ergodicity breaking parameter is maintained. In addition, the Einstein relation is still valid. On the other hand, the Fokker-Planck equation corresponding to the former model (44) includes the Riemann-Liouville fractional derivative and the Fokker-Planck operator while the latter model (54) involves the fractional substantial derivative.

V Subdiffusive dynamics in time-dependent potential

This section focuses on the influence of the time-dependent periodic oscillation force acting on the original process in (2) with white Gaussian noise and on the subordinated process in (3) with compound noise. Similar to another two kinds of forces discussed in Sec. III and Sec. IV, the essential difference is still that the time-dependent force here only acts at the moments of jump for the former case, but keeps acting on the system all the time for the latter case.

V.1 Force acting on original process $x(s)$

Consider the Langevin system with the time-dependent force acting on the original process $x(s)$, and it is expressed as Magdziarz et al. (2008); Weron and Magdziarz (2008)

[TABLE]

Noting the subordination of the original process $x(s)$, the force term $F(t(s(t)))=F(t)$ has reasonable physical meaning, since a physical force should act on a system at physical time $t$ not internal time $s$ Magdziarz et al. (2008). Its corresponding single Langevin equation describing the subordinated process $y(t)=x(s(t))$ in physical time is

[TABLE]

It is obvious that the time-dependent force $F(t)$ acts on the system only at the moments of jump, and the corresponding Fokker-Planck equation is Sokolov and Klafter (2006); Magdziarz et al. (2008); Eule and Friedrich (2009); Cairoli and Baule (2017); Sokolov and Klafter (2006); Magdziarz (2009)

[TABLE]

where the Riemann-Liouville fractional derivative $D_{t}^{1-\alpha}$ cannot be interchanged with $-\frac{\partial}{\partial y}F(t)+\sigma\frac{\partial^{2}}{\partial y^{2}}$.

Using the Fokker-Planck equation (61) derived from CTRW model, Sokolov et al. Sokolov and Klafter (2006) obtained the recursive relation of the moments $r_{n}(t):=\langle y^{n}(t)\rangle$

[TABLE]

with $r_{0}(t)=1,r_{-1}(t)=0,$ and $n\in\mathbb{N}$. Two years later, Magdziarz et al. Magdziarz et al. (2008) derived the same recursive relation of the moments through the Langevin equation (60). Hence, the correspondence between the Fokker-Planck equation and the Langevin equation with a time-dependent force is established.

Here we take an oscillating external force $F(t)=f_{0}\sin(\omega t)$. Although the ensemble averaged MSD of the stochastic process $y(t)$ can be obtained from the recursive relation (62), here we also present the precise results of the position autocorrelation function and time averaged MSD by using the Laplace transform method. Integrating (60) with respect to time $t$ leads to

[TABLE]

The position autocorrelation function of stochastic process $y(t)$ for $t_{2}\geq t_{1}$ is

[TABLE]

the detailed derivation of which is presented in Appendix A; and it shows the non-stationary property of the stochastic process $y(t)$. After long times, the asymptotic expression of position autocorrelation function can be obtained by taking $\lambda_{1},\lambda_{2}\rightarrow 0$ in (80) and making inverse Laplace transform,

[TABLE]

for $t_{2}\geq t_{1}$. Taking $t_{1}=t_{2}$ in (64), one obtains the second moment of $y(t)$,

[TABLE]

In order to see the fluctuation of the second moment more clearly, we simulate the first term denoted as $D_{1}(t)$ in the second moment (66). The accordance between the simulation result and the analytical result could be found in Fig. 8. For long times,

[TABLE]

which exhibits subdiffusion behavior, consistent with the asymptotic expression in Sokolov and Klafter (2006). The field-dependent contribution, which comes from the first term in (66), makes oscillation and additional dispersion of the particle position compared with the free subdiffusion case (11). And for long times, this additional dispersion grows sublinearly with time $t$ in (67). In addition, the first moment of stochastic process $y(t)$ for long times can be easily obtained as

[TABLE]

The oscillation of the mean value of process $y(t)$ tends to a constant for long times, which means the response to the external perturbation dies out for long times. It is also one of the manifestations of aging Sokolov et al. (2001).

Another equivalent expression of stochastic process $y(t)$ in (60) is

[TABLE]

Based on this expression, the position autocorrelation function of $y(t)$ can also be obtained. Detailed derivation is shown in Appendix A.

Combining the first two moments and the position autocorrelation function of the stochastic process $y(t)$, one obtains the correlation coefficient for fixed $t_{1}$ and $t_{2}\rightarrow\infty$

[TABLE]

which is the same as the free particle in (23). The essential reason for this interesting finding is that the oscillating external force changes the coefficient of position autocorrelation function and ensemble averaged MSD of the free subdiffusive process at the same degree.

The ensemble averaged integrand in (25) for long times is

[TABLE]

the dependence of which on time $t$ implies the aging phenomenon of this Langevin system. The ensemble-time averaged MSD for $\Delta\ll T$ is

[TABLE]

See the simulation results in Fig. 9. Comparing with the case of free particle in (12), the oscillating external force here adds an additional contribution on the time averaged MSD, which grows linearly with the lag time $\Delta$. The disagreement between the time and ensemble averaged MSD, which scale as $\Delta$ and $\Delta^{\alpha}$ respectively, indicates the non-ergodicity behavior of the stochastic process. The ergodicity breaking parameter here is as same as the case of free particle due to the same degree of the changes on the coefficients of ensemble and time averaged MSD.

V.2 Force acting on subordinated process $y(t)$

Next, we consider the case in which the time-dependent external force acting on the system all the time, i.e.,

[TABLE]

The corresponding coupled Langevin equation is Eule and Friedrich (2009)

[TABLE]

In order to compare with the previous model (60), we also consider the oscillating force $F(t)=f_{0}\sin(\omega t)$ here. The firstly moment of stochastic process $y(t)$ is

[TABLE]

which shows a significant difference with the constant mean value (68) of the model (60). The mean value here keeps oscillation at a fixed frequency $\omega$ with the evolution of time since the oscillating external force $F(t)=f_{0}\sin(\omega t)$ influences this system for the whole time. The MSD of this model can be easily obtained as

[TABLE]

which is identical with the case of free particle in (11). It means that the time-dependent external force field here acts as a decoupled force, independent of the diffusion behavior. It is Galilean invariant model while another model (59) breaks Galilean invariance. The Fokker-Planck equation corresponding to the Langevin system (72) is

[TABLE]

with the operator in Fourier space

[TABLE]

See the detailed derivations in Appendix B. Taking the constant force $F(t)=F$, the operator $\mathcal{A}_{t}^{1-\alpha}$ reduces to the fractional substantial derivative $\mathcal{D}_{t}^{1-\alpha}$ and the Fokker-Planck equation goes back to (57).

As for the time averaged MSD, using the first two moments (74) and (75), together with the position autocorrelation function

[TABLE]

for $t_{2}>t_{1}$, one could obtain the time averaged MSD for $\Delta\ll T$

[TABLE]

which indicates the non-ergodic behavior of this Langevin system. It is also consistent with the case of free particle in (12). The simulation results could be found in Fig. 10. In addition, the ergodicity breaking parameter and the correlation coefficient of this model are also unchanged, compared with the case of free particle.

In conclusion, the time-dependent force in (72) is a decoupled force, just as the decoupled constant force in (54). However, the time-dependent force in (59) plays as a biasing force, just as the biased constant force in (44). They all change the ensemble and time averaged MSD of the free particle. One difference is that the time-dependent force keeps the same correlation coefficient as the free particle case, while the constant force weakens the correlation. In addition, the Fokker-Planck equation corresponding to the former model (59) includes the Riemann-Liouville fractional derivative while the latter model (72) involves a novel fractional derivative similar to the fractional substantial derivative.

VI Summary

This paper focuses on the subdiffusion in an external force field. We mainly consider the influences of different patterns of external forces acting on the original process $x(s)$ or on the subordinated process $y(t)$. For this, we choose three kinds of common external forces — linear force, constant force, and time-dependent oscillation force, and investigate some important statistical quantities depending on one-point or two-point PDF, such as ensemble and time averaged MSD, position autocorrelation function, correlation coefficient, and so on. There are obvious discrepancies between position-dependent and position-independent forces.

One example of the position-dependent forces is the linear force (harmonic potential) in this paper. If it acts on original process $x(s)$, the ensemble averaged MSD tends to a non-zero constant for long times, while it tends to zero at power-law rate if this force acts on the subordinated process $y(t)$. The essential reason is that the external force drags the particle to zero position for all times in the latter case. These two stochastic processes are both non-ergodicity, non-stationary, and aging. However, the position at two different times in the former Langevin system is long-range dependent while in the latter it is not due to the continuous effects of external force. The position-dependent force acting on the original process $x(s)$ or on the process $y(t)$ does not affect the power of ergodicity breaking parameter, compared with the case of free particle.

As for the position-independent force, such as, constant force and periodic oscillation force in this paper, we find that it behaves as a biasing force if acting on the original process $x(s)$ and does change the ensemble and time averaged MSD, as well as the position autocorrelation function. One unexpected finding is the correlation coefficient — it becomes weaker at the effects of constant force, while it remains unchanged in the case of the periodic oscillation force. On the other hand, if the position-independent force acts on the subordinated process $y(t)$, it behaves as a decoupled force and does not make change of the statistical quantities we study.

Especially, the Fokker-Planck equations of the processes in different force field are different, being embodied by fractional derivatives; it is Riemann-Liouville type if the external force affects the process only at the moments of jump, but fractional substantial derivative or a similar novel fractional derivative if the external force keeps acting on the process all the time.

Collecting the properties of the statistical quantities and the discrepancies between different forces or different patterns the forces act on the system, we have a global knowledge of the motion of the subdiffusive particles in an external force field. This will help us to distinguish a large amount of processes with some similar features.

Acknowledgments

This work was supported by the National Natural Science Foundation of China under grant no. 11671182, and the Fundamental Research Funds for the Central Universities under grant no. lzujbky-2018-ot03.

Appendix A Derivation of the autocorrelation function (64)

The position autocorrelation function of process $y(t)$ is

[TABLE]

where the cross-terms are zero due to the independence of $B(s)$ and $s(t)$ and could be omitted. For the first term in (77), denoted as $\langle y(t_{1})y(t_{2})\rangle_{1}$, it can be dealt with in Laplace space. Using the expression of

[TABLE]

and the correlation function of inverse subordinator $s(t)$ in Laplace space Baule and Friedrich (2005)

[TABLE]

one get

[TABLE]

where $\lambda_{j}^{\pm}=\lambda_{j}\pm iw$ with $j=1,2$. After performing the inverse Laplace transform, the first term of position autocorrelation function is

[TABLE]

where

[TABLE]

for $t_{2}\geq t_{1}$. The second term of $\langle y(t_{1})y(t_{2})\rangle$ is

[TABLE]

Here we use the independence of $B(s)$ and $s(t)$, as well as $\langle s(t)\rangle=\frac{1}{\Gamma(1+\alpha)}t^{\alpha}$ Baule and Friedrich (2005). Finally, the position autocorrelation function of stochastic process $y(t)$ for $t_{2}\geq t_{1}$ is (64).

We now present another exact expression of $y(t)=x(s(t))$. Similar to the method in Cairoli and Baule (2015b),

[TABLE]

where the last equality is obtained by using $\delta(s-s(t))=-\frac{\partial}{\partial s}\Theta(t-t(s))$ Baule and Friedrich (2005) and integration by parts. Then the differential expression of process $y(t)$ is

[TABLE]

The first term of the autocorrelation function of $y(t)$ (80) in Laplace space could also be obtained from the first term in (83):

[TABLE]

where $p(t(s_{1}),t(s_{2}),s_{1},s_{2})$ is the two-point joint PDF of subordinator $t(s)$. Using the expression of $p(t_{1},t_{2},s_{1},s_{2})$ in Laplace space Baule and Friedrich (2005), i.e.,

[TABLE]

and (78), we can also obtain the term (80) of autocorrelation function of $y(t)$ in Laplace space.

Appendix B Derivation of the Fokker-Planck equation (76)

By the similar method shown in Cairoli and Baule (2017, 2015a), we now derive the Fokker-Planck equation corresponding to the Langevin equation $\dot{y}(t)=F(t)+\sqrt{2\sigma}\overline{\xi}(t)$, equivalently, $y(t)=\int_{0}^{t}F(t^{\prime})dt^{\prime}+\sqrt{2\sigma}B(s(t))$. As we all know, $\int_{0}^{t}F(t^{\prime})dt^{\prime}$ is a process with finite variation and $B(s(t))$ is a martingale, which lead process $y(t)$ to be a semi-martingale Kunita (1997). In addition, $y(t)$ has continuous path. Hence we can use its Itô formula as follows Kunita (1997)

[TABLE]

where $[y,y]_{t}=\sum_{i}|y(t_{i})-y(t_{i-1})|^{2}=2\sigma\int_{0}^{t}ds(\tau)$ is the quadratic variation of process $y(t)$ Cairoli and Baule (2017); Øksendal (2005), and it could be gotten by $dtdt=dtdB_{t}=dB_{t}dt=0$. The PDF of process $y(t)$ in Fourier space is $p(k,t)=\mathcal{F}_{y\rightarrow k}[\langle\delta(y-y(t))\rangle]=\langle e^{iky(t)}\rangle$. So we take $f(y(t))=e^{iky(t)}$. Then

[TABLE]

Taking the ensemble average of (85), making inverse Fourier transform, and taking partial derivative with respect to $t$, we obtain

[TABLE]

The last term $\langle\delta(y-y(t))\dot{s}(t)\rangle$ could be dealt with in Fourier space

[TABLE]

Here we use the fact $\int_{0}^{t}\delta(s-s(\tau))ds(\tau)=\Theta(t-t(s))$ since $\Theta(t-t(s))=1-\Theta(s-s(t))$ Baule and Friedrich (2005, 2007). For simplicity of notation, let us define $G_{1}(k,t)=\int_{0}^{\infty}\langle e^{ik\sqrt{2\sigma}B(s)}\rangle\langle\delta(t-t(s))\rangle ds$. Taking Laplace transform, one has

[TABLE]

Similarly, the PDF of $y(t)$ in Fourier space is

[TABLE]

We define $G_{2}(k,t)=\int_{0}^{\infty}\langle e^{ik\sqrt{2\sigma}B(s)}\rangle\langle\delta(s-s(t))\rangle ds$, the Laplace transform of which is

[TABLE]

Here we use the PDF of inverse subordinator $s(t)$ in Laplace space (5). Therefore, $G_{1}(k,\lambda)=\lambda^{1-\alpha}G_{2}(k,\lambda)$, i.e.,

[TABLE]

Then one has

[TABLE]

Finally, the Fokker-Planck equation corresponding to the Langevin system $\dot{y}(t)=F(t)+\sqrt{2\sigma}\overline{\xi}(t)$ is

[TABLE]

where the operator $\mathcal{A}_{t}$ in Fourier space is $\mathcal{F}_{y\rightarrow k}[\mathcal{A}_{t}p(y,t)]=e^{ik\int_{0}^{t}F(t^{\prime})dt^{\prime}}D_{t}^{1-\alpha}e^{-ik\int_{0}^{t}F(t^{\prime})dt^{\prime}}p(k,t)$.