Ergodic properties of heterogeneous diffusion processes in a potential well

TL;DR

This paper studies the ergodic and non-ergodic behaviors of heterogeneous diffusion processes in potential wells, revealing how the interplay between diffusivity and potential shape determines long-term statistical properties.

Contribution

It provides a comprehensive analysis of ergodic properties in space-dependent diffusion processes with power-law potentials and diffusivities, including new classifications based on parameter regimes.

Findings

System is ergodic when potential dominates diffusivity ($eta>\alpha$).

System is non-ergodic when diffusivity dominates potential ($eta<\alpha$).

The $eta=\alpha$ case shows delicate ergodic behavior depending on parameters.

Abstract

Heterogeneous diffusion processes can be well described by an overdamped Langevin equation with space-dependent diffusivity $D(x)$. We investigate the ergodic and non-ergodic behavior of these processes in an arbitrary potential well $U(x)$ in terms of the observable---occupation time. Since our main concern is the large-$x$ behavior for long times, the diffusivity and potential are, respectively, assumed as the power-law forms $D(x)=D_0|x|^\alpha$ and $U(x)=U_0|x|^\beta$ for simplicity. Based on the competition roles played by $D(x)$ and $U(x)$, three different cases, $\beta>\alpha$, $\beta=\alpha$, and $\beta<\alpha$, are discussed. The system is ergodic for the first case $\beta>\alpha$, where the time average agrees with the ensemble average, being both determined by the steady solution for long times. In contrast, the system is non-ergodic for $\beta<\alpha$, where the relation between time average and ensemble average is uncovered by infinite-ergodic theory. For the middle case $\beta=\alpha$, the ergodic property, depending on the prefactors $D_0$ and $U_0$, becomes more delicate. The probability density distribution of the time averaged occupation time for three different cases are also evaluated from Monte Carlo simulations.

Equations

$$⟨x^2(t)⟩\simeq 2K_{\gamma }{t}^{\gamma },$$

$$\frac{dx}{dt}=\sqrt{2D(x)}\xi (t),$$

$$D(x)=D_0|x{|}^{\alpha },\alpha \le 2,$$

$$U(x)=U_0|x{|}^{\beta }.$$

$$\frac{dx}{dt}=\sqrt{2D(x)}\xi (t),$$

$$\begin{array}{rl}y& ={\int }^x\frac{1}{\sqrt{2D(x^{\prime })}}d{x}^{\prime }\\ & =\frac{\sqrt{2/D_0}}{2-\alpha }|x{|}^{\frac{2-\alpha }{2}}\text{sign}(x),\end{array}$$

$$p(x,t)=\frac{|x{|}^{-\alpha /2}}{\sqrt{4\pi D_{0}t}}\exp \left(-\frac{|x{|}^{2-\alpha }}{(2-\alpha {)}^2{D}_{0}t}\right)$$

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

$$y(t)=(2D_0{)}^{-1/2}\text{sign}(x)\ln |x|.$$

$$p(x,t)=\frac{1}{\sqrt{4\pi D_{0}t}}|x{|}^{-1-\frac{\ln |x|}{4D_{0}t}}$$

$$⟨x^2(t)⟩\propto e^{4D_{0}t},$$

$$\frac{dx}{dt}=f(x)+\sqrt{2D(x)}\xi (t),$$

$$\frac{\partial p(x,t)}{\partial t}={\mathcal{L}}_{\text{FP}}p(x,t)$$

$${\mathcal{L}}_{\text{FP}}=-\frac{\partial }{\partial x}f(x)+\frac{\partial }{\partial x}\sqrt{D(x)}\frac{\partial }{\partial x}\sqrt{D(x)}.$$

$$J(x,t)=f(x)p(x,t)-\sqrt{D(x)}\frac{\partial }{\partial x}\sqrt{D(x)}p(x,t).$$

$$p^{st}(x)=\frac{N}{\sqrt{D(x)}}\exp \left({\int }^x\frac{f(x^{\prime })}{D(x^{\prime })}d{x}^{\prime }\right),$$

