Probability density of the fractional Langevin equation with reflecting walls

TL;DR

This paper studies how anomalous diffusion described by the fractional Langevin equation behaves in confined spaces with reflecting walls, revealing unique distribution patterns influenced by force correlations.

Contribution

It provides a detailed analysis of the probability density evolution for fractional Langevin processes with reflecting boundaries, highlighting differences from normal diffusion and fractional Brownian motion.

Findings

Probability density converges to uniform distribution on finite intervals.

Persistent correlations cause particle accumulation at the wall.

Anti-persistent correlations lead to particle depletion near the wall.

Abstract

We investigate anomalous diffusion processes governed by the fractional Langevin equation and confined to a finite or semi-infinite interval by reflecting potential barriers. As the random and damping forces in the fractional Langevin equation fulfill the appropriate fluctuation-dissipation relation, the probability density on a finite interval converges for long times towards the expected uniform distribution prescribed by thermal equilibrium. In contrast, on a semi-infinite interval with a reflecting wall at the origin, the probability density shows pronounced deviations from the Gaussian behavior observed for normal diffusion. If the correlations of the random force are persistent (positive), particles accumulate at the reflecting wall while anti-persistent (negative) correlations lead to a depletion of particles near the wall. We compare and contrast these results with the strong accumulation and depletion effects recently observed for non-thermal fractional Brownian motion with reflecting walls, and we discuss broader implications.

Figures

Equations

$$m\frac{d^2}{d{t}^2}x(t)=-\bar{\gamma }\frac{d}{dt}x(t)+{\xi }_w(t).$$

$$m\frac{d^2}{d{t}^2}x(t)=-\bar{\gamma }{\int }_0^{t}d{t}^{\prime }\mathcal{K}(t-t^{\prime })\frac{d}{d{t}^{\prime }}x(t^{\prime })+\xi (t).$$

$$⟨\xi (t)\xi (t^{\prime })⟩=k_{B}T\bar{\gamma }\mathcal{K}(t-t^{\prime }).$$

$$⟨\xi (t)\xi (t^{\prime })⟩\sim \alpha (\alpha -1)K_{\alpha }|t-t^{\prime }{|}^{\alpha -2}$$

$$⟨x^2(t)⟩\sim t^{2-\alpha },$$

$$V(x)=V_0\exp [\mp \lambda (x-x_0)].$$

$$F(x)=-dV(x)/dx=\pm F_0\exp [\mp \lambda (x-x_0)]$$

$$v_{n+1}$$

$$x_{n+1}$$

$$⟨{\xi }_m{\xi }_{m+n}⟩=K_{\alpha }(\Delta t{)}^{\alpha -2}\left(|n-1{|}^{\alpha }-2|n{|}^{\alpha }+|n+1{|}^{\alpha }\right).$$

$$⟨{\xi }_m{\xi }_{m+n}⟩=T{\mathcal{K}}_n.$$

$$P(x,t)=\frac{1}{\sigma (t)}Y\left[\frac{x}{\sigma (t)}\right]=\frac{1}{b{t}^{1-\alpha /2}}Y\left[\frac{x}{b{t}^{1-\alpha /2}}\right]$$

$$x_n\to 2x_0-x_n,v_n\to -v_n$$

$$S_n=\sum _{m=0}^n{\mathcal{K}}_{n-m}{v}_m,$$

$$I_n(i,j)=\sum _{m=i-j}^{i+j}{\mathcal{K}}_{n-m}{v}_m.$$

$${\mathcal{K}}_{n-m}={\mathcal{K}}_{n-i}-(m-i){\mathcal{K}}_{n-i}^{{}^{\prime }}\pm \dots$$

$$I_n(i,j)={\mathcal{K}}_{n-i}\sum _{m=i-j}^{i+j}{v}_m-{\mathcal{K}}_{n-i}^{{}^{\prime }}\sum _{m=i-j}^{i+j}(m-i)v_m\pm \dots .$$

Full text

I Introduction

Normal diffusion is an almost omnipresent transport process that can be found in biological, chemical, as well as physical systems, and beyond. It is characterized by a linear relation between the mean-square displacement $\langle x^{2}\rangle$ of the diffusing object and the elapsed time $t$. According to Einstein Einstein (1956), diffusive behavior results from stochastic motion that is local in both time and space. This means individual random displacements (step lengths) have a finite correlation time and a finite variance.

In recent years, considerable attention has been attracted by anomalous diffusion, i.e., random motion that does not fulfill the dependence $\langle x^{2}\rangle\sim t$. Both subdiffusion (for which $\langle x^{2}\rangle$ grows slower than $t$) and superdiffusion (where $\langle x^{2}\rangle$ grows faster than $t$) have been observed in numerous experiments and described mathematically (for reviews see, e.g., Refs. Metzler and Klafter (2000); Höfling and Franosch (2013); Bressloff and Newby (2013); Metzler et al. (2014); Meroz and Sokolov (2015); Metzler et al. (2016); Nørregaard et al. (2017) and references therein).

Anomalous diffusion can arise via several different mechanisms if the condition of locality in space and time is violated. For example, if the step lengths are broadly distributed such that their variance does not exist, the motion can become superdiffusive. Conversely, a broad distribution of waiting times with a diverging mean may lead to subdiffusive motion. Another important mechanism that leads to anomalous diffusion is the presence of long-range power-law correlations in time between individual random displacements (steps). Such correlations can produce subdiffusive or superdiffusive behavior even if the step lengths and waiting times are narrowly distributed. Fractional Brownian motion (FBM) and the fractional Langevin equation (FLE) are two prototypical mathematical models for this situation.

FBM was introduced by Kolmogorov Kolmogorov (1940) and further investigated by Mandelbrot and van Ness Mandelbrot and Ness (1968). It has been studied in the mathematical literature quite extensively (see, e.g., Refs. Kahane (1985); Yaglom (1987); Beran (1994); Biagini et al. (2008)), and it has been applied to a wide variety of problems including polymer dynamics Chakravarti and Sebastian (1997); Panja (2010), diffusion inside living cells Szymanski and Weiss (2009); Jeon et al. (2011a), traffic in electronic networks Mikosch et al. (2002), as well as stock markets Comte and Renault (1998); Rostek and Schöbel (2013). FBM is a self-similar, non-Markovian Gaussian process with stationary increments $\xi$ that feature long-range power-law correlations in time. If the increments are positively correlated (persistent), the resulting motion is superdiffusive whereas anticorrelated (antipersistent) increments produce subdiffusive motion.

Recently, large-scale computer simulations of FBM confined by reflecting walls have demonstrated that the interplay between the long-time correlations and the confinement strongly affects the probability density function $P(x,t)$ of the diffusing particle. Specifically, if the motion is restricted to a semi-infinite interval by a reflecting wall at the origin, the probability density becomes highly non-Gaussian and develops a power-law singularity, $P\sim x^{\kappa}$ at the wall Wada and Vojta (2018); Wada et al. (2019). For persistent noise, particles accumulate at the wall, $\kappa<0$. In contrast, particles are depleted close to the wall, $\kappa>0$ for anti-persistent noise. Analogous simulations of FBM on a finite interval, with reflecting walls at both ends, have established that the stationary probability density reached for long times strongly deviates from the uniform distribution found for normal diffusion Guggenberger et al. (2019). Particles accumulate at the walls and are depleted in the middle of the interval for persistent noise whereas the opposite is true for anti-persistent noise.

FBM can be understood as random motion governed by external noise Klimontovich (1995). It does not obey the fluctuation-dissipation theorem and generally does not reach a thermal equilibrium state. It is thus interesting to ask whether the unusual features of the probability density close to reflecting walls survive or disappear if the long-range correlated noise fulfills the fluctuation-dissipation theorem.

To answer this question, we study in this paper the fractional Langevin equation (FLE) in the presence of reflecting potentials that restrict the motion to finite or semi-infinite intervals. Our results can be summarized as follows. As the random and damping forces fulfill the appropriate fluctuation-dissipation relation, the FLE motion confined to a finite interval is expected to approach the thermal equilibrium state for long times. We show here that this is indeed the case. Consequently, the probability density of the particle position approaches a uniform distribution. In contrast, the probability density on a semi-infinite interval shows pronounced deviations from the Gaussian behavior observed for normal diffusion. If the correlations of the random force are persistent, particles accumulate at the reflecting wall while anti-persistent correlations lead to a depletion of particles near the wall, compared to the naively expected Gaussian distribution.

In addition, we analyze a generalized Langevin equation featuring long-range correlated fractional noise but a conventional, instantaneous damping term. For this equation, which violates the fluctuation-dissipation theorem, we find that the probability density on a finite interval is non-uniform and resembles the result for FBM.

Finally, we compare reflecting potentials (as used in the bulk of our work) with hard reflecting boundary conditions for which the velocity simply changes sign when the particle hits the wall, $v\to-v$. Such reflecting boundary conditions were recently employed by Holmes Holmes (2019) who reported a non-uniform stationary probability density for the FLE on a finite interval. Our simulations demonstrate that the FLE with hard reflecting boundary conditions is very sensitive to the time step in the numerical integration. A uniform distribution, as required by thermal equilibrium, is only recovered for extremely small time steps.

The remainder of our paper is organized as follows. We introduce the FLE in Sec. II where we also describe our simulation method. Section III is devoted to the results for the FLE on a finite interval, whereas Sec. IV presents the behavior for a semi-infinite interval. Anti-persistent noise and damping forces that violate the fluctuation dissipation theorem are discussed in Sec. V. We conclude in Sec. VI. The appendix is devoted to a discussion of the boundary conditions as well as an algorithm to speed up the solution of the FLE.

II Fractional Langevin equation

II.1 Definition

Our starting point is the well-known Langevin equation Langevin (1908),

[TABLE]

It describes the time evolution of the position $x$ of a particle of mass $m$ moving under the influence of an uncorrelated random force (Gaussian white noise) $\xi_{w}(t)$ and a linear damping force with damping coefficient $\bar{\gamma}$.

If the random force is correlated (nonwhite noise) and the damping force is nonlocal in time, the motion follows the generalized Langevin equation Zwanzig (2001); Hänggi (1978); Goychuk (2012)

[TABLE]

In thermal equilibrium at temperature $T$, the fluctuation-dissipation theorem Kubo (1966) requires that the noise correlation function and the damping (memory) kernel $\mathcal{K}$ are related via

[TABLE]

We are interested in the case of $\xi(t)$ being a fractional Gaussian noise Qian (2003), i.e., a Gaussian process of zero mean, $\langle\xi(t)\rangle=0$ and a power-law correlation function

[TABLE]