$$\begin{array}{r}{p}^{st}(x)=\{\begin{array}{ll}\frac{N}{\sqrt{D_0}}|x{|}^{-\alpha /2}\exp \left(-\frac{U_0\beta }{D_0(\beta -\alpha )}|x{|}^{\beta -\alpha }\right)& \beta \ne \alpha \\ \frac{N}{\sqrt{D_0}}|x{|}^{-\alpha /2-U_0\beta /D_0}& \beta =\alpha .\end{array}\end{array}$$

$$N=\frac{\sqrt{D_0}(\beta -\alpha )}{2\Gamma (\varrho )}{\left(\frac{U_0\beta }{D_0(\beta -\alpha )}\right)}^{\varrho }$$

$${\bar{p}}_{\mathcal{D}}=\frac{T_{\mathcal{D}}(t)}{t}$$

$$f({\bar{p}}_{\mathcal{D}})=\delta ({\bar{p}}_{\mathcal{D}}-P_{\mathcal{D}}^{st}),$$

$$\begin{array}{rl}& f_{\text{CT}}({\bar{p}}_{\mathcal{D}})=\frac{\sin (\mu \pi )}{\pi }\times \\ & \frac{{\mathcal{R}}_{\mathcal{D}}{\bar{p}}_{\mathcal{D}}^{\mu -1}(1-{\bar{p}}_{\mathcal{D}}{)}^{\mu -1}}{{\mathcal{R}}_{\mathcal{D}}^2(1-{\bar{p}}_{\mathcal{D}}{)}^{2\mu }+{\bar{p}}_{\mathcal{D}}^{2\mu }+2{\mathcal{R}}_{\mathcal{D}}(1-{\bar{p}}_{\mathcal{D}}{)}^{\mu }{\bar{p}}_{\mathcal{D}}^{\mu }\cos (\mu \pi )}\end{array}$$

$${\mathcal{R}}_{\mathcal{D}}=\frac{P_{\mathcal{D}}^{st}}{1-P_{\mathcal{D}}^{st}}.$$

$${\hat{\psi }}_{\pm }(\lambda )={\int }_0^{\infty }{e}^{-\lambda t}{\psi }_{\pm }(t)dt.$$

$${\mathcal{R}}_{+}=\frac{A_{+}}{A_{-}}$$

$$f({\bar{p}}_{+})=\delta \left({\bar{p}}_{+}-\frac{{\mathcal{R}}_{+}}{1+{\mathcal{R}}_{+}}\right),$$

$$\frac{\partial P(x_0,x,t)}{\partial t}={\mathcal{L}}_{\text{FP}}^{\ast }P(x_0,x,t),$$

$${\mathcal{L}}_{\text{FP}}^{\ast }=f(x_0)\frac{\partial }{\partial x_0}+\sqrt{D(x_0)}\frac{\partial }{\partial x_0}\sqrt{D(x_0)}\frac{\partial }{\partial x_0}.$$

$$\text{Prob}\{t\ge {\tau }_{+}\}=S(x_0,{\tau }_{+})={\int }_0^{\infty }P(x_0,x,{\tau }_{+})dx,$$

$$Q(x_0,{\tau }_{+})=-\frac{\partial }{\partial {\tau }_{+}}S(x_0,{\tau }_{+});$$

$$\frac{\partial Q(x_0,{\tau }_{+})}{\partial {\tau }_{+}}={\mathcal{L}}_{\text{FP}}^{\ast }Q(x_0,{\tau }_{+})$$

Full text

††preprint: AIP/123-QED

I Introduction

The anomalous diffusion phenomena are ubiquitous in the natural world, especially in numerous microscopic systems. Over the last two decades, much effort has been devoted to the study of anomalous diffusion phenomena Haus and Kehr (1987); Bouchaud and Georges (1990); Metzler and Klafter (2000), characterized by nonlinear time dependence of mean squared displacement (MSD)

[TABLE]

where $\gamma\neq 1$ is the anomalous diffusion exponent and $K_{\gamma}$ is the anomalous diffusion coefficient with physical dimension $[\textrm{cm}^{2}/\textrm{s}^{\gamma}]$; it is called subdiffusion for $0<\gamma<1$ and superdiffusion for $\gamma>1$.

However, the MSD is not the only criterion to distinguish the processes, i.e., some fundamentally different processes may have the same MSD. For example, the subdiffusion was first observed in amorphous semiconductors described by continuous-time random walk (CTRW) models Scher and Montroll (1975), being Non-Gaussian process, while it was also observed in the class of viscoelastic system described by the generalized Langevin equation with (tempered) power-law friction kernel Lutz (2001); Goychuk (2012); Ślȩzak, Metzler, and Magdziarz (2018); M.-Garcia et al. (2018); Deng and Barkai (2009) and of (tempered) fractional Brownian motion Mandelbrot and Ness (1968); Meerschaert and Sabzikar (2013); Chen, Wang, and Deng (2017), both of which are Gaussian processes. There are also many different superdiffusion processes, such as, Lévy flight Shlesinger, Zaslavsky, and Frisch (1995), Lévy walk Zaburdaev, Denisov, and Klafter (2015), and those in the underdamped Langevin systems Eule et al. (2007); Wang, Chen, and Deng (2018a). In recent years, more and more dynamics in heterogeneous systems can be well described by the overdamped Langevin equation with space-dependent diffusivity Cherstvy, Chechkin, and Metzler (2013); Cherstvy and Metzler (2013, 2014)

[TABLE]

such as the Richardson diffusion in turbulence Richardson (1926), mesoscopic approaches to transport in heterogenous porous media Haggerty and Gorelick (1995); Dentz et al. (2012) and on random fractals Loverdo et al. (2009); O’Shaughnessy and Procaccia (1985). The heterogeneous dynamical behavior is particularly remarkable in biological systems. The cytoplasm of biological cells is always crowded with various obstacles, including proteins, nucleic acids, ribosomes, the cytoskeleton, as well as internal membranes compartmentalizing the cell Zimmerman and Minton (1993); Zhou, Rivas, and Minton (2008). The nonuniform distribution of crowders in the cytoplasm provides the heterogenous media for tracer particles of different sizes in it. The motion of the tracer particles might be imposed by the cage effects if the size of the particle is comparable to the local mesh size in the crowded environments Wong et al. (2004); Godec, Bauer, and Metzler (2014); Lee et al. (2014). This kind of heterogeneous diffusion process also can be realized in the experiments on eukaryotic cells Kühn et al. (2011) and a local variation of the temperature in thermophoresis experiments Mast et al. (2013); Maeda, Tlusty, and Libchaber (2012).

Nowadays, single particle tracking has become a powerful tool to study transport processes in cellular membranes Saxton and Jacobson (1997) and probe the microrheology of the cytoplasm Wirtz (2009); Yao et al. (2009). It can be used to evaluate the time averaged observables in cells through video microscopy of fluorescently labeled molecules. The time average is equal to the ensemble average for an ergodic system. But for anomalous diffusion, especially for the molecules diffusing in living cells, the time average often becomes a random variable and irreproducible for the individual trajectories. For the heterogeneous diffusion process described by Langevin equation (2) with the diffusivity $D(x)$ in the power-law form Gardiner (1983); Srokowski (2007):

[TABLE]

it is subdiffusion for $\alpha<0$ and superdiffusion for $\alpha>0$; the weak non-ergodic behavior has also been observed in Cherstvy, Chechkin, and Metzler (2013); Cherstvy and Metzler (2013, 2014), where the time averaged MSD remains linear in the lag time and thus differs from the corresponding ensemble average.

Here we consider ergodic/non-ergodic properties of the heterogeneous diffusion processes (2) in a potential well

[TABLE]

The non-ergodic behavior of subdiffusive CTRW model in a confined potential has been investigated in Bel and Barkai (2006). Even for the ergodic system where normal (anomalous) diffusion modeled by (fractional) Brownian motion and the (fractional) Langevin equation is confined by a harmonic potential, the time averaged MSDs behave differently from the ensemble averaged counterparts Jeon and Metzler (2012). We here investigate the ergodic behavior in terms of the occupation time $T_{\mathcal{D}}(t)$ for a particle staying in a certain domain $\mathcal{D}$. It is found that the physical mechanism of the heterogeneous diffusion process in (2) is different from CTRW model and fractional Langevin equation, although they may have the same anomalous diffusion exponents defined in (1). If the potential well (4) is deep enough, the system might reach a steady state and the ergodic hypothesis is shown to be valid.

However, for a potential well that is not very deep, the steady solution is found to be non-normalizable, implying that the non-trival steady solution may not exist. Non-normalizable steady solutions have been previously discussed in the literature on physical issues, e.g., Lévy walk and laser-cooled atoms Kessler and Barkai (2010); Lutz and Renzoni (2013); Rebenshtok et al. (2014); Aghion, Kessler, and Barkai (2017), the weakly-chaotic maps Korabel and Barkai (2009); Akimoto and Miyaguchi (2010), nonlinear oscillators Meyer and Kantz (2017), the thermodynamics Aghion, Kessler, and Barkai (2018), and multiplicative noise Leibovich and Barkai (2018). From these discussions, one can note that the (non-normalizable) steady solution is determined by the competition role played by the potential $U(x)$ and the space-dependent diffusivity $D(x)$ in our concerned model. A generalized form of ergodicity, namely infinite-ergodic theory Aaronson (1997); J. Aaronson and Zweimuller (2005); Thaler and Zweimüller (2006), applies to the non-normalizable case, which uncovers the relation between time averages and ensemble averages and helps to characterize the distribution of time averages.

This paper is organized as follows. In Sec. II, we briefly review the heterogeneous diffusion processes, especially their probability density functions (PDFs), and then turn to their steady solutions under an arbitrary potential well. In Sec. III, the system with a deep potential (i.e., $\beta>\alpha$) is explicitly discussed; in this case, the steady solution is normalizable, and the distribution of time averaged occupation time is a $\delta$-function centered on the value of the ensemble average, which is determined by the steady solution for long times. In the case of comparable potential well (i.e., $\beta=\alpha$) and weak potential well (i.e., $\beta<\alpha$), where the steady solution is non-normalizable, the infinite-ergodic theory and the distribution of the time averages are discussed in Sec. IV and Sec. V, respectively. A summary of the key results is provided in Sec. VI. In the appendices some mathematical details are collected.

II Heterogeneous diffusion model

II.1 No potential well

We firstly review the stochastic Langevin model in a heterogeneous medium with space-dependent diffusivity $D(x)$, namely,

[TABLE]

where $\xi(t)$ is a Gaussian white noise with zero mean $\langle\xi(t)\rangle$=0 and $\delta$-correlation $\langle\xi(t)\xi(t^{\prime})\rangle=\delta(t-t^{\prime})$. Here the multiplicative noise term $\sqrt{2D(x)}$ is interpreted in the Stratonovich sense Risken (1989). For a general stochastic differential equation, the Stratonovich prescription has the advantage of ordinary chain rule formulas under a transformation Risken (1989); Øksendal (2005) and of yielding physically correct results especially for the noise with infinitely short correlation times West et al. (1979), while the Itô prescription is mathematically and technically the most satisfactory Gardiner (1983) and easy for simulations.

During the analytical calculations, the specific form of $D(x)$ in (3) is taken. It is actually regularized at $x=0$ in physical systems and in simulations Cherstvy, Chechkin, and Metzler (2013). For $\alpha>0$, a modified form $D(x)=D_{0}(A_{0}+|x|^{\alpha})$ with a positive constant $A_{0}$ is to prevent the particles trapping at $x=0$. For $\alpha<0$, we take $D(x)=\frac{D_{0}A_{0}}{A_{0}+|x|^{-\alpha}}$ to avoid the divergences at $x=0$. On the other hand, the regularization at $x=0$ makes it satisfy the Lipschitz condition and the request $\alpha\leq 2$ ensures the growth condition for existence and uniqueness of the solution of the Markovian stochastic differential equation (5) Øksendal (2005). By the substitution Cherstvy, Chechkin, and Metzler (2013)

[TABLE]

where $y(t)$ is the standard Brownian motion Risken (1989), the PDF and ensemble averaged MSD of $x(t)$ for $\alpha<2$ are obtained in Cherstvy, Chechkin, and Metzler (2013)

[TABLE]

and

[TABLE]

which is obviously superdiffusion for $\alpha>0$ and subdiffusion for $\alpha<0$. For the critical value $\alpha=2$, by the same substitution (6), one can obtain the relation between $y(t)$ and $x(t)$:

[TABLE]

Therefore, the PDF and ensemble averaged MSD of $x(t)$ for $\alpha=2$ are, respectively,

[TABLE]

and

[TABLE]

which grows exponentially and much faster than (8). The case of $\alpha=2$ is common in applications; more often it emerges together with a harmonic potential. In this case, the Langevin system can be interpreted as the model for a massless particle connected to a spring with the spring coefficient varying randomly in time Deutsch (1994), which also is used to characterize the Lévy-type stochastic processes in terms of continuous trajectories of walker motion Lubashevsky, Friedrich, and Heuer (2009). Note that the PDFs (7) and (10) decay faster than any power-law form, implying that the condition $\alpha\leq 2$ indeed guarantees the existence of the solution of Langevin equation (5).

II.2 Steady solution in a potential well

We assume that the potential well $U(x)=U_{0}|x|^{\beta}$ with $U_{0}>0$, in which the exponent $\beta$ is usually positive; especially $\beta=2$ is for the potential between the monomers in a bead spring polymer chain Shin et al. (2015). Sometimes, the exponent $\beta$ might also be negative (e.g., the Lennard-Jones potential Shin et al. (2015); Shin, Cherstvy, and Metzler (2015) approximating the interaction between a pair of neutral atoms or molecules) and in this case we regularize $U(x)$ at the origin by the way similar to $D(x)$. Then the force is $f(x)=-dU(x)/dx=-U_{0}\beta|x|^{\beta-1}\textrm{sign}(x)$ and the Langevin equation for a particle in this potential well is presented as

[TABLE]

the Fokker-Planck equation Risken (1989); Cherstvy, Chechkin, and Metzler (2013) of which for the PDF $p(x,t)$ is

[TABLE]

with the Fokker-Planck operator

[TABLE]

This equation can be represented through probability current $J(x,t)$ due to mass conservation Risken (1989) as $\partial_{t}p(x,t)=-\partial_{x}J(x,t)$ with

[TABLE]

In many circumstances, the steady solution $p^{st}(x)$ of (12) can be directly obtained by setting $J(x,t)=0$, which yields

[TABLE]

with $N$ being the normalization constant. It can be seen that the strengths of $f(x)$ and $D(x)$ dominate the shape of $p^{st}(x)$ in (16). For different exponents of the power-law form of the potential $U(x)$ and diffusivity $D(x)$, the solution $p^{st}(x)$ behaves as

[TABLE]

Here it should be emphasized that the real diffusivity $D(x)$ in (3) is regular at origin. So the PDF $p^{st}(x)$ should also be regular at origin. What we pay attention to in this paper is the large-$x$ behavior of the solution $p^{st}(x)$ in (17).

For the case of $\beta>\alpha$ and $\beta>0$, $p^{st}(x)$ decays exponentially for large $x$ and it can be normalized whatever other parameters are. After some calculations, the normalization constant $N$ can be obtained as

[TABLE]

with $\varrho=(2-\alpha)/(2(\beta-\alpha))$. For the other two cases $\beta=\alpha$ and $\beta<\alpha$, however, when $\alpha>0$, the PDF $p^{st}(x)$ decays as power-law $|x|^{-\alpha/2-U_{0}\beta/D_{0}}$ and $|x|^{-\alpha/2}$, respectively, which even can not be normalized for some parameters. The divergence at infinity means that the steady state may not exist in the usual sense and a time-dependent solution should be considered to measure some observables.

In the following, we will use the steady solution (17) to investigate the time and ensemble averages of some observables as well as the ergodic properties of the Langevin system (12). It has been shown that the steady solution is completely different for different magnitudes of $\beta$ and $\alpha$. For the first case of $\beta>\alpha$ and $\beta>0$, the steady solution can be normalized even for $\alpha<0$. For the other two cases $\beta=\alpha$ and $\beta<\alpha$, the condition $\alpha>0$ is assumed; otherwise, the steady solution grows with time for large $x$, which is impossible for a PDF. We mainly explore the occupation time statistic $T_{\mathcal{D}}(t)$ in a given domain $\mathcal{D}$ for three different cases, respectively, and then make some comparisons. The concept of occupation time has been studied by probabilists for a long time Lamperti (1958); in recent years, it attracts interests of physicists together with some related observables, such as first-passage time Redner (2001), persistence Majumdar, Rosso, and Zoia (2010), and the area under Brownian and non-Brownian path Grebenkov (2007); Majumdar and Comtet (2005); Barkai, Aghion, and Kessler (2014).

III Ergodic phase for $\beta>\alpha$

In the case of $\beta>\alpha$, we claim that for large $t$ the PDF of the occupation fraction

[TABLE]

satisfying

[TABLE]

where $P^{st}_{\mathcal{D}}=\int_{\mathcal{D}}p^{st}(x)dx$ is the ensemble average in the steady state; this kind of phenomenon has been ever detected in Bel and Barkai (2006), where the PDF of $\bar{p}_{\mathcal{D}}$ for a certain CTRW model is

[TABLE]

with $\mu$ being the exponent of the power-law distributed waiting times and

[TABLE]

When $\mu=1$, the PDF of $\bar{p}_{\mathcal{D}}$ in (20) reduces to the ergodic phase (19).

Now we want to derive the PDF $f(\bar{p}_{\mathcal{D}})$ in (19) for our model (in Langevin picture) based on the CTRW framework in Bel and Barkai (2006). It was proposed by Fogedby Fogedby (1994) that an overdamped Langevin equation (with additive Gaussian noise) coupled with a subordinator can describe the process defined by CTRW model, being the same one after taking scaling limit. But it seems not very clear what the multiplicative noise in our model (12) underlies in a CTRW model. Some papers Fedotov and Falconer (2012); Srokowski (2009); Kamińska and Srokowski (2018) derived the generalized Fokker-Planck equation from CTRW model with variable jumping rate, which shows that the space-dependent diffusivity may come from space-dependent distribution of waiting time. Here, we present an alternative method to clarify a clear relation between the space-dependent diffusivity in a Langevin equation and the distribution of waiting times in CTRW framework (see the detailed derivation in Appendix A). The multiplicative noise term $\sqrt{2D(x)}$ in (12) implies the exponential distributed waiting time with rate $2D(x)$ in CTRW model. The bigger $D(x)$ is, the larger number of renewals per unit time gets and the faster the diffusion becomes. Note that even the exponent of the power-law form of $D(x)$ decides the diffusion behavior of the particle in (5) (superdiffusion for $\alpha>0$ and subdiffusion for $\alpha<0$), the mechanism is, respectively, completely different from the superdiffusive Lévy flight Shlesinger, Zaslavsky, and Frisch (1995) and subdiffusive CTRW with power-law distributed waiting times Haus and Kehr (1987); Bouchaud and Georges (1990); Metzler and Klafter (2000). For the system (5), the second moment of the jump lengths and the mean of waiting times are all finite, moreover the latter is space-dependent; this interpretation from CTRW viewpoint plays an important role in deriving and understanding the result (19), which will be shown immediately.

The distribution of the occupation time of a Brownian motion or non-Brownian motion can be derived by probability theory Lévy (1939) or by backward Feynmann-Kac approach Kac (1949); Majumdar (2005); Turgeman, Carmi, and Barkai (2009); Wang, Chen, and Deng (2018b). The form of the PDF $f(\bar{p}_{\mathcal{D}})$ in (20) has been derived by Lamperti Lamperti (1958) based on the mathematical theory of occupation times and rederived by Bel and Barkai Bel and Barkai (2006) in the statistical physics community for a two-state process. Here, we also consider a two-state process with the two states: the occupation times $T_{+}$ and $T_{-}$ in the positive and negative half-space, respectively. We take advantage of the method in Bel and Barkai (2006) and find that for the PDF $f(\bar{p}_{\mathcal{D}})$, the dominant role is played by the distributions of the sojourn times in the two states, denoted as $\psi_{+}(\tau)$ and $\psi_{-}(\tau)$, respectively, for states in positive and negative half-spaces. As usual, we deal with the PDFs by using the techniques of Laplace transform, defined as

[TABLE]

According to the above discussions that the distribution of waiting times is space-dependent exponential form, we know that the mean of waiting time is finite and so is $\psi_{\pm}(t)$. Therefore, the expansions of $\psi_{\pm}(t)$ for small $\lambda$ are $\hat{\psi}_{\pm}(\lambda)=1-A_{\pm}\lambda$, where $A_{\pm}$ are their first moments, respectively. For the occupation fraction $\bar{p}_{+}=T_{+}/t$, the result in Bel and Barkai (2006) is

[TABLE]

in (20) and the PDF of $\bar{p}_{+}$ is

[TABLE]

a $\delta$-function due to the finite mean of $\hat{\psi}_{\pm}(\lambda)$ (i.e., $\mu=1$ in (20)).

III.1 Calculations of $A_{+}$ and $A_{-}$

The most important thing remaining is to calculate the two values $A_{+}$ and $A_{-}$, i.e., the means of the sojourn times $T_{\pm}$. We will focus on the state in positive half-space while the case of the negative one is similar. The sojourn time $T_{+}$ in positive half-space can be viewed as the first-passage time of a particle starting at the position $\epsilon>0$ and reaching the origin $x=0$ for the first time. We eventually obtain the mean first-passage time $A_{+}$ by taking $\epsilon\rightarrow 0$. Intuitively, the value $A_{+}$ must tend to zero when $\epsilon\rightarrow 0$ since the particle has arrived at $x=0$ at initial time. This is completely right, but it could still provide us with a correct result of $\mathcal{R}_{+}$ in (23). This technique has also been used to calculate the area under the Brownian excursion Majumdar and Comtet (2005); Majumdar (2005) and that under the Bessel excursion Barkai, Aghion, and Kessler (2014); Kessler and Barkai (2012).

Firstly, let us pay attention to the PDF $P(x_{0},x,t)$ of finding a particle at position $x$ at time $t$ starting from initial position $x_{0}$. It satisfies the backward Kolmogorov equation Arecchi, Politi, and Ulivi (1982); Risken (1989),

[TABLE]

where $\mathcal{L}_{\textrm{FP}}^{\ast}$ is the adjoint operator of $\mathcal{L}_{\textrm{FP}}$ in (13), defined as

[TABLE]

Since the first-passage time that a particle starts at $x_{0}=\epsilon$ and reaches $x=0$ for the first time is considering, for (25) we set the absorbing boundary condition at $x=0$. With this boundary condition, the survival probability $S(x_{0},\tau_{+})$, namely, the probability that a particle starts at $x_{0}$ and keeps in the positive half-space before the time $\tau_{+}$, satisfies Gardiner (1983); Deng, Wu, and Wang (2015)

[TABLE]

where $t$ is the exact time that the particle leaves the positive half-space. From the second equality in (27), we know that $S(x_{0},\tau_{+})$ also satisfies the backward Kolmogorov equation (25) with the same boundary condition. From the first equality in (27), one can obtain the relation between the PDF of first-passage time $Q(x_{0},\tau_{+})$ and the survival probability

[TABLE]

thus $Q(x_{0},\tau_{+})$ also satisfies the backward Kolmogorov equation (25), i.e.,

[TABLE]

with the normalization $\int_{0}^{\infty}Q(x_{0},\tau_{+})d\tau_{+}=1$. Multiplying $\tau_{+}$ and integrating over $\tau_{+}$ in both sides of (29), we obtain the equation for the first moment of first-passage time $\langle\tau_{+}\rangle$:

[TABLE]

The boundary condition of (30) is the same as those of $S(x_{0},\tau_{+})$ (i.e., absorbing at $x=0$), due to its relation with the survival probability $S(x_{0},\tau_{+})$:

[TABLE]

With the specified boundary condition, Eq. (30) has a unique solution Gardiner (1983); Arecchi, Politi, and Ulivi (1982)

[TABLE]

where $W(x)=\sqrt{D(x)}\exp[\int^{x}f(x^{\prime})/D(x^{\prime})]dx^{\prime}$.

Now, we come back to the initial task of calculating $A_{\pm}$: the means of the sojourn time $\psi_{\pm}(t)$. For this, taking $x_{0}=\epsilon$ in (32) and letting $\epsilon\rightarrow 0$ result in

[TABLE]

On the other hand, we calculate $A_{-}$ by assuming that the particle starts at $x_{0}=-\epsilon$ with the boundary condition absorbing at $x=0$. The mean first-passage time $\langle\tau_{-}\rangle$ for a particle reaching $x=0$ for the first time also satisfies (30) and its solution is

[TABLE]

Being similar to (33), there exists

[TABLE]

and then $\mathcal{R}_{+}$ in (23) is obtained as

[TABLE]

where the steady solution is given in (16). Though $A_{\pm}\rightarrow 0$ as $\epsilon\rightarrow 0$, the quantity of $\mathcal{R}_{+}$ of interest only depends on the speed of their tendencies to zero. Substituting this result into (24) leads to the PDF of the occupation fraction in positive half-space

[TABLE]

where $P_{+}^{st}=\int_{0}^{\infty}p^{st}(x)dx$ is the ensemble average of the occupation time in positive half-space in the steady state. This equation shows that the ergodic hypothesis is valid with respect to the occupation fraction in positive half-space.

More generally, for any domain $\mathcal{D}$, the result (19) is also valid, such as a disconnected domain $\mathcal{D}=\{|x|>a\}$. In this case, the whole domain is divided into three sub-domains $\mathcal{D}_{1}=\{x<-a\}$, $\mathcal{D}_{2}=\{-a\leq x\leq a\}$, and $\mathcal{D}_{3}=\{x>a\}$. The sojourn time in these three sub-domains are denoted as $\psi_{1}(\tau),\,\psi_{2}(\tau),\,\psi_{3}(\tau)$ with finite means $A_{1},\,A_{2},\,A_{3}$, respectively. Following the above discussions, the key is to calculate $(A_{1}+A_{3})/A_{2}$. Taking advantage of the result (36) and considering the domains $\mathcal{D}_{1}$ and $\mathcal{D}_{2}$ as a whole, we have

[TABLE]

similarly, considering the domains $\mathcal{D}_{2}$ and $\mathcal{D}_{3}$ as a whole, there exists

[TABLE]

Combining (38) and (39) results in

[TABLE]

Substituting (40) into (24) leads to the PDF of the occupation fraction $\bar{p}_{\mathcal{D}}$ in disconnected domains $\mathcal{D}_{1}$ and $\mathcal{D}_{3}$:

[TABLE]

where $P_{dis}^{st}=\int_{x>a}p^{st}(x)dx$. So, it implies that the ergodic hypothesis is valid for the occupation fraction in any domain $\mathcal{D}$.

III.2 Simulations

We first simulate the steady solution $p^{st}(x)$ in (17) for the case of $\beta>\alpha$. In Fig. 1, they are specified as $\beta=2$ and $\alpha=1$. In a regularized form, the space-dependent diffusivity and the potential well are $D(x)=|x|+A_{0}$ and $U(x)=(|x|+A_{0})^{2}/2$, respectively, where $A_{0}=0.01$. Then the theoretical solution is $p^{st}(x)\propto(|x|+A_{0})^{-1/2}\exp(-|x|)$, which is confirmed by the simulations in Fig. 1 with trajectories up to $10^{4}$ and the measurement time up to $10^{3}$. In Fig. 2, we present the PDFs of the occupation fraction in domain $\mathcal{D}=\{x<a\}$ with three different values $a=-0.5$, [math], and $0.5$, respectively. In simulations, the same $D(x)$ and $U(x)$ are used as those in Fig. 1. We use trajectories up to $10^{3}$ and the measurement time up to $10^{4}$. Similarly, we also simulate the PDFs of the occupation fraction for the case of $\alpha<0$ in Fig. 3. In this case, we take $D(x)=1/(1+|x|^{2})$ and $U(x)=|x|$ (i.e., $\beta=1$ and $\alpha=-2$). The simulations in Fig. 3 are also consistent to the theoretical result (19). The space-dependent diffusivity $D(x)$ with $\alpha=-2$ is common in biological system, e.g., the small fluorescently labelled proteins in the cytoplasm of mammalian NLFK and HeLa cells Kühn et al. (2011).

It deserves to be mentioned that the exponent $\alpha$ can be bigger than $2$ when reflecting boundary condition is applied. Since, in some sense, the reflecting boundary condition at $\pm L$ can be regarded as the potential well $U(x)=|x/L|^{\beta}/\beta$ with $\beta\rightarrow\infty$; for any given $\alpha$, the condition $\beta>\alpha$ is always satisfied and the result of the PDF of the occupation fraction is also valid. In this case, the steady solution is $p^{st}(x)\propto|x|^{-\alpha/2}$ in (17), which only depends on the diffusivity $D(x)$. The steady solution in bounded domain varies for different processes. Some interesting results have been found for Lévy flight and Lévy walk Dybiec et al. (2017), fractional Brownian motion Wada and Vojta (2018) in a bounded domain.

One specific example for the process driven by a Langevin equation with multiplicative noise in a bounded domain has been discussed in Wang, Chen, and Deng (2018b), where the PDF of the occupation fraction in positive half-space was derived to be a $\delta$-function based on the backward Feynman-Kac equation. The multiplicative noise in Wang, Chen, and Deng (2018b) was interpreted in the Itô sense, different from the Stratonovich sense considered in this paper. The results in this paper can be easily extended to other different interpretations of the multiplicative noises, including Itô sense Ito (1950) and Hänggi-Klimontovich sense Hänggi (1982); Klimontovich (1990). For different interpretations, one just needs to modify the external force term in the Fokker-Planck equation (13) and obtain a similar steady solution $p^{st}(x)$, which will not bring any difficulties for the discussions in this paper.

IV Infinite ergodic theory for $\beta=\alpha$

Now, we turn to the case of $\beta=\alpha>0$, where the steady solution obtained in (17) is

[TABLE]

It can be seen that this PDF decays as power-law and its tail now relies heavily on the magnitudes $D_{0}$ in the diffusivity $D(x)$ and $U_{0}$ in the potential well $U(x)$. The PDF (42) can be normalized only when $\alpha/2+\alpha U_{0}/D_{0}>1$; under this condition, the ensemble average of the observable occupation fraction in the positive half-space is finite, namely, $\int_{0}^{\infty}p^{st}(x)dx<\infty$ (note that $p^{st}(x)$ is regular at origin as mentioned in Sec. II). But if the observable of interest is second moment, it is only finite under a stronger request: $\alpha/2+\alpha U_{0}/D_{0}>3$, which shows that a relatively stronger external potential well or a larger radio $U_{0}/D_{0}$ is responsible for a finite higher-order moment with respect to the steady solution. What happens when these requests are not satisfied? As we all know, the higher-order moments mainly depend on the tail of the PDF. But the steady solution $p^{st}(x)$ do not provide a quite correct information at its tail when we measure the higher-order moments. In turn, $p^{st}(x)$ is correct enough to calculate some suitable statistics. For example, we can get the correct ensemble averaged occupation fraction rather than second moment under the condition $1<\alpha/2+\alpha U_{0}/D_{0}<3$.

What about the ensemble averaged occupation fraction under the condition $\alpha/2+\alpha U_{0}/D_{0}<1$? Since $p^{st}(x)$ fails to describe the information at its tail, instead, we consider a time-dependent solution to obtain a correct decay rate. The common methods to obtain the solution are the eigenfunction expansion Risken (1989); Metzler, Barkai, and Klafter (1999) and scaling ansatz Kessler and Barkai (2010); Barkai, Aghion, and Kessler (2014) of the Fokker-Planck equation (13). We will apply the latter to solve (13) with the general $U(x)$ and $D(x)$. Due to the complexity of the system for any $\alpha<2$ (when $\alpha=2$, $\alpha/2+\alpha U_{0}/D_{0}<1$ does not hold), we will simplify the problem firstly through the variable substitution (6). Then, Eq. (12) reduces to the Langevin equation of $y(t)$:

[TABLE]

where

[TABLE]

This system describes the motion of a Brownian particle in a logarithmic potential, which has been discussed in detail in Lutz (2004); Dechant et al. (2011a). The solution $y(t)$ in (43) is also named as Bessel process, being related to a free Brownian particle starting from the origin in high dimensions. The exact form of its PDF $p(y,t)$ can be found in Bray (2000); Martin, Behn, and Germano (2011); Dechant et al. (2011b); Leibovich and Barkai (2018). On the other hand, to avoid complicated calculations, it is easy to do scaling ansatz for Eq. (43), which has been mentioned in Kessler and Barkai (2010), and we present the details in Appendix B for completeness. For $|y|\ll t^{1/2}$, we obtain the asymptotic behavior of the solution $p(y,t)$ in (75)

[TABLE]

with the normalization parameter $N_{g}=2^{A-1/2}[\Gamma(1/2-A)]^{-1}$. Considering the relation between $y(t)$ and $x(t)$ in (6), i.e., $|y|=\frac{\sqrt{2/D_{0}}}{2-\alpha}|x|^{\frac{2-\alpha}{2}}$, we obtain the asymptotic behavior of the solution $p(x,t)$: for $|x|\ll t^{1/(2-\alpha)}$,

[TABLE]

where $N_{x}=(2-\alpha)^{2A}D_{0}^{A-1/2}[2\Gamma(1/2-A)]^{-1}$. Note that the exponent of $x$ is consistent with the steady solution $p^{st}(x)$ in (42), which shows the potential value of the non-normalized steady solution. This consistency also explains why we consider the asymptotic behavior in $|x|\ll t^{1/(2-\alpha)}$ (i.e., the small $z$ in scaling ansatz in Appendix B). Moreover, there is a relation between them

[TABLE]

From (47), the parameter $N$ in (42) can be determined as $N=(2-\alpha)^{2A}D_{0}^{A}[2\Gamma(1/2-A)]^{-1}$. It looks natural that $p^{st}(x)$ is non-normalized, since the term $t^{\frac{1}{2}-A}$ tends to infinity for $A<1/2$. On the other hand, though $p^{st}(x)$ is non-normalized, it is still useful since it only differs from the real $p(x,t)$ by a known term of $t$. With this relation, the infinite-ergodic theory can be built up, uncovering the relation between the time average and ensemble average of an observable.

Considering an observable $\mathcal{O}(x)$, its ensemble average is defined as

[TABLE]

which can be represented through the non-normalized steady solution in (47) for long times

[TABLE]

if the integral on the right-hand side exists. This condition is not rigid for many observables, especially for those decaying for large $x$ which cure the non-integrability of $p^{st}(x)$ at infinity. It also confirms the explanation that the steady solution $p^{st}(x)$ is valid except its tail-part. A common one is the pulse function $\Theta_{\mathcal{D}}(x)$ with respect to a finite interval $\mathcal{D}$ (e.g., $\mathcal{D}=[-a,a]$ in the following), which is defined as $\Theta_{\mathcal{D}}(x)=1$ for $x\in\mathcal{D}$ and $\Theta_{\mathcal{D}}(x)=0$ otherwise. So it vanishes at infinity and it is integrable with respect to $p^{st}(x)$.

Now, we turn to the time average of the observable $\mathcal{O}(x)$, which is defined as

[TABLE]

and the corresponding ensemble-time average is

[TABLE]

Taking the limit $t\rightarrow\infty$ and using Eq. (49), we obtain

[TABLE]

where the constant factor $2/(2A+1)$ comes from the time-integral of the term $t^{A-1/2}$ in (49). Compared to the result (19) in the case of $\beta>\alpha$, where the normalized steady solution $p^{st}(x)$ directly decides the equivalent time and ensemble average, Eq. (52) shows that time and ensemble averages are still related but only differ by a constant factor for long times. In addition, the two kinds of averages are both found to be dependent on the non-normalized steady state $p^{st}(x)$ and decay as $t^{A-\frac{1}{2}}$, which can be seen from equations (49) and (52). It means that the number of the particles in the finite domain $\mathcal{D}$ decays as $t^{A-\frac{1}{2}}$ and the diffusion behavior goes on even for long times. In the case of $A<1/2$, i.e., $U_{0}/D_{0}<(2-\alpha)/(2\alpha)$, the potential well is too weak, compared with the space-dependent diffusivity, to suppress the diffusion; more and more particles go to infinity and thus they are not counted by the observable $\mathcal{O}(x)$.

IV.1 PDF of the time averaged observable

The distribution of the time averaged occupation fraction $f(\bar{p}_{\mathcal{D}})$ has been shown to be a $\delta$-function in (19), implying the ergodic phase for $\beta>\alpha$; on the contrary, the distribution for the case $\beta=\alpha$ becomes different, which will be shown to be broad and asymmetric in the following. The scatter of the amplitude of the time average for a set of individual trajectories is a useful indicator for the degree of non-ergodicity and the classification of different processes Metzler et al. (2014). We are interested in the distribution of the time averaged observable $\overline{\mathcal{O}(x)}$ around its mean for long times, so we define the random variable of dimensionless

[TABLE]

Especially, we consider the observable $\mathcal{O}(x)$ to be the pulse function $\Theta_{\mathcal{D}}(x)$ with $\mathcal{D}=[-a,a]$. Then the time averaged observable $\overline{\mathcal{O}(x)}$ becomes the occupation fraction $p_{\mathcal{D}}$ in domain $\mathcal{D}$.

When $\beta=\alpha$ and $A<1/2$ in (44), the diffusion of the particles cannot be controlled by the potential well. More and more particles leave the bounded domain $\mathcal{D}$ and never come back, resulting in the decay of the number of the particles in $\mathcal{D}$ as $t^{A-\frac{1}{2}}$. The trajectory of a particle can be described by two states: leaving the domain $\mathcal{D}$ and returning to the domain. We show a simple sketch map in Fig. 4, where the negative half-space is omitted due to symmetry. The two states $x<a$ and $x>a$ have completely different properties. The mean of the sojourn time in the bounded domain $\mathcal{D}$ (on the left side $x<a$) is finite, being similar to the case of $\beta>\alpha$. But for the sojourn time outside the domain $\mathcal{D}$ (on the right side $x>a$), its PDF is heavy-tailed $\phi(\tau)\sim\tau^{-1-\kappa},\,0<\kappa<1$ for large $\tau$ with infinite mean value due to a large fluctuation. The events that the particle outside (inside) $\mathcal{D}$ crosses the division point $x=a$ are recurrent. These events are not correlated and the PDF of the number of renewals $N_{o}(t)\,(N_{i}(t))$ before time $t$ can be solved by renewal theory Cox (1962). Importantly, the occupation time $T_{\mathcal{D}}(t)$ inside $\mathcal{D}$ is proportional to $N_{o}(t)$ due to the finite mean of the sojourn time inside $\mathcal{D}$. Therefore, the distribution of the random variable $\eta$ is related to the distribution of $N_{o}(t)$:

[TABLE]

Now only the exponent $\kappa$ in the heavy-tailed law $\phi(\tau)$ remains to be determined to obtain the PDF of $N_{o}(t)/\langle N_{o}(t)\rangle$. Actually, it only depends on the parameter $A$. Since the mean of the sojourn time inside $\mathcal{D}$ is finite, the average number of renewals is $\langle N_{i}(t)\rangle\propto t$. Similarly, considering the PDF of the sojourn time outside the domain is heavy-tailed, we have $\langle N_{o}(t)\rangle\propto t^{\kappa}$ Cox (1962); Godrèche and Luck (2001); Wang et al. (2018). Then the net number of the particles inside $\mathcal{D}$ decays as $t^{\kappa-1}$, which should be consistent to the result previously obtained $t^{A-1/2}$, and thus it leads to

[TABLE]

Then the PDF of $\eta$ is known as the Mittag-Leffler distribution Aaronson (1997); Feller (1971):

[TABLE]

where $L_{\kappa}[\cdot]$ is the one-sided Lévy density with order $\kappa$. Note that $f(\eta)$ in (56) is obtained for large $t$ and it finally becomes independent on $t$. Though the total time $t$ in the renewal process consists of two parts: the time particles spend inside and outside $\mathcal{D}$, the first part inside $\mathcal{D}$ can be statistically omitted compared to the second part.

We have used $A$ to characterize the radio between the strength of the potential well $U(x)$ and space-dependent diffusivity $D(x)$ in (44), and then $A$ determines the parameter $\kappa$ in (55) and the distribution of the time average $f(\eta)$. It is known that $f(\eta)\rightarrow\delta(\eta-1)$, returning to the ergodic phase when $\kappa\rightarrow 1$. On the contrary, when $\kappa$ takes its minimum value $1/2$, the distribution $f(\eta)$ becomes the broadest, implying the most significant non-ergodic behavior. Therefore, for the value of $A$ turning from [math] to $1/2$, $\kappa$ is from $1/2$ to $1$, and the system changes from non-ergodicity to the ergodic phrase.

IV.2 Simulations

We first simulate the solution $p(x,t)$ in Fig. 5 for the case of $\beta=\alpha=1$. The regularized factor is taken to be $A_{0}=1$ and other parameters are specified as $D_{0}=1,U_{0}=1/4$, and then $A=1/4$ in (44). The theoretical solution in (46) is $p(x,t)\simeq N_{x}t^{-1/4}|x|^{-3/4}$ for $|x|\ll t$. The auxiliary solid line in Fig. 5 has the same slope as theoretical result $N_{x}|x|^{-3/4}$, which is also the non-normalized steady solution $p^{st}(x)$ in (42). The prefactor of the auxiliary solid line is a little different from $N_{x}$ due to the regularized factor $A_{0}$. As the time goes by, the non-normalized solution $p^{st}(x)$ (solid line) is approached.

In the upper panel of Fig. 6, we demonstrate the result of the radio $\langle\overline{\mathcal{O}(x)}\rangle/\langle\mathcal{O}(x)\rangle$ in (52); due to its dependence on $A$, we choose parameters $\alpha=\beta=1$, $A_{0}=1$, $D_{0}=1$, and different $U_{0}$ ($U_{0}=1/6,\,1/4$) corresponding to different $A=1/6,\,1/4$ and different theoretical radio $\langle\overline{\mathcal{O}(x)}\rangle/\langle\mathcal{O}(x)\rangle=3/2,\,4/3$ for long times. The observed domain is $\mathcal{D}=[-0.5,\,0.5]$, i.e., $a=0.5$. The agreement between theoretical and simulated results can be seen in the upper panel with the measurement time up to $10^{5}$. In the lower panel, we verify the PDF $f(\eta)$ of the time averaged occupation time $\eta$ with the same parameters as above. So the theoretical PDF is the Mittag-Leffler distribution with order $\kappa=2/3$ and $3/4$ in (56), respectively, which is shown as the solid line and consistent to the simulation results in the lower panel. This theoretical lines can be obtained by numerical inverse Laplace transform, since the one-sided Lévy density is $L_{\kappa}(z)=\mathcal{L}_{s\rightarrow z}^{-1}[\exp(-s^{\kappa})]$.

V Infinite ergodic theory for $\beta<\alpha$

In the case of $\beta<\alpha$ and $\alpha>0$, the steady solution has been given in (17) as

[TABLE]

which can not be non-normalized for any parameters due to its large-$x$ behavior $|x|^{-\alpha/2}$. Similar to the case of $\beta=\alpha$, we need to solve the Fokker-Planck equation (13) to obtain a time-dependent solution. Since the potential well $U(x)$ in this case is too weak to change the shape of the solution $p(x,t)$ in (7) for a free particle, we employ a different scaling ansatz for the time-dependent solution, and obtain

[TABLE]

for long times, where $\mathcal{Z}(t)=1/\sqrt{4\pi D_{0}t}$. The details of the derivation of (58) and $N=\sqrt{D_{0}}$ in (57) are provided in Appendix C. Similar to (47), Eq. (58) shows the relation between the non-normalized steady solution and the time-dependent solution for long times. That the term $\mathcal{Z}(t)$ decays as $t^{-1/2}$ shows the particles tend to accumulate at infinity. For some observable $\mathcal{O}(x)$, the relation between ensemble-time average and ensemble average of the observable can be obtained

[TABLE]

and they both decay as $t^{-1/2}$. This result is consistent to the case of $A=0$ ($\kappa=1/2$) in the previous section $\beta=\alpha$, which means that the minimum value of $\kappa$ is $1/2$ and the non-ergodic behavior will not be stronger than the case of $\kappa=1/2$. This specific value $1/2$ is associated with the Gaussian white noise $\xi(t)$. More precisely, the diffusion behavior of a free Brownian particle is characterized by the scaling $x\sim t^{1/2}$, implying that the number of the particles in a bounded domain decays as $t^{-1/2}$. Even in a heterogeneous media, where the particle undergoes subdiffusion or superdiffusion, the asymptotic behavior of its PDF for large $t$ is also $t^{-1/2}$ in (7) and (10). Similarly, we define the random variable $\eta$ to characterize the distribution of the time averaged observable, and obtain that its PDF $f(\eta)$ also obeys Mittag-Leffler distribution with the order $\kappa=1/2$.

V.1 Simulations

We first simulate the solution $p(x,t)$ in Fig. 7 for the case of $\alpha=1.5$ and $\beta=0.5$ with other parameters specified as $D_{0}=1$ and $U_{0}=2$. The theoretical solution is the non-normalized steady solution in (57) $p^{st}(x)\propto x^{-3/4}$ for large $x$, shown as the solid line in Fig. 7. As the time goes by, the non-normalized solution $p^{st}(x)$ (solid line) is approached.

In the upper and middle panels of Fig. 8, we demonstrate the results of the radio $\langle\overline{\mathcal{O}(x)}\rangle/\langle\mathcal{O}(x)\rangle$ in (52) with two sets of different parameters: $\alpha=1.5,\,\beta=0.3$; and $\alpha=1,\,\beta=0.3$. Other parameters are $D_{0}=U_{0}=1$. The observed domain is $\mathcal{D}=[-0.5,0.5]$. The agreement between theoretical and simulated results can be seen for total measurement time up to $10^{4}$ in the upper panel and $10^{5}$ in the middle panel. It implies that for stronger diffusivity $D(x)$ (larger $\alpha$), the system displays a more enhanced diffusion behavior and reaches the state (52) using less time. In the lower panel, we verify the PDF $f(\eta)$ of the time averaged occupation time $\eta$ with the same parameters as above and a set of new parameter with negative $\beta$: $\alpha=1,\beta=-5$. For the three sets of different parameters, all the theoretical PDFs are the Mittag-Leffler distribution with order $\kappa=1/2$ in (56), which is shown as the solid line, being consistent to the simulation results in the lower panel.

VI Summary and conclusions

As a conclusion, we condense the main results for three different cases into a diagram presented in Fig. 9. The horizontal and vertical axes are the parameters $D_{0}$ and $U_{0}$ in diffusivity $D(x)$ and potential well $U(x)$, respectively. It has been fully discussed that the competition role played by $D_{0}$ and $U_{0}$ for $\alpha=\beta$ decides whether the system is ergodic or not in terms of the observable—occupation time in Sec. IV. It can be noted that the dividing line shown in Fig. 9 is

[TABLE]

from the condition $A<1/2$ in (44). In the case of $\alpha=\beta$, the system is ergodic when $U_{0}/D_{0}$ is larger than the critical value $(2-\alpha)/(2\alpha)$ and non-ergodic otherwise. The distribution of the time averaged observable is a $\delta$-function in ergodic phase, i.e., $\kappa=1$ in the Mittag-Leffler distribution. In non-ergodic phase, the time average is a random variable and its distribution is the Mittag-Leffler distribution with order $\kappa$ depending on $U_{0}$, $D_{0}$, and $\alpha$, i.e., $\kappa=A+1/2$ with $A$ defined in (44). Compared to $\beta=\alpha$, the results are relatively simple for another two cases $\beta>\alpha$ and $\beta<\alpha$ in Sec. III and Sec. V, respectively. For any positive $U_{0}$ and $D_{0}$, the former is ergodic with $\kappa=1$ and the latter is non-ergodic with $\kappa=1/2$ as shown in Fig. 9.

For the case of $\beta=\alpha$, there is an interesting phenomenon with the dividing line in Fig. 9. Since $\alpha>0$ has been assumed in this case, we consider two extreme values [math] and $2$ for $\alpha$. If $\beta=\alpha\rightarrow 2$, the dividing line with slope $(2-\alpha)/(2\alpha)$ approaches the horizontal line ($D_{0}$-axis), implying the ergodic phase for any value of $U_{0}\neq 0$. In contrast, for $\beta=\alpha\rightarrow 0$, the dividing line approaches the vertical line ($U_{0}$-axis), implying the non-ergodic phase for any value of $D_{0}\neq 0$. One example of the latter is a free Brownian particle exhibiting normal diffusion.

We take the potential well $U(x)$ and space-dependent diffusivity $D(x)$ to be the power-law form, due to its simple way to analyze the large-$x$ behavior. For general $U(x)$ and $D(x)$, e.g., the asymmetric ones and the bistable potential, the methods used and the results obtained in this paper are still valid. We find that a steady state can be reached when the potential well is stronger than diffusivity, i.e., $\beta>\alpha$. Therefore, a harmonic potential ($\beta=2$) is strong enough to control the enhanced diffusion for any $\alpha\leq 2$. In this case, the system is ergodic in terms of the observable occupation fraction. But for a more energetic particle, such as Lévy flight with divergent MSD, a deeper potential well is needed Dybiec, Sokolov, and Chechkin (2011). If the potential well is not deep enough, the steady solution can not be normalized; for some discussions on the observables of interest, the time-dependent solution should be considered. In addition, infinite-ergodic theorem helps us to find the relation between time averages and ensemble averages, as well as the distribution of time averages.

For the case of $\beta=\alpha$, the Langevin equation with multiplicative noise is turned into a new Langevin equation with additive noise (i.e., the Brownian particle in a logarithmic potential) through the variable substitution. This way looks similar to the treatment of a Langevin equation with multiplicative noise $D(x)$ in non-Stratonovich sense (e.g., Itô or Hänggi-Klimontovich sense) in a recent paper Leibovich and Barkai (2018). Compared to it, we regularize $D(x)$ at the origin and investigate the effects of the large-$x$ behavior of $D(x)$ and $U(x)$ on the ergodic dynamics. For different cases: $\beta=\alpha$ and $\beta<\alpha$, we use different kinds of scaling ansatz to obtain the scaling properties of the time-dependent solution $p(x,t)$ and even an analytic solution.

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 Relation between the space-dependent diffusivity and the distribution of waiting time

We first consider the coupled Langevin equation

[TABLE]

where $\xi(t)$ and $\eta(t)$ are independent; it can describe a specific CTRW model Fogedby (1994), and the subordinator $t(s)$ decides the distribution of the waiting times of CTRW model, of which the characteristic function is assumed to be

[TABLE]

By combining the two sub-equations in (61) Cairoli and Baule (2015), the subordinated process $y(t):=x(s(t))$ can be rewritten as a new Langevin equation (64). In fact,

[TABLE]

where we technically add the $\delta$-function in the second line, use a relation: $\Theta(s-s(t))=1-\Theta(t-t(s))$ between inverse subordinator $s(t)$ and subordinator $t(s)$ Baule and Friedrich (2005) in the third line, and make integration by parts in the last line. Therefore, $y(t)$ satisfies

[TABLE]

where $\bar{\xi}(t)$ can be regarded as a new noise of the subordinated process $y(t)$; and letting $\xi(t)$ be white noise leads to the correlation function of $\bar{\xi}(t)$ as

[TABLE]

the concrete form of which can be obtained from its Laplace transform ($t_{1}\rightarrow\lambda_{1},t_{2}\rightarrow\lambda_{2}$) by using the characteristic function of $t(s)$ in (62), i.e.,

[TABLE]

Taking the inverse Laplace transform in the above equation, when $\Phi(\lambda)=\lambda/(2D(x))$, we obtain the same form of the correlation function of multiplicative white noise

[TABLE]

This specific $\Phi(\lambda)$ in (62) signifies that the waiting time in the associated CTRW model obeys exponential distribution with rate $2D(x)$, i.e.,

[TABLE]

The space-dependent rate $2D(x)$ denotes the number of renewals per unit time, which means that for large $D(x)$ the mean of waiting time is small and thus the diffusion is enhanced and vice versa. Note that no external fields are considered here for simplicity, since we aim to derive the new noise $\bar{\xi}(t)$ in (67) with the space-dependent correlation. Actually, the Langevin equation (12) is equivalent to the first equation in (64) with an external field $f(y)$.

Appendix B Scaling ansatz for $\beta=\alpha$

We employ the scaling ansatz with the scaling function $g(z)$ for $p(y,t)$ of (43) as:

[TABLE]

The Fokker-Planck equation corresponding to the Langevin equation of $y(t)$ (43) is

[TABLE]

Substituting (69) into the Fokker-Planck equation (70), we get the equation of the scaling function $g(z)$ after replacing the variable $y$ with $z$:

[TABLE]

Comparing the terms with $t$, we find the scaling solution can be obtained by eliminating the terms with $t$ when $\gamma=1/2$. Then for small $z$ the left-hand side of Eq. (71) tends to zero compared to the right-hand side, and it becomes

[TABLE]

and the solution of (72) is

[TABLE]

The latter one is rejected since a PDF must decay at the infinity. To obtain the concrete coefficient in front of $|z|^{-2A}$, the complete form of $g(z)$ including the scaling for large $z$ should be derived firstly. For this, we assume $g(z)=N_{g}|z|^{-2A}h(z)$, and obtain

[TABLE]

the solution of which is $\exp(-z^{2}/2)$. Then the complete form of $g(z)$ is

[TABLE]

with the normalized parameter $N_{g}$ being

[TABLE]

Substituting (74) into (69) yields

[TABLE]

when $|y|\ll t^{1/2}$.

Appendix C Scaling ansatz for $\beta<\alpha$

Considering the exact solution in (7) for a free particle in a heterogeneous media, we use the scaling ansatz as

[TABLE]

to keep the same shape of the solution in (7) for long times. Note that we remain the terms depending on time $t$ but hide the term $|x|^{-\alpha/2}$ from (7) by putting it into the auxiliary function $p_{0}(x)$. Substituting the $p(x,t)$ (76) into the Fokker-Planck equation (13), we obtain

[TABLE]

For large $t$, the terms containing $\partial q/\partial t$ and $\partial q/\partial x$ can be omitted, compared with $q$. The leading terms remained form the equation

[TABLE]

the solution of which is exactly the steady solution in (57), i.e.,

[TABLE]

So we finally get the time-dependent solution

[TABLE]

Since for large $x$, Eq. (57) shows that

[TABLE]

the coefficient $N$ can be determined through the normalization of $\lim_{t\rightarrow\infty}p(x,t)$ in (80), and

[TABLE]