for $t\neq t^{\prime}$. In this case, eq. (2) is called the FLE Lutz (2001). Here, $K_{\alpha}$ denotes the noise amplitude. The correlation exponent $\alpha$ is restricted to the range $1<\alpha<2$ because the damping integral in (2) diverges at $t=t^{\prime}$ for $\alpha<1$, and $\alpha>2$ is unphysical as it corresponds to correlations that increase with time. (The Hurst exponent $H$, often used in the mathematical literature, is given by $H=\alpha/2$.) In the range $1<\alpha<2$, the fractional Gaussian noise is persistent (positively correlated). The limiting case, $\alpha=1$, corresponds to normal diffusion.

The behavior of the FLE in the absence of confining potentials is well understood (see, e.g., Ref. Metzler et al. (2014) and references therein). Consider, e.g., a particle starting from rest at the origin at time $t=0$. The probability densities of both its velocity and its position are Gaussian. For long times, the mean-square velocity approaches the value $k_{B}T/m$ prescribed by the classical equipartition theorem whereas the mean-square displacement crosses over from ballistic behavior, $\langle x^{2}(t)\rangle\sim t^{2}$, at short times to anomalous diffusion,

[TABLE]

at long times. This means that the FLE with persistent noise, $1<\alpha<2$, leads to subdiffusion while FBM with the same noise produces superdiffusion, $\langle x^{2}(t)\rangle\sim t^{\alpha}$. This is caused by the fact that the fluctuation-dissipation relation (3) couples persistent noise with long-range memory in the damping term.

To confine the motion to either a finite or a semi-infinite interval, we introduce reflecting walls via repulsive potentials. Specifically, to model a reflecting wall at position $x_{0}$, we introduce an external potential

[TABLE]

The $\pm$ sign in the exponent distinguishes walls reflecting from the right or from the left. This potential creates an additional external force

[TABLE]

with $F_{0}=\lambda V_{0}$ on the r.h.s. of the FLE (2). We are interested in the limit of large $\lambda$ for which the decay length $\lambda^{-1}$ of the wall force is small compared to the considered interval lengths $L$. In this limit, particles in the interior of the interval will not experience an external force. The slight “softening” of the reflecting wall due to a finite $\lambda$ is needed for the consistent numerical integration of the FLE as described in Sec. II.2. For more details and a comparison with alternative implementations of reflecting walls in the FLE, see Appendix A.

II.2 Simulation method

To simulate the FLE (2), we set the mass $m$, the damping coefficient $\bar{\gamma}$, and the Boltzmann constant $k_{B}$ to unity. We then discretize time $t_{n}=n\Delta t$ with $n=0,1,\ldots,N_{t}$. After replacing the time derivatives with finite-difference expressions, the FLE turns into the recursion relations

[TABLE]

The $\xi_{n}$ are a discrete fractional Gaussian noise Qian (2003), i.e., identical Gaussian random numbers of zero mean and correlation function

[TABLE]

These correlated random numbers are calculated using the Fourier-filtering method Makse et al. (1996). It starts from a sequence of independent Gaussian random numbers $\chi_{n}$. The Fourier transform $\tilde{\chi}_{\omega}$ of these numbers is then converted via ${\tilde{\xi}_{\omega}}=[\tilde{C}(\omega)]^{1/2}\tilde{\chi}_{\omega}$, where $\tilde{C}(\omega)$ is the Fourier transform of the correlation function (10). The inverse Fourier transformation of the ${\tilde{\xi}_{\omega}}$ gives the desired fractional Gaussian noise.

The noise and the damping kernel fulfill the discrete version of the fluctuation-dissipation relation

[TABLE]

A naive implementation of the recursion (8) is numerically expensive because the evaluation of a single damping integral scales linearly with the number of time steps. Thus the total effort grows quadratically with the number of time steps. We have devised an improved algorithm that speeds up the evaluation of the damping integrals by several several orders of magnitude. This algorithm is discussed in Appendix B.

Choosing a suitable value for the time step $\Delta t$ is crucial for the performance of the simulation. On the one hand, $\Delta t$ needs to be small enough to limit the error due to the time discretization 111The choice of $\Delta t$ is affected by the magnitude and steepness of the wall force (7) with larger values of $F_{0}$ and $\lambda$ requiring smaller $\Delta t$.. On the other hand, a small $\Delta t$ increases the numerical effort to reach long simulation times. To optimize the time step, we study the dependence on $\Delta t$ of the mean-square velocity $\langle v^{2}\rangle$ at long times of the FLE restricted to the semi-infinite interval $(0,\infty)$. The inset of Fig. 1 shows $\langle v^{2}\rangle$ vs. $\Delta t$ for $\alpha=1.5$ at temperature $T=1$.

The data show that the deviation of $\langle v^{2}\rangle$ from the value of 1 (required by the fluctuation-dissipation theorem) decreases with decreasing $\Delta t$, as expected. Based on these results and analogous calculations for other values of the correlation exponent $\alpha$, we have chosen a time step of $\Delta t=0.01$ for the majority of our simulations. As a further test, we determine how the probability density $P(x,t)$ is affected by the time step $\Delta t$. The main panel of Fig. 1 presents results for $\alpha=1.5$ and $t=20972$. They show that the probability densities for time steps $\Delta t=0.01$, 0.005, and 0.0025 agree very well, giving us further confidence in using $\Delta t=0.01$.

Using the above numerical method, we perform simulations for several values of the correlation exponent $\alpha=1.0,1.1,1.2,1.4,1.5,1.6$ and 1.7. We use up to $N_{t}=2^{27}\approx 134$ million time steps. All data are averages over a large number of trajectories ranging between $10^{4}$ and $5\times 10^{6}$. Unless otherwise noted, our simulations are performed for noise amplitude $K_{\alpha}=[\alpha(\alpha-1)]^{-1}$ and a temperature $T=1$. To model the reflecting walls, we use parameters $F_{0}=5$ and $\lambda=5$ in eq. (7).

III Finite interval

We now turn to the results, starting with simulations of the FLE confined to an interval $(-L,L)$ by reflecting potentials (6) with the appropriate signs located at $x_{0}=\pm L$. In these simulations, the particles start from rest at the origin, $x=0$, at time $t=0$. We then follow the time evolution until a steady state is reached.

Figure 2 presents the time dependence of the mean-square velocity $\langle v^{2}\rangle$ and the mean-square displacement $\langle x^{2}\rangle$ for several values of the correlation exponent $\alpha$ including the uncorrelated case, $\alpha=1$.

For all $\alpha$, the mean-square velocity, shown in the inset, quickly settles on the value 1, as expected from the fluctuation-dissipation theorem at unit temperature. The mean-square displacement, shown in the main panel of the figure, features more interesting properties. After transient ballistic behavior at very short times, it grows following the same power law, $\langle x^{2}\rangle\sim t^{2-\alpha}$, as the mean-square displacement of the free FLE, see eq. (5). At long times, it saturates at a constant value that is independent of the correlation exponent $\alpha$ and coincides with the expectation $\langle x^{2}\rangle=L^{2}/3$ for a uniform distribution of particles over the interval $(-L,L)$. These results agree with earlier findings by Jeon et al. Jeon and Metzler (2010).

The time evolution of the probability density $P(x,t)$ is shown in Fig. 3 for $\alpha=1.5$.

At short times, $P(x,t)$ is a Gaussian in $x$ centered at the starting point $x=0$ that broadens with time. Initially, this Gaussian is practically unaffected by the walls. For longer times, $P(x,t)$ crosses over to a stationary distribution. In thermal equilibrium, the stationary probability density is expected to follow the Boltzmann distribution, i.e., $P(x)$ should be proportional to $\exp(-V(x)/T)$. As $V(x)$ vanishes in the interior of the interval, we expect a uniform distribution. Very close to the wall (at distances of the order of $\lambda^{-1}$) we expect $P(x)$ to be “rounded” due to the wall potential (6). Specifically, the potential energy of a particle increases upon approaching the wall, suppressing $P(x)$. As $\lambda$ is a fixed system parameter that does not scale with the interval length $L$, the rounding becomes negligible in the scaling limit $\lambda L\gg 1$. The data in Fig. 3 are in perfect agreement with these expectations. For other values of the correlation exponent $\alpha$, the probability density $P(x,t)$ behaves in exactly the same way.

We emphasize that the behavior of the FLE on a finite interval observed here is very different from the behavior of FBM on a finite interval reported in Ref. Guggenberger et al. (2019). Whereas the probability density of the FLE reaches a uniform distribution for long times for all values of the correlation exponent $\alpha$, the stationary probability density of FBM is not uniform and depends on the value of $\alpha$. (The FBM probability density is uniform only for uncorrelated noise, $\alpha=1$). This difference is a direct consequence of the fact that the FLE fulfills the fluctuation-dissipation theorem and can thus reach thermal equilibrium whereas FBM does not fulfill the fluctuation dissipation theorem (and the fractional Gaussian noise is interpreted as external noise Klimontovich (1995)). To test this further, we also perform simulations of a generalized Langevin equation (2) with a damping force that does not fulfill the fluctuation-dissipation relation (3). This is discussed in Sec. V.2.

IV Semi-infinite interval

We now turn to the results for the FLE confined to the semi-infinite interval $(0,\infty)$ by a single reflecting wall at the origin. Because the motion is unbounded, we expect the mean-square displacement of a particle that starts at the origin to grow without limit with time. Figure 4 presents the time evolution of the mean-square displacement $\langle x^{2}\rangle$ and mean-square velocity $\langle v^{2}\rangle$ for several values of the correlation exponent $\alpha$.

The figure shows that the mean-square velocity quickly settles on the value $\langle v^{2}\rangle=1$ expected from the classical equipartition theorem. The mean-square displacement follows the same power law $\langle x^{2}\rangle\sim t^{2-\alpha}$ as is observed for the free (unconfined) FLE.

Because the FLE on a semi-infinite interval does not reach a steady state and thus not thermal equilibrium, the functional form of $P(x,t)$ is not constrained to follow the Boltzmann distribution. In fact, our simulations show that $P(x,t)$ has quite unexpected features. For reference, let us first discuss the case of normal diffusion (which corresponds to $\alpha=1$). The probability density of normal diffusion with a reflecting wall at the origin can be found by solving the diffusion equation, $\partial_{t}P=D\partial_{x}^{2}P$ under the flux-free boundary condition $\partial_{x}P=0$ at $x=0$. This yields a Gaussian distribution of the same width as in the unconfined case but with twice the amplitude (because $P(x,t)$ is restricted to nonnegative values).

Figure 5 presents simulation results for the probability densities $P(x,t)$ comparing a free, unconfined FLE to the FLE on our semi-infinite interval $(0,\infty)$ for correlation exponent $\alpha=1.5$. (In both cases, the particles start at the origin at $t=0$.)

For the free FLE, the simulations confirm the Gaussian functional form of $P(x,t)$. The probability density on the semi-infinite interval, however, displays significant deviations from the naively expected half-Gaussian form. We observe analogous non-Gaussian behavior for all studied $\alpha>1$.

The left panel of Fig.  6 presents the time evolution of the probability density for $\alpha=1.5$.

It shows that the non-Gaussian character persists to the longest times. In fact, if $x$ is scaled with the root-mean-square displacement $\sigma(t)=\sqrt{\langle x^{2}(t)\rangle}$ at each time $t$, the probability densities for all times collapse onto a common master curve, as is demonstrated in the right panel of Fig.  6. This implies that the probability density fulfills the scaling form

[TABLE]

where $Y$ is a dimensionless scaling function and $b$ is a constant. The small deviations from a perfect scaling collapse that can be seen for small $x$ at short times $t$ can be attributed to the “soft” wall potential that rounds $P(x,t)$ over a distance of order $\lambda^{-1}$ close to the wall. This rounding becomes negligible for long times, i.e., in the scaling regime $\lambda\sigma(t)\gg 1$.

The large-$x$ tail of the probability density takes a Gaussian form as can be clearly seen in Fig. 7 where we replot the data of Fig. 6 as $\log P(x,t)$ vs. $x^{2}$ such that a Gaussian leads to a straight line.

For small $x$, in contrast, the probability density is increased compared to the Gaussian, i.e., particles accumulate at the wall.

We observe analogous behavior for all values of the correlation exponent $\alpha$ studied in the range $1<\alpha<2$. The magnitude of the particle accumulation close to the wall depends on the $\alpha$ value . This is illustrated in Fig. 8 which shows the scaled probability density at long times for several $\alpha$ values.

For the case of normal diffusion, (uncorrelated random forces, $\alpha=1$) our simulations yield a Gaussian probability density as expected from the solution of the diffusion equation with a flux-free boundary condition at the origin. With increasing $\alpha$, i.e., as the correlations become more long-ranged, the deviations of $P(x,t)$ from Gaussian behavior increase. The exact functional form of $P(x,t)$ at the reflecting wall is hard to determine from the available numerical data as the rounding due to the soft wall force limits the resolution close to the wall. However, we see no indications of a divergence for $x/\sigma(t)\to 0$. Instead, the data appear compatible with a cusp at $x=0$.

It is interesting to compare the non-Gaussian behavior of the probability density of the FLE on a semi-infinite-interval with the non-Gaussian behavior of reflected FBM studied in Refs. Wada and Vojta (2018); Wada et al. (2019). In both cases, persistent, positively correlated noise ($\alpha>1$) leads to an accumulation of particles at the wall. Moreover, in both cases, the probability density fulfills the scaling form (12). However, the accumulation is qualitatively stronger for FBM for which the scaled probability density $\sigma(t)P(x,t)$ diverges at the wall. This divergence is well described by the power law $P(x,t)\sim x^{\kappa}$ with $\kappa=2/\alpha-2$; it is crucially important in applications that are dominated by rare events (see Ref. Wada et al. (2018) for an example). For the FLE, the accumulation of particles at the wall is weaker, and the scaled probability density remains finite at $x=0$. We emphasize that persistent noise leads to an accumulation at the wall for both FLE and FBM despite the fact that the motion is subdiffusive for the FLE and superdiffusive for FBM. Subdiffusive FBM which is caused by antipersistent noise, in contrast, leads to a depletion of particles near the wall. This suggests that the sign of the correlation determines whether or not particles are accumulated or depleted and not the type of anomalous diffusion (superdiffusion vs. subdiffusion). We note that FBM in a steep confining potential also shows an accumulation of probability away from the origin (i.e., close to the walls) in the persistent noise case, as compared to the corresponding Boltzmann distribution Guggenberger et al. (2019).

V Extensions

V.1 Anti-persistent noise

In the FLE defined by eqs. (2), (3), and (4), the correlation exponent $\alpha$ is restricted to $1<\alpha<2$, i.e., to the case of persistent noise. Anti-persistent noise (corresponding to $\alpha$ values below unity) is impossible because the damping integral diverges for $\alpha<1$. However, as the divergence arises from the short-time behavior at $t=t^{\prime}$, it can be cut off without changing the physically important long-time behavior of the noise correlations and the damping kernel. In fact, in the discretized version of the FLE defined in eqs. (8), (9), (10), and (11), the singularity is already cut off as the damping sum in (8) remains finite for all correlation exponents in the physically interesting interval $0<\alpha<2$.

In the regime $\alpha<1$, the FLE describes a peculiar physical situation. As the noise is anti-persistent, the fluctuation-dissipation theorem implies that the damping kernel ${\mathcal{K}}_{n-m}$ has negative values for $n\neq m$. Instead of damping, these terms thus provide anti-damping, i.e., a positive feedback for the velocity. Consequently, the mean-square displacement in the free unconfined case grows super-diffusively, $\langle x^{2}\rangle\sim t^{2-\alpha}$, for $\alpha<1$.

To investigate the case of anti-persistent noise, we perform simulations of the discretized FLE given by (8), (9) on the semi-infinite interval $(0,\infty)$ for $\alpha=0.8$ and $\alpha=0.5$. We find that the mean-square displacement and the mean-square velocity show the same qualitative behavior as in the free, unconfined case. This means the mean-square velocity quickly settles on the value $\langle v^{2}\rangle=1$ prescribed by the fluctuation-dissipation theorem, and the mean-square displacement grows as $\langle x^{2}\rangle\sim t^{2-\alpha}$.

The (scaled) probability density of particles starting at the origin at time $t=0$ is shown in Fig. 9 for $\alpha=0.5$.

As in the case of persistent noise, the probability densities taken at different times collapse when $x$ is scaled with the root-mean-square displacement $\sigma(t)$. This confirms that the system has reached the asymptotic (scaling) regime at these times. The functional form of $P(x,t)$ is again non-Gaussian but in contrast to the case of persistent noise discussed in Sec. IV, particles are depleted in the region close to the wall. We observe a similar depletion for $\alpha=0.8$.

This depletion of particles close to the reflecting wall for the FLE with anti-persistent noise is analogous to the behavior of FBM with anti-persistent noise Wada and Vojta (2018); Wada et al. (2019). However, for FBM, the probability density actually vanishes right at the wall whereas Fig. 9 shows that it remains nonzero for the FLE. As anti-persistent noise leads to sub-diffusive FBM but super-diffusive motion in the FLE, these results provide further support for the notion that the sign of the correlation determines whether or not particles are accumulated or depleted and not the type of anomalous diffusion (superdiffusion vs. subdiffusion).

V.2 Generalized Langevin equation without fluctuation-dissipation theorem

In Sec. III, we found that the probability density of the FLE on a finite interval reaches a uniform distribution for long times if the fluctuation-dissipation theorem is fulfilled. In contrast, the stationary distribution of FBM on a finite interval is non-uniform Guggenberger et al. (2019). To further test the role that the fluctuation-dissipation theorem plays in these different situations, we perform simulations of a generalized Langevin equation (2) that violates the fluctuation-dissipation relation (3). Specifically, we employ the same long-time correlated fractional Gaussian noise (4) as before but combine it with an instantaneous damping force $-\bar{\gamma}dx(t)/dt$. (This corresponds to a $\delta$ function kernel, ${\mathcal{K}}(t-t^{\prime})=\delta(t-t^{\prime})$, in the damping integral in (2).

Figure 10 presents a comparison of the stationary probability density of this generalized Langevin equation with the uniform distribution obtained in Sec. III for the FLE that fulfills the fluctuation-dissipation theorem ($\alpha=1.2$).

The figure demonstrates that the stationary probability density strongly deviates from a uniform distribution if the fluctuation dissipation theorem is not fulfilled. Instead, it develops singularities close to the walls and resembles the probability density of reflected FBM found in Ref. Guggenberger et al. (2019). We obtained similar results for $\alpha=1.5$.

In fact, FBM can be understood as the overdamped limit of the generalized Langevin equation (2) with long-time correlated random forces but instantaneous damping. In the overdamped limit, the acceleration term can be neglected, and the Langevin equation turns into $\bar{\gamma}dx(t)/dt=\xi(t)$ which is equivalent to FBM.

VI Conclusions

To summarize, motivated by recent observations of a non-Gaussian probability density for reflected FBM Wada and Vojta (2018); Wada et al. (2019) as well as corresponding deviations from a uniform stationary probability density for FBM on a finite interval Guggenberger et al. (2019), we studied FLE motion confined by reflecting potential barriers. The main part of our work focuses on the case when the random and damping forces fulfill the fluctuation-dissipation theorem.

If confined to a finite interval, the FLE therefore reaches thermal equilibrium in the long-time limit. Correspondingly, our simulations yield a stationary probability density $P(x)$ that is uniform on the interval, independent of the value of the correlation exponent $\alpha$. This is strikingly different from the behavior of FBM (which contains exactly the same type of long-range correlated noise) on a finite interval. For FBM, the stationary probability density is non-uniform and depends on $\alpha$ Guggenberger et al. (2019).

If the FLE is confined to a semi-infinite interval by a single reflecting wall, it will not reach an equilibrium state because the probability density broadens without limit (while the mean-square velocity settles on the thermal equilibrium value). Here, the probability density shows unexpected deviations from the Gaussian form observed for normal Brownian motion. Specifically, particles accumulate at the reflecting wall for persistent noise ($1<\alpha<2$) while they are depleted close to the wall for anti-persistent noise ($\alpha<1$). The probability density fulfills the scaling form $P(x,t)=Y\left[x/\sigma(t)\right]/\sigma(t)$ where $\sigma(t)$ is the root-mean-square displacement. This demonstrates that the non-Gaussian behavior is not a finite-time effect but a feature of the asymptotic long-time behavior.

In other words, the probability density of reflected FBM Wada and Vojta (2018); Guggenberger et al. (2019) shows the same type of singular behavior on a finite interval and on a semi-infinite interval, whereas the FLE shows the usual uniform probability density on a finite interval but unconventional behavior on a semi-infinite interval.

Let us contrast the behaviors of FBM and FLE motion on a semi-infinite interval. For both processes, the probability density $P(x,t)$ shows deviations from the Gaussian behavior observed for normal diffusion. Moreover, for both processes, persistent noise causes an accumulation of particles at the wall whereas anti-persistent noise leads to a depletion. This happens even though persistent noise leads to superdiffusive FBM but subdiffusive motion in the FLE (and vice verse for anti-persistent noise). This suggests that the behavior of $P(x,t)$ at the wall is controlled by the character of the noise rather than the type of anomalous diffusion. However, the accumulation and depletion effects are qualitatively stronger for FBM (where they cause the scaled probability density $\sigma(t)P(x,t)$ to either vanish or diverge at the wall) than for the FLE (where $\sigma(t)P(x,t)$ remains finite at the wall).

Our simulations have focused on the unbiased case for which the random force has zero mean, $\langle\xi(t)\rangle=0$. However, the expected behavior in the biased case, $\langle\xi(t)\rangle\neq 0$, is easily discussed. For definiteness, consider the semi-infinite interval $(0,\infty)$. A bias towards the wall, $\langle\xi(t)\rangle<0$ corresponds to a potential energy $V(x)$ that increases linearly with the distance from the wall. In this case, the system is expected to reach thermal equilibrium in the long-time limit. The stationary probability density is thus given by the equilibrium density $P(x)\sim\exp[-V(x)/(k_{B}T)]$ independent of the value of $\alpha$. If the bias is away from the wall, $\langle\xi(t)\rangle>0$, the peak of the probability density will move away from the wall with constant speed whereas its width increases more slowly. For long times, we will therefore recover the Gaussian probability density of the unconfined FLE. This means neither of the biased cases is expected to feature an unusual $\alpha$-dependent probability density in the long-time limit.

To clarify the role that the fluctuation-dissipation theorem plays in establishing the probability density, we also performed simulations for a generalized Langevin equation for which the random force $\xi(t)$ is a (long-time correlated) fractional Gaussian noise but the damping force has no memory and is proportional to the instantaneous velocity. For this equation, which violates the fluctuation-dissipation theorem, our simulations show that the stationary probability density on a finite interval is not uniform but resembles the $\alpha$ dependent probability density observed for FBM.

Recently, Holmes Holmes (2019) performed simulations of an overdamped FLE on a finite interval. In contrast to our results and in violation of the Boltzmann distribution expected in thermal equilibrium, he reports that the steady-state density is not uniform but develops depletion zones close to the boundaries. We believe that this discrepancy may stem from the fact that in Ref. Holmes (2019) the simple reflection condition $x_{n}\to 2x_{0}-x_{n},v_{n}\to-v_{n}$ is used instead of reflecting potentials to implement the walls. This effect, especially the role of the simulation time step, is discussed in more detail in Appendix A.

We also note that free FBM and FLE motion are ergodic processes Deng and Barkai (2009); Schwarzl et al. (2017). FBM and FLE motion were, however, shown to exhibit transient non-ergodic behavior in a harmonic external potential Jeon and Metzler (2012); Jeon et al. (2013) as well as transient ageing Kursawe et al. (2013).

Our results show that the probability density functions of FLE motion in a semi-infinite domain, as well as FBM in confinement and a semi-infinite domain have non-zero slopes at reflecting boundaries. This is a strong indication for why the traditional method of images Redner (2001) fails for processes fueled by correlated fractional Gaussian noise. Indeed, for FBM, the first passage behavior is only known analytically in the long time limit in a semi-infinite domain Molchan (1999), in terms of perturbation theory Wiese (2019), or through conjectures Jeon et al. (2011b); Sanders and Ambjörnsson (2012).

It will therefore be interesting to compare the reflecting walls considered in this paper with absorbing walls. Absorbing walls prevent the FLE from reaching a thermal equilibrium state, even on a finite interval. The probability density is thus not simply given by the appropriate Boltzmann distribution. Instead, similar to the results of Sec. IV, the probability density may be expected to show nontrivial $\alpha$-dependent behavior. Note that absorbing walls in FBM are known to produce a power-law singularity of the probability density, $P(x,t)\sim x^{\kappa}$ similar to reflecting walls but with an exponent $\kappa=2/\alpha-1$ Chatelain et al. (2008); Zoia et al. (2009); Wiese et al. (2011) (whereas the exponent takes the value $\kappa=2/\alpha-2$ at a reflecting wall).

Both the FLE and FBM are used to model a wide variety of anomalous diffusion processes in complex systems ranging from electronic networks to the motion inside biological cells. Our results indicate that thermal equilibrium plays a crucial role in the accumulation or depletion of particles due to the interplay of the long-time memory and the reflecting walls. In a true equilibrium state, no such accumulation/depletion zones occur. Instead the probability density follows the Boltzmann distribution, as expected. If, on the other hand, the system violates the fluctuation-dissipation relation (as is the case for FBM or the generalized Langevin equation of Sec. V.2), the probability density develops power-law singularities close to the walls and becomes $\alpha$-dependent. Our simulations also show that correlation-driven accumulation or depletion of particles can occur even in systems that fulfill the fluctuation dissipation theorem if they are not in an equilibrium state (as is the case on a semi-infinite interval). It is interesting to consider this question for biological systems. Some anomalous diffusion processes in cells are likely well described by thermal equilibrium, others may be dominated by the active motion of some of the constituents. Working out in detail how this affects the probability densities remains a task for the future.

Acknowledgements.

This work was supported in part by the NSF under Grant Nos. DMR-1506152 and DMR-1828489 (T.V.) and by the DFG under Grant No. ME 1535/7-1 (R.M.). R.M. also acknowledges support from the Foundation of Polish Science (Fundacja na rzecz Nauki Polskiej) in terms of an Alexander von Humboldt Polish Honorary Research Scholarship.

Appendix A Alternative implementation of reflecting walls

In the main part of this paper, the reflecting walls have been modeled as “soft” repulsive potentials (6), $V(x)=V_{0}\exp[\mp\lambda(x-x_{0})]$, as visualized in Fig. 11.

Our simulations in Sec. III demonstrated that particles governed by the fractional Langevin equation and confined to a finite interval by such potentials reach a true thermodynamic equilibrium state. Their kinetic energy fulfills the equipartition theorem and the spatial probability density follows the Boltzmann distribution, leading to a uniform stationary $P(x)$.

With increasing $\lambda$ the repulsive potential becomes steeper, reaching an infinitely high “hard” wall at position $x_{0}$ in the limit of $\lambda\to\infty$. Such a wall is expected to reflect the particle elastically, i.e., the velocity simply changes sign when the particle hits the wall. In the simulations, it is thus tempting to replace the external force $F(x_{n})$ in the recursion relations (8) and (9) by the simple reflection condition

[TABLE]

if the particle finds itself behind a wall located at $x_{0}$.

Even though this reflection step cannot be interpreted as a normal Euler-like time-step in the numerical approximation of the fractional Langevin equation (because the change of velocity is not small), one might still expect it to capture the relevant physics. In fact, for the normal Langevin equation (1), this implementation of the reflecting walls leads to the expected behavior.

To test the effects of such hard reflecting boundary conditions on the fractional Langevin equation, we have analyzed the time evolution of the probability density $P(x,t)$ on a finite interval $(-L,L)$, combining the recursion relations (8) and (9) with the reflection condition (13) for the two walls at $x=\pm L$. The results for $L=15$ and a correlation exponent $\alpha=1.5$ are presented in Fig. 12.

The figure shows that the probability density reaches a stationary state for long times, but the stationary $P(x)$ is not uniform, in contrast to what is expected in thermal equilibrium. For comparison, the corresponding system with reflecting potentials (6) reaches a uniform stationary density, as shown in Fig. 3.

The data in Fig. 12 were obtained using the same time step $\Delta t=0.01$ that was used successfully in the main part of the paper. To gain additional insight, we have performed simulations of the FLE with reflection condition (13) with several smaller time steps. The resulting probability densities at a fixed time $t=167772$ are presented in Fig. 13.

The figure shows that the probability density continues to change with $\Delta t$ down to the smallest studied value of $\Delta t=3.125\times 10^{-4}$. This needs to be contrasted with the results for the reflecting potential shown in Fig. 1 which are essentially independent of the time step for $\Delta t\lesssim 0.01$. The deviations in Fig. 13 from a uniform distribution decrease with decreasing $\Delta t$. The inset of Fig. 13 shows that the probability density $P(0,t)$ in the center of the interval and $P(-15,t)$ at the wall vary linearly with $(\Delta t)^{1/2}$. We have used this relation to extrapolate $P(0,t)$ and $P(-15,t)$ to zero time step. Both probability densities extrapolate to the same value which is very close to 0.0333 expected for a uniform distribution. This suggests that the system actually reaches the equilibrium state in the limit of $\Delta t\to 0$ even if the reflection condition (13) is used. However, extremely small $\Delta t$ appear to be necessary to approach that limit. This would significantly reduce practical applications.

We have also performed simulations on a semi-infinite interval using the reflection condition (13) and observed similar deviations from the expected behavior. A detailed microscopic analysis of this phenomenon will be presented elsewhere 222T. Vojta, S. Skinner, and R. Metzler, unpublished..

Appendix B Evaluation of the damping integral

A straight-forward implementation of the recursion relations (8) and (9) is not very efficient numerically, as the total computational effort scales quadratically with the number $N_{t}$ of time steps. This is caused by the damping term,

[TABLE]

which requires summing over all time steps preceding the current one in each evaluation of the damping force. The unfavorable quadratic scaling severely limits the maximum possible simulation times. (Note that the creation of the correlated random numbers via the Fourier filtering method shows a much more favorable $N_{t}\log N_{t}$ scaling with the number of time steps.)

We have developed an improved method that speeds up the evaluation of the damping sums by several orders of magnitude. It is based on the observation that, for large time lag $(n-m)$, the kernel ${\mathcal{K}}_{n-m}$ is small and slowly varying. This allows us to approximate it via Taylor expansion. Consider a subsequence of $2j+1$ consecutive terms in the damping sum, centered at term $i$,

[TABLE]

We now Taylor-expand the kernel about $m=i$ giving

[TABLE]

Inserting this into the partial sum (15) yields

[TABLE]

The key insight is that the kernel ${\mathcal{K}}$ and its derivatives can be precalculated, and the sums in the expression (17) need to be computed only once and can then be used in every damping sum at later times. This reduces the numerical effort for computing the partial sum $I_{n}(i,j)$ from $2j+1$ multiplications and additions to only one multiplication and addition if only the leading term in the expansion is kept, or to two, or perhaps three multiplications and additions if higher order terms are kept. The accuracy of this approximation can be easily controlled by varying the number of terms kept in the Taylor expansion and the width $2j+1$ of the interval used in each partial sum.

In our simulations, we vary the interval width with increasing lag time $n-m$. Specifically, when calculating the damping sum $S_{n}$, we calculate the first $M_{0}$ terms ($n-m$ from 1 to $M_{0}$) exactly. We then use intervals of width $w_{1}$ until $n-m$ reaches $M_{1}$, then intervals of width $w_{2}$ until $n-m$ reaches $M_{2}$, and so on. Both the widths $w_{k}$ and the cutoffs $M_{k}$ increase in a geometric fashion. For our largest simulations involving up to $2^{27}\approx 134$ million time steps, we use 4 different interval widths. (This effectively reduces the scaling of the numerical effort from $N_{t}^{2}$ to $N_{t}^{1.2}$, albeit with a large prefactor.) Careful optimization of the interval widths $w_{k}$ and cutoff points $M_{k}$ keeps the relative approximation error for position and velocity at about $10^{-5}$ while improving the performance by several orders of magnitude. We have performed systematic benchmarks for a simulation with $2^{21}$ time steps. On a standard core-i5 PC, the calculation of a single trajectory takes about 1800 seconds if the damping sums are evaluated term by term. Our fully optimized version takes only 4 seconds for the same calculation. For larger numbers of time steps the gains are even bigger.