Brownian motion with alternately fluctuating diffusivity: Stretched-exponential and power-law relaxation

TL;DR

This paper develops a theoretical framework for Brownian motion with fluctuating diffusivity, revealing different relaxation behaviors such as stretched-exponential and power-law decay depending on the timescale and initial conditions.

Contribution

It introduces a novel theory of alternating renewal processes to analyze relaxation functions in Brownian motion with fluctuating diffusivity, covering exponential and power-law sojourn-time distributions.

Findings

Short-time relaxation shows exponential or stretched-exponential behavior.

Long-time relaxation exhibits power-law decay with an exponential cutoff.

Results are consistent with numerical simulations.

Abstract

We investigate Brownian motion with diffusivity alternately fluctuating between fast and slow states. We assume that sojourn-time distributions of these two states are given by exponential or power-law distributions. We develop a theory of alternating renewal processes to study a relaxation function which is expressed with an integral of the diffusivity over time. This relaxation function can be related to a position correlation function if the particle is in a harmonic potential, and to the self-intermediate scattering function if the potential force is absent. It is theoretically shown that, at short times, the exponential relaxation or the stretched-exponential relaxation are observed depending on the power law index of the sojourn-time distributions. In contrast, at long times, a power law decay with an exponential cutoff is observed. The dependencies on the initial ensembles (i.e., equilibrium or non-equilibrium initial ensembles) are also elucidated. These theoretical results are consistent with numerical simulations.

Equations

$$\frac{d\mathbf{r}(t)}{dt}=-\frac{\tilde{u}\mathbf{r}(t)}{\gamma (t)}+\sqrt{2D(t)}\mathbf{\xi }(t),$$

$${⟨\mathbf{\xi }(t)⟩}_{\mathbf{\xi }}=0,{⟨\mathbf{\xi }(t)\mathbf{\xi }(t^{\prime })⟩}_{\mathbf{\xi }}=\mathbf{I}\delta (t-t^{\prime }),$$

$${⟨\sqrt{D(t)}\mathbf{\xi }(t)\sqrt{D(t^{\prime })}\mathbf{\xi }(t^{\prime })⟩}_{\mathbf{\xi }}=D(t)\mathbf{I}\delta (t-t^{\prime }).$$

$$\frac{k_{B}T}{\gamma (t)}=D(t).$$

$$\Phi (t):=\frac{{⟨\mathbf{r}(t)\cdot \mathbf{r}(0)⟩}_{\mathbf{\xi },\mathrm{I},D}}{{⟨{\mathbf{r}}^2(0)⟩}_{\mathrm{I}}}={⟨\exp \left[-u{\int }_0^{t}d{t}^{\prime }D(t^{\prime })\right]⟩}_D,$$

$${⟨|\mathbf{r}(t)-\mathbf{r}(0){|}^2⟩}_{\mathbf{\xi },\mathrm{I},D}=2{⟨{\mathbf{r}}^2(0)⟩}_{\mathrm{I}}\left[1-\Phi (t)\right].$$

$$\frac{d\mathbf{r}(t)}{dt}=\sqrt{2D(t)}\mathbf{\xi }(t),$$

$$G_s(\mathbf{r},t):={⟨\delta \left(\mathbf{r}-\left[\mathbf{r}(t)-\mathbf{r}(0)\right]\right)⟩}_{\mathbf{\xi },\mathrm{I},D},$$

$$G_s[\mathbf{r},t|D]:={⟨\delta \left(\mathbf{r}-\left[\mathbf{r}(t)-\mathbf{r}(0)\right]\right)⟩}_{\mathbf{\xi },\mathrm{I}},$$

$$G_s(\mathbf{r},t)={⟨G_s\left[\mathbf{r},t|D\right]⟩}_D=\int \mathcal{D}D{G}_s\left[\mathbf{r},t|D\right]\mathcal{P}\left[\mathcal{D}\right],$$

$$\frac{\partial G_s[\mathbf{r},t|D]}{\partial t}=D(t){\nabla }^2{G}_s[\mathbf{r},t|D],$$

$$F_s[\mathbf{k},t|D]:=\int d\mathbf{r}{e}^{-i\mathbf{k}\cdot \mathbf{r}}{G}_s[\mathbf{r},t|D].$$

$$\frac{\partial F_s[\mathbf{k},t|D]}{\partial t}=-D(t)k^2{F}_s[\mathbf{k},t|D],$$

$$F_s[\mathbf{k},t|D]=\exp \left[-k^2{\int }_0^{t}D(t^{\prime })d{t}^{\prime }\right].$$

$$F_s(\mathbf{k},t)={⟨\exp \left[-k^2{\int }_0^{t}D(t^{\prime })d{t}^{\prime }\right]⟩}_D.$$

$$D(t):=\{\begin{array}{l}{D}_{+}(+\text{state; fast state}),\\ D_{-}(-\text{state; slow sate}).\end{array}$$

$$p_{\pm }^{\mathrm{eq}}=\frac{{\mu }_{\pm }}{{\mu }_{+}+{\mu }_{-}},$$

$$w_{+}^{\mathrm{eq}}({\tau }_1)+w_{-}^{\mathrm{eq}}({\tau }_1)=p_{+}^{\mathrm{eq}}{\rho }^{+,\mathrm{eq}}({\tau }_1)+p_{-}^{\mathrm{eq}}{\rho }^{-,\mathrm{eq}}({\tau }_1),$$

$$w_{+}^{\mathrm{neq}}({\tau }_1)+w_{-}^{\mathrm{neq}}({\tau }_1)=p_{+}^0{\rho }^{+}({\tau }_1)+p_{-}^0{\rho }^{-}({\tau }_1),$$

$$w_{+}^0({\tau }_1)+w_{-}^0({\tau }_1)=p_{+}^0{\rho }^{+,0}({\tau }_1)+p_{-}^0{\rho }^{-,0}({\tau }_1),$$

$${\hat{f}}^D(u;t|w^0):=⟨\exp \left[-u{\int }_0^{t}d{t}^{\prime }D(t^{\prime })\right]⟩.$$

$$\bar{D}(t):={\int }_0^{t}D(t^{\prime })d{t}^{\prime }$$

$$=\left(D_{+}-D_{-}\right)T_{+}(t)+D_{-}t,$$

$$f^{T_{-}}(t_{-};t|w^0)=f^{T_{+}}(t-t_{-};t|w^0).$$

$$f^D(\bar{d};t|w^0)=\frac{1}{D_{+}-D_{-}}{f}^{T_{+}}\left(\frac{\bar{d}-D_{-}t}{D_{+}-D_{-}};t|w^0\right)$$

$${\breve{f}}^D(u;s|w^0)={\breve{f}}^{T_{+}}\left((D_{+}-D_{-})u;s+D_{-}u|w^0\right),$$

$${\breve{f}}^{T_{+}}(u;s|w^0)=\sum _{n=0}^{\infty }$$

$$+$$

$$f_{\pm ,n}^{T_{+}}(t_{+};t|w^0)=⟨{\delta }_{\pm }\delta \left(t_{+}-T_{+}(t)\right)I(t_n\le t<t_{n+1})⟩,$$

$${\delta }_{\pm }=\{\begin{array}{ll}1,& \text{if the state is ± at t=0},\\ 0,& \text{if the state is ∓ at t=0}.\end{array}$$

$${\breve{f}}_{\pm ,n}^{T_{+}}(u;s|w^0)$$

Full text

I Introduction

The Brownian motion, which is random motion of microscopic particles suspended in fluids, is observed in various systems. The examples are colloid dynamics in suspensions Dhont (1996), beads and center-of-mass motion of polymer chains Doi and Edwards (1986), and macromolecular diffusion in intracellular environments Bressloff and Newby (2013). The Brownian motion is usually characterized by mean square displacements (MSDs). In fact, the MSD for a simple Brownian particle increases linearly with time, and the proportionality factor is the diffusion coefficient multiplied by $2d$ with $d$ being the spatial dimension of the system in which the Brownian particle is immersed. The diffusion coefficient of a rigid particle depends on its shape and the viscosity of the surrounding fluid Landau and Lifshitz (1987); Kim and Karrila (2013). For example, the diffusion coefficient of a spherical particle with a small Reynolds number is given by the Stokes law, $D=k_{B}T/6\pi\eta a$, with $k_{B}$ the Boltzmann constant, $T$ the absolute temperature, $\eta$ the viscosity of the fluids, and $a$ the radius of the particle.

The diffusion coefficient is usually considered as constant and is independent of time. However, this is not always true; in crowded systems such as supercooled liquids Yamamoto and Onuki (1998a, b) and intracellular environments Parry et al. (2014), the diffusion coefficient is not a constant, and often behaves like a random variable. For such systems, the diffusion coefficient should be interpreted as a time-dependent and fluctuating quantity. The center-of-mass motion of polymer chains actually exhibits time-dependent and fluctuating diffusivity Uneyama et al. (2012, 2015); Miyaguchi (2017). A prominent feature of such Brownian motion with fluctuating diffusivity is that its displacement distribution becomes non-Gaussian, though the MSD shows normal diffusion (i.e., the linear time dependence of the MSD mentioned above) Chubynsky and Slater (2014); Uneyama et al. (2015).

Recently, Brownian motion with fluctuating diffusivity has been intensively studied in the absence of the external force Chubynsky and Slater (2014); Uneyama et al. (2015); Manzo et al. (2015); Chechkin et al. (2017); Miyaguchi (2017); Jain and Sebastian (2017), but it would be also important to investigate the Brownian dynamics with some external forces, such as the forces by some confinements and external potentials. For example, the diffusion coefficient observed in single-particle-tracking experiments in bacterial cytoplasms is typically of the order of $10^{-2}\,[\text{\textmu}\mathrm{m}^{2}/\mathrm{sec}]$ Golding and Cox (2006); Parry et al. (2014), whereas the size of the bacterial cell is of the order of $1\text{\textmu}\mathrm{m}$. Thus, the time for diffusion across a cell is of the order of $10\,\mathrm{sec}$ Milo and Phillips (2015). On the other hand, the measurement times of single-particle-tracking experiments are often longer than $10\,\mathrm{sec}$ Golding and Cox (2006); Parry et al. (2014), and thus effects of confinement inside the cells cannot be ignored.

In Ref. Uneyama et al. (2019), the authors analyzed the Brownian particle with fluctuating diffusivity confined in a harmonic potential, that is, the Ornstein-Uhlenbeck process with fluctuating diffusivity (OUFD). In particular, a relaxation function $\Phi(t)$ of the position coordinates [See Eq. (5) below] was studied instead of the MSD [the MSD is easily derived from $\Phi(t)$, if the system is in equilibrium]. In Ref. Uneyama et al. (2019), the diffusion coefficient $D(t)$ is assumed to be Markovian stochastic processes, and an eigenmode expansion for the relaxation function $\Phi(t)$ is obtained through the eigenmode analysis of a transfer operator. Formally, these results can be applicable to any systems with the Markovian diffusivity $D(t)$. But it is also important to study systems with non-Markovian diffusivity $D(t)$, because some properties (such as a large scatter of a time-averaged MSD and weak ergodicity breaking) observed in single-particle-tracking experiments Golding and Cox (2006); Weber et al. (2010); Weigel et al. (2011); Jeon et al. (2011); Burov et al. (2011); Tabei et al. (2013); Yamamoto et al. (2014) are believed to originate from non-Markovian properties of the system He et al. (2008); Meroz et al. (2010); Jeon et al. (2011); Burov et al. (2011); Tabei et al. (2013); Miyaguchi and Akimoto (2011).

In this paper, we study the Brownian motion with $D(t)$ being a non-Markovian stochastic process. Because general analysis for the non-Markovian diffusivity $D(t)$ is extremely difficult, we assume that $D(t)$ takes only two values $D_{+}$ and $D_{-}$, and randomly switches between these two states (the fast and slow states). It would be worth noting that such two-state diffusion dynamics is actually a good approximation to describe protein diffusion along DNA strands Leith et al. (2012), large colloidal diffusion in bacterial cytoplasm Parry et al. (2014), and movement patterns of several kinds of bacteria Detcheverry (2017). To theoretically analyze such two-state dynamics, we use the alternating renewal theory Cox (1962); Goychuk and Hänggi (2003); Akimoto and Seki (2015); Miyaguchi et al. (2016), in which the system is assumed to switch randomly between the two states. Similar theoretical framework for the two-state dynamics is developed by utilizing non-commutable operator calculus in Ref. Detcheverry (2017).

Furthermore, sojourn-time distributions for the two states are assumed to be a power law in this paper, which is a typical non-Markovian stochastic process. As discussed in Ref. Miyaguchi et al. (2016), such power-law sojourn-time distributions can originate from the cage effects in crowded systems Odagaki and Hiwatari (1990); Doliwa and Heuer (2003); Helfferich et al. (2014) or spatial heterogeneity. Due to the power-law sojourn-time distribution, the relaxation function $\Phi(t)$ shows exponential or stretched-exponential decay at short times, whereas power-law decay with an exponential cutoff at long times.

This paper is organized as follows. In Sec. II, the OUFD and the Langevin equation with fluctuating diffusivity (LEFD) are defined. In Sec. II, the sojourn-time distribution and the initial ensembles are also introduced. In Sec. III, by utilizing an alternating renewal theory, the relaxation function is represented with the sojourn-time distributions. In Sec. IV, the relaxation function is derived for the exponential and the power-law sojourn-time distributions. Finally, Sec. V is devoted to a discussion. In Appendices, we summarize some technical matters, including simulation details.

II Definition of the models

In this section, we define the OUFD, and introduce a position correlation function $\Phi(t)$. In addition, the LEFD is also introduced, and it is shown that the self-intermediate scattering function is given by a formula similar to $\Phi(t)$.

II.1 Ornstein-Uhlenbeck process with fluctuating diffusivity

The Ornstein-Uhlenbeck process with fluctuating diffusivity is defined by the following equation Uneyama et al. (2015, 2019):

[TABLE]

where $\bm{r}(t)$ is the $d$-dimensional position vector of the Brownian particle at time $t$, $\gamma(t)$ is a friction coefficient, and $\tilde{u}$ is a constant which characterizes the strength of the restoring force [$\tilde{u}$ can be related to the spring constant of a harmonic potential $U(\bm{r})$, as $U(\bm{r})=\tilde{u}\bm{r}^{2}/2$].

The diffusion coefficient $D(t)$ is time dependent and fluctuating, i.e., $D(t)$ is a stochastic process. Moreover, $\bm{\xi}(t)$ is a white Gaussian noise which satisfies

[TABLE]

where $\bm{I}$ is the identity matrix, and we use the notation $\left\langle\dots\right\rangle_{\bm{\xi}}$ to represent an ensemble average over the noise history $\bm{\xi}(t)$. It is also assumed that $\bm{\xi}(t)$ and $D(t)$ are mutually independent.

Because of Eq. (2), the thermal noise is delta-correlated as

[TABLE]

We require that the system can reach the equilibrium state if the process $D(t)$ can reach the equilibrium state. This is equivalent to require that the detailed-balance condition (or the local equilibrium condition) is satisfied for the friction and noise coefficients Sekimoto (2010).

[TABLE]

Multiplying both sides of Eq. (1) by $\bm{r}(0)$, and averaging over a noise history $\bm{\xi}(t)$, we have an ordinary differential equation $d\left\langle\bm{r}(t)\cdot\bm{r}(0)\right\rangle_{\bm{\xi}}/dt=-uD(t)\left\langle\bm{r}(t)\cdot\bm{r}(0)\right\rangle_{\bm{\xi}}$, where we define $u$ as $u:=\tilde{u}/k_{B}T$ and used $\left\langle\bm{\xi}(t)\cdot\bm{r}(0)\right\rangle_{\bm{\xi}}=0$. Integrating this equation and averaging over realizations of $D(t)$ as well as initial distributions of $\bm{r}(0)$, we have a relaxation function

[TABLE]

where the subscripts $D$ and $\mathrm{I}$ represent the ensemble averages over $D(t)$ and $\bm{r}(0)$, respectively.

It follows from Eq. (5) that the relaxation function $\Phi(t)$ is independent of the initial distributions of $\bm{r}(0)$; namely, $\Phi(t)$ is characterized only by the stochastic process $D(t)$. In numerical simulations, we use the local equilibrium distribution $(u/2\pi)^{d/2}\exp\left[-u\bm{r}^{2}(0)/2\right]$ as an initial distribution of $\bm{r}(0)$, but the results are independent of the choice of the initial distribution. Therefore, whether $\bm{r}(0)$ follows the local equilibrium distribution or not is unimportant concerning the relaxation function $\Phi(t)$. From Eq. (5), on the other hand, the relaxation function clearly depends on $D(t)$ and thus the statistical properties of the stochastic process $D(t)$ is important. In this work, we consider rather general stochastic processes for $D(t)$ including the one that may not reach the equilibrium state (See Sec. IV.3). We use the words equilibrium and non-equilibrium to specify the initial ensembles of the stochastic process $D(t)$. We will explain these initial ensembles in Sec. II.3 and Appendix A.

Moreover, if the system is in a stationary state—i.e., if $\bm{r}(0)$ follows the local equilibrium distribution and $D(t)$ is a stationary stochastic process—the MSD is given by Uneyama et al. (2019)

[TABLE]

Thus, the MSD is easily obtained from the relaxation function $\Phi(t)$, if the system is stationary.

II.2 Langevin equation with fluctuating diffusivity

In this subsection, we show that the relaxation function $\Phi(t)$ given by Eq. (5) is related to the self-intermediate scattering function for the LEFD

[TABLE]

where the potential term is absent, and thus it is translationally invariant in space. To study dynamical behavior of such a system, a density correlation function $\left\langle\delta(\bm{r}^{\prime\prime}-\bm{r}(t^{\prime\prime}))\delta(\bm{r}^{\prime}-\bm{r}(t^{\prime}))\right\rangle_{\bm{\xi},D}$ is convenient. Due to the translational symmetry, this correlation function depends on the position coordinates $\bm{r}^{\prime\prime}$ and $\bm{r}^{\prime}$ only through the displacement $\bm{r}:=\bm{r}^{\prime\prime}-\bm{r}^{\prime}$. By also using a lag time $t:=t^{\prime\prime}-t^{\prime}$, the correlation function is rewritten as $\left\langle\delta(\bm{r}+\bm{r}^{\prime}-\bm{r}(t+t^{\prime}))\delta(\bm{r}^{\prime}-\bm{r}(t^{\prime}))\right\rangle_{\bm{\xi},D}$. Integrating with $\bm{r}^{\prime}$, we obtain Hansen and McDonald (1990)

[TABLE]

where we set $t^{\prime}=0$ and take the average over $\bm{r}(0)$. Note that $D(t)$ is not stationary in general, thus $G_{s}(\bm{r},t)$ depends on the choice of the origin of time. This distribution function for the particle displacement $G_{s}(\bm{r},t)$ is called the van Hove self-correlation function and it is widely utilized to study the dynamical behavior of glass formers [From the spherical symmetry, the angular integration of $G_{s}(\bm{r},t)$ gives $4\pi r^{2}G_{s}(r,t)$ with $r=|\bm{r}|$, and this latter function is usually utilized Kob and Andersen (1995)].

Here, let us consider the same quantity with a given $D(t)$.

[TABLE]

where the ensemble average for $D(t)$ is not taken.

Note that this correlation function is a functional of $D(t)$ and is related to $G_{s}(\bm{r},t)$ by

[TABLE]

where $\int\mathcal{D}\!D$ represents the path integral (functional integral) over the diffusivity $D(t)$ and $\mathcal{P}[\mathcal{D}]$ is a path probability of $D(t)$.

The functional $G_{s}[\bm{r},t|D]$ is equivalent to a propagator of the Brownian particle for a given $D(t)$, and thus it follows the Fokker-Planck equation Berne and Pecora (2000)

[TABLE]

with an initial condition $G_{s}[\bm{r},0|D]=\delta(\bm{r})$. Here, $\nabla^{2}$ is the Laplacian operator. This Fokker-Planck equation is solved by using a Fourier transform

[TABLE]

This functional is related to the self-intermediate scattering function as shown below. The Fourier transform of Eq. (11) is given by

[TABLE]

with an initial condition $F_{s}[\bm{k},0|D]=1$ and $k=|\bm{k}|$. This differential equation is easily solved, and we obtain

[TABLE]

The Fourier transform $F_{s}(\bm{k},t)$ of the correlation function $G_{s}(\bm{r},t)$ is usually referred to as the self-intermediate scattering function and widely utilized to study dynamical behavior of glass formers. By using Eq. (14), the function $F_{s}(\bm{k},t)$ is given by

[TABLE]

Note that this function $F_{s}(\bm{k},t)$ has the same form as $\Phi(t)$ in Eq. (5). A similar result was also derived in Ref. Jain and Sebastian (2017) in a slightly different way.

II.3 Definition of initial ensemble

In this subsection, the stochastic process $D(t)$ is defined as a two-state process, and sojourn-time distributions for these states are introduced. If the stochastic process $D(t)$ is a Markovian process, the system can reach the equilibrium state if the detailed blance is satisfied (for the Markovian processes, whether the detailed blance is satisfied or not can be checked rather straightforwardly). On the other hand, in this work, we consider more general processes for $D(t)$ which can be non-Markovian. Thus the situation is (at least apparently) not intuitive. Here we briefly explain the defintions of equilibrium and non-equilibrium initial states. The detailed derivations for them are presented in Appendix A.

In this paper, the diffusion coefficient $D(t)$ is assumed to be a dichotomous process as illustrated in Fig. 1. Namely, at each time $t$, the value of $D(t)$ is either of the two values $D_{+}$ or $D_{-}$:

[TABLE]

In addition, we assume that $0\leq D_{-}<D_{+}$. The condition $D_{-}\geq 0$ is used for calculating double unilateral Laplace transforms in Eq. (25) below. Contrastingly, in Ref. Godrèche and Luck (2001), a similar calculation is carried out by using a bilateral Laplace transform, thereby obtaining an equation corresponding to Eq. (25). This difference in computation originates from the fact that, in Ref. Godrèche and Luck (2001), a value corresponding to $D_{-}$ is set as $-1$ (In Ref. Godrèche and Luck (2001), this corresponding variable is not the diffusion coefficient, but it is related to the position of a diffusing particle. Thus, the negative value is physically relevant).

We define the dichotomous process $D(t)$ by using the times at which the transitions between the two states occur, $t_{n}\,(n=1,2,\dots)$; these times are referred to as renewal times [we also set $t_{0}=0$ for convenience, though $t_{0}$ is not a renewal time in general. See Fig. 1(a)]. The renewal time $t_{n}\,(n=1,2,\dots)$ can be written by a sum of successive sojourn times in the two states, $\tau_{k}$ ($k=1,2,\dots$), as $t_{n}\equiv\sum_{k=1}^{n}\tau_{k}$. A switching rule between the two states can then be given by sojourn-time probability density functions (PDFs) of these states. Namely, the sojourn times $\tau_{k}\,(k=2,3,\dots)$ for $+$ and $-$ states are random variables that follow sojourn-time PDFs, $\rho^{+}(\tau)$ and $\rho^{-}(\tau)$, respectively.

To fully specify the process $D(t)$, we also need to define the PDF for the first sojourn time $\tau_{1}$ (See Fig. 1). We mainly study two kinds of the first sojourn-time PDFs: the equilibrium and non-equilibrium ensembles.

II.3.1 Equilibrium ensemble

For the equilibrium ensemble, the process starts at $t=-\infty$, and thus the system is in equilibrium at $t=0$. The diffusivity $D(0)$ is therefore at $+$ or $-$ state with equilibrium fractions $p_{+}^{\mathrm{eq}}$ and $p_{-}^{\mathrm{eq}}$. These fractions $p_{\pm}^{\mathrm{eq}}$ are simply given by

[TABLE]

where $\mu_{\pm}$ are mean sojourn times for $\rho^{\pm}(\tau)$. From this expression, it is clear that the equilibrium ensemble exists only when the mean sojourn times $\mu_{\pm}$ are finite; if $\mu_{\pm}$ diverges, the system cannot reach the equilibrium state Miyaguchi and Akimoto (2013).

Moreover, the first sojourn time $\tau_{1}$ follows an equilibrium sojourn-time PDF $\rho^{\pm,\mathrm{eq}}(\tau_{1})$, not $\rho^{\pm}(\tau_{1})$. This is because, in the equilibrium ensemble, $t=0$ is not a renewal time as illustrated in Fig. 1(a). Thus, the PDF of $\tau_{1}$ in the equilibrium ensemble is given by

[TABLE]

where $w^{\mathrm{eq}}_{\pm}(\tau_{1})$ is defined by $w^{\mathrm{eq}}_{\pm}(\tau_{1}):=p_{\pm}^{\mathrm{eq}}\rho^{\pm,\mathrm{eq}}(\tau_{1})$. Namely, $\rho^{\pm,eq}(\tau_{1})$ is a conditional PDF of the first sojourn time $\tau_{1}$ given that the process starts at $t=-\infty$ and the state is $\pm$ at $t=0$. For a derivation and an explicit expression of $\rho^{\pm,\mathrm{eq}}(\tau)$, see Appendix A.

II.3.2 Non-equilibrium ensemble

We also study a non-equilibrium ensemble. Here, we limit ourselves that a transition occurs exactly at the inital time $t=0$ [See Fig. 1(b)], and call such an ensemble as non-equilibrium. Note that this is a typical non-equilibrium ensemble employed in the renewal theory Godrèche and Luck (2001) and CTRW theory Miyaguchi and Akimoto (2013). In this ensemble, the diffusivity $D(0)$ is at $+$ or $-$ state with some inital fractions $p_{+}^{0}$ and $p_{-}^{0}$; one of these fractions is arbitrary, and the other is given by the normalization $p_{+}^{0}+p_{-}^{0}=1$.

Moreover, if $D(0)=D_{\pm}$, the first sojourn time $\tau_{1}$ follows the sojourn-time PDF $\rho^{\pm}(\tau_{1})$; i.e., the first sojourn-time PDF is the same as those for the following sojourn times $\tau_{n}\,(n=3,5,\dots)$. This is because we assume that $t=0$ is a renewal time for the non-equilibrium ensemble. It follows that, the first sojourn-time PDF is given by

[TABLE]

where $w^{\mathrm{neq}}_{\pm}(\tau_{1})$ is defined by $w^{\mathrm{neq}}_{\pm}(\tau_{1}):=p_{\pm}^{0}\rho^{\pm}(\tau_{1})$.

II.3.3 General ensemble

In the next section, we present a theory for a more general inital ensemble rather than the specific ensembles defined above. Let us denote this general inital ensemble as

[TABLE]

where $w^{0}_{\pm}(\tau_{1}):=p_{\pm}^{0}\rho^{\pm,0}(\tau_{1})$, and $\rho^{\pm,0}(\tau)$ is the first sojourn-time PDF given that the initial state is $\pm$. Therefore, if we replace $p_{\pm}^{0}$ and $\rho^{\pm,0}(\tau_{1})$ in Eq. (20) with $p_{\pm}^{\mathrm{eq}}$ and $\rho^{\pm,\mathrm{eq}}(\tau_{1})$, we obtain the equilibrium ensemble defined in Eq. (18); similarly, if we replace $\rho^{\pm,0}(\tau_{1})$ in Eq. (20) with $\rho^{\pm}(\tau_{1})$, we obtain the non-equilibrium ensemble defined in Eq. (19). These replacement rules are used in the following sections to obtain results for the equilibrium and non-equilibrium ensembles.

We also use notations such as $w^{0}(\tau_{1}):=w_{+}^{0}(\tau_{1})+w_{-}^{0}(\tau_{1})$, $w^{\mathrm{eq}}(\tau_{1}):=w_{+}^{\mathrm{eq}}(\tau_{1})+w_{-}^{\mathrm{eq}}(\tau_{1})$ and $w^{\mathrm{neq}}(\tau_{1}):=w_{+}^{\mathrm{neq}}(\tau_{1})+w_{-}^{\mathrm{neq}}(\tau_{1})$, with which we can completely specify the initial ensemble for $D(t)$. Moreover, instead of the notation $\left\langle\dots\right\rangle_{D}$ used in the previous subsections, hereafter we employ the following notations for the ensemble average over $D(t)$. We denote the average in terms of the equilibrium ensembles $w^{\mathrm{eq}}(\tau_{1})$ as $\left\langle\dots\right\rangle_{\mathrm{eq}}$, and the average in terms of the genelal initial ensemble $w^{0}(\tau_{1})$ as $\left\langle\dots\right\rangle$ (the bracket without a subscript).

III Alternating renewal theory: general analysis

In this section, we derive a general formula for the relaxation function $\Phi(t)$ by using the alternating renewal theory Cox (1962); Miyaguchi et al. (2016), in which the system has two states with different sojourn-time PDFs in general. In fact, we show that $\Phi(t)$ is closely related to a quantity which is called occupation times in the renewal theory Godrèche and Luck (2001).

III.1 Definition of occupation time

To explicitly represent the variables and initial conditions, let us denote the relaxation function $\Phi(t)$ as $\hat{f}^{D}(u;t|w^{0})$; namely [see Eq. (5)]

[TABLE]

Also, we define an integral quantity $\bar{D}(t)$ as

[TABLE]

where the occupation times $T_{+}(t)$ and $T_{-}(t)$ are the total times spent by the system in the $+$ and $-$ states up to time $t$, respectively. Then, $\hat{f}^{D}(u;t|w^{0})=\left\langle\exp(-u\bar{D}(t))\right\rangle$ is the Laplace transform of the PDF $f^{D}(\bar{d};t|w^{0})$ with the initial ensemble $w^{0}(\tau_{1})$. Note that if $D_{+}=1$ and $D_{-}=-1$, $\bar{D}(t)$ is often referred to as a mean magnetization by analogy between the two-state trajectory $D(t)=\pm 1$ and a spin configuration Godrèche and Luck (2001).

Let us denote the PDF of $T_{\pm}(t)$ with an initial ensemble $w^{0}(\tau_{1})$ as $f^{T_{\pm}}(t_{\pm};t|w^{0})$. Here, we focus on the PDF $f^{T_{+}}(t_{+};t|w_{0})$, because $f^{T_{-}}(t_{-};t|w^{0})$ can be readily obtained from $f^{T_{+}}(t_{+};t|w_{0})$ by

[TABLE]

Similarly, from Eq. (22), the PDF $f^{D}(\bar{d};t|w^{0})$ is also given by $f^{T_{+}}(t_{+};t|w^{0})$ as

[TABLE]

Carrying out double Laplace transforms of Eq. (24) in terms of $\bar{d}\leftrightarrow u$ and then $t\leftrightarrow s$, we have

[TABLE]

where we used the assumption $D_{-}\geq 0$ in the Laplace transform with respect to $\bar{d}$.

III.2 PDF of occupation time

To derive an explicit form of $\breve{f}^{T_{+}}\left(u;s|w^{0}\right)$, we define $f^{T_{+}}_{\pm,n}(t_{+};t|w^{0})$ as a conditional joint PDF that (i) the state is $\pm$ at time $t=0$, (ii) the number of the renewals up to time $t$ is $n$, and (iii) the occupation time of the $+$ state up to time $t$ is $t_{+}$ under the condition that the initial ensemble is $w_{0}$. The function $\breve{f}^{T_{+}}(u;s|w^{0})$ can be expressed by using $\breve{f}^{T_{+}}_{\pm,n}(u;s|w^{0})$ as

[TABLE]

Now, we derive each term in the right-hand side of Eq. (26). The joint PDF $f^{T_{+}}_{\pm,n}(t_{+};t|w^{0})$ can be written as Godrèche and Luck (2001)

[TABLE]

where $I(\dots)$ and $\delta_{\pm}$ are quantities that take $1$ or [math]. $I(\dots)$ becomes $1$ if the inside of the bracket is satisfied, and otherwise becomes 0. Moreover, $\delta_{\pm}:=\delta_{D(0),D_{\pm}}$ is a random variable indicating the initial state, i.e.,

[TABLE]

By taking the Laplace transforms of Eq. (27) in terms of $t_{+}\leftrightarrow u$ and $t\leftrightarrow s$, we have

[TABLE]

where $\left\langle\dots\right\rangle_{\pm}$ is a conditional average under the condition that the initial state is $\pm$.

The ensemble average in Eq. (29) can be calculated as follows. If the initial state is the $+$ state [i.e., $D(0)=D_{+}$], the occupation time $T_{+}(t)$ can be expressed as

[TABLE]

when $2n$ transitions occur up to time $t$, or

[TABLE]

when $2n-1$ transitions occur up to time $t$. Thus, from Eq. (29), for $n=0$ we have

[TABLE]

where we used $p_{+}^{0}\left\langle e^{-(u+s)\tau_{1}}\right\rangle=p_{+}^{0}\hat{\rho}^{+,0}(u+s)=\hat{w}_{+}^{0}(u+s)$. For $n=1,2,\dots$, we can obtain $\breve{f}^{T_{+}}_{+,n}(u;s|w^{0})$ with a calculation similar to Eq. (69) in Appendix A. From Eq. (29), we have

[TABLE]

for $n=1,2,\dots$. Here, the function $\hat{\rho}(u,s)$ is defined as $\hat{\rho}(u,s):=\hat{\rho}^{+}(u+s)\hat{\rho}^{-}(s)$.

Similarly, if the initial state is the $-$ state [i.e., $D(0)=D_{-}$], we have

[TABLE]

when $2n$ transitions occur up to time $t$, or

[TABLE]

when $2n-1$ transitions occur up to time $t$. Thus, from Eq. (29), we have

[TABLE]

for $n=1,2,\dots$.

Substituting Eqs. (32)–(34) and Eqs. (37)–(39) into (26), we have

[TABLE]

Now we can express the PDF of the occupation time $T^{+}(t)$ by using the initial fractions $p_{\pm}^{0}$, and the first and subsequent sojourn-time PDFs, $\hat{w}^{0}_{\pm}(s)$ and $\hat{\rho}^{\pm}(s)$.

III.3 Explicit form of relaxation function

Finally, we derive explicit forms of the relaxation function. From Eqs. (25) and (40), we have

[TABLE]

where $z_{\pm}$ is defined by $z_{\pm}:=D_{\pm}u+s$. Here, let us define functions $\tilde{\rho}^{\pm}(s)$ as

[TABLE]

By using these auxiliary functions, Eq. (III.3) can be rewritten as

[TABLE]

where we used $\hat{w}_{\pm}^{0}(s)=p_{\pm}^{0}\hat{\rho}^{0,\pm}(s)$. If the two diffusion coefficients are the same $D_{+}=D_{-}=:D$, we have $z_{+}=z_{-}$; and thus the second and fourth of terms in Eq. (III.3) vanish. The first and third terms give $\breve{f}^{D}(u;s|w^{0})=1/(s+Du)$, and thus we obtain an exponential decay: $\hat{f}^{D}(u;t|w^{0})=e^{-Dut}$. This result can be also obtained more directly from Eq. (21).

The double Laplace transform $\breve{f}^{D}(u;s|w^{\mathrm{eq}})$ for the equilibrium ensemble $w^{\mathrm{eq}}$ can be obtained simply by replacing $\hat{\rho}^{\pm,0}$ and $p_{\pm}^{0}$ with $\hat{\rho}^{\pm,\mathrm{eq}}$ and $p_{\pm}^{\mathrm{eq}}$ in Eq. (III.3), respectively. Similarly, $\breve{f}^{D}(u;s|w^{\mathrm{neq}})$ for the non-equilibrium ensemble is obtained by replacing $\hat{\rho}^{\pm,0}$ with $\hat{\rho}^{\pm}$.

IV Case studies

In this section, we study the relaxation function $\hat{f}^{D}(u;t|w^{0})$ for the cases in which the sojourn-time PDFs are given by the exponential distributions or the power-law distributions. For the case where the mean sojourn times exist (i.e., $\mu_{\pm}<\infty$), we assume that the switching between the two states is fast compared with the diffusive time scale—i.e., $\mu_{\pm}\ll{1}/{D_{+}u}$. Moreover, we focus on the time regime $t\gtrsim{1}/{D_{+}u}$; in the Laplace domain, this corresponds to $s\lesssim D_{+}u$. These fast switching assumptions can be rewritten as

[TABLE]

Thus, we obtain $\mu_{\pm}\ll 1/z_{\pm}$ or equivalently $z_{\pm}\ll k_{\pm}$, where $k_{\pm}$ is the transition rate and given by $k_{\pm}=1/\mu_{\pm}$.

IV.1 Exponential distribution

If both sojourn-time PDFs, $\rho^{\pm}(\tau)$, are given by the exponential distributions, the first sojourn-time PDF is given by

[TABLE]

Therefore, in this case, the general initial ensemble $w^{0}(\tau_{1})$ is equivalent to the typical non-equilibrium ensemble $w^{\text{neq}}(\tau_{1})$. This equivalence is due to the fact that the exponential distributions are Markovian; a direct proof is possible as follows. Generally, $\rho^{\pm,0}(\tau_{1})$ should be given by a forward-recurrence-time PDF with respect to $\rho^{\pm}(\tau)$, but this forward-recurrence-time PDF can be shown to be equivalent to $\rho^{\pm}(\tau)$ if $\rho^{\pm}(\tau)$ is an exponential distribution due to its Markovian property [See Eq. (81) in Appendix A]. By using Eqs. (44), (78), and (79), $\breve{f}^{D}(u;s|w^{0})$ in Eq. (III.3) can be rewritten as

[TABLE]

This function is equivalent to that for the typical non-equilibrium ensemble $\breve{f}^{D}(u;s|w^{\text{neq}})$.

For the case of the equilibrium ensemble, we have $w^{\mathrm{eq}}(\tau_{1})=p^{\mathrm{eq}}_{+}\rho^{+}(\tau_{1})+p^{\mathrm{eq}}_{-}\rho^{-}(\tau_{1})$, and Eq. (45) can be further simplified by replacing $p_{\pm}^{0}$ with $p_{\pm}^{\mathrm{eq}}$ as

[TABLE]

This equation is consistent with the result obtained in Ref. Uneyama et al. (2019), in which a transfer operator method is utilized. In fact, the relaxation function $\breve{f}^{D}(u;s|w^{\mathrm{eq}})$ has two simple poles in terms of $s$, which correspond to two different relaxation modes Uneyama et al. (2019).

If we assume that the switching between the two-state is fast, i.e., $z_{\pm}\ll k_{\pm}$, we have

[TABLE]

where $\left\langle D\right\rangle_{\mathrm{eq}}$ is defined by $\left\langle D\right\rangle_{\mathrm{eq}}:=D_{+}p_{+}^{\mathrm{eq}}+D_{-}p_{-}^{\mathrm{eq}}$. Therefore, in this limit, we obtain a single mode relaxation, as if the particle diffuses with the diffusion coefficient $\left\langle D\right\rangle_{\mathrm{eq}}$.

The same result is obtained from the relaxation function with the non-equilibrium initial ensemble [Eq. (45)]. This is because, under the assumption of the fast switching $\left\langle D\right\rangle_{\mathrm{eq}}u\ll k_{\pm}=1/\mu_{\pm}$, the relaxation of the fractions [Eq. (82)] is much faster than that given in Eq. (IV.1).

IV.2 Power-law distribution: $1<\alpha_{-}<2$

Let us study the case where the sojourn-time PDFs are given by the power-law distributions $\rho^{\pm}(\tau)\sim 1/\tau^{\alpha_{\pm}}$ with $1<\alpha_{\pm}\leq 2$. In this case, the mean sojourn times exist, and thus the equilibrium ensemble can be defined. In this and the next subsections, we assume that $D_{-}\ll D_{+}$.

IV.2.1 Equilibrium ensemble

If the sojourn-time distributions $\rho^{\pm}(\tau)$ follow the power-law distributions with finite means $\mu_{\pm}$, their Laplace transforms ($\tau\leftrightarrow\sigma$) can be expressed as

[TABLE]

where $1<\alpha_{-}<\alpha_{+}\leq 2$, and $o(\sigma^{\alpha_{\pm}})$ is Landau’s notation. When $\alpha_{+}=2$, $\rho^{+}(\tau)$ is the exponential distribution. Note also that the fast switching assumption between the two states, $z_{\pm}\ll k_{\pm}=1/\mu_{\pm}$, is necessary to use the expansion in Eq. (48) with $\sigma=z_{\pm}$. In addition to this condition, we assume that $\mu_{\pm}z_{\pm}\gg a_{\pm}z_{\pm}^{\alpha_{\pm}}$. The auxiliary function $\tilde{\rho}^{\pm}(\sigma)$ [Eq. (42)] is then given by

[TABLE]

The function $\breve{f}^{D}(u;s|w^{\mathrm{eq}})$ can be given in a form similar to Eq. (III.3). In fact, $\breve{f}^{D}(u;s|w^{\mathrm{eq}})$ can be obtained just by replacing $\hat{\rho}^{\pm,0}$ and $p_{\pm}^{0}$ in Eq. (III.3) with $\hat{\rho}^{\pm,\mathrm{eq}}$ and $p_{\pm}^{\mathrm{eq}}$, respectively. Furthermore, by using $\hat{\rho}^{\pm,\mathrm{eq}}(z_{\pm})=\tilde{\rho}^{\pm}(z_{\pm})/(\mu_{\pm}z_{\pm})$ [see Eq. (76)] and Eq. (17), we rewrite Eq. (III.3) as

[TABLE]

This expression for the equilibrium ensemble is valid for any sojourn-time PDFs $\rho^{\pm}(\tau)$ with finite means $\mu_{\pm}$.

Now, substituting $\tilde{\rho}^{\pm}(z_{\pm})\simeq\mu_{\pm}z_{\pm}-a_{\pm}z_{\pm}^{\alpha_{\pm}}$ into Eq. (50), we have a formula for the power-law sojourn-time PDFs as

[TABLE]

where $\mu$ and $\Delta D$ are defined by $\mu:=\mu_{+}+\mu_{-}$, and $\Delta D:=D_{+}-D_{-}$, respectively.

For $s\sim D_{+}u$, the first term in Eq. (IV.2.1) is dominant over the second, and thus the latter can be neglected. Therefore, we have

[TABLE]

or, by the Laplace inversion, we obtain

[TABLE]

Thus, at the short time scale $t\sim 1/D_{+}u$, the exponential relaxation is observed.

However, for $s\ll D_{+}u$, we have $s+\left\langle D\right\rangle_{\mathrm{eq}}\sim\left\langle D\right\rangle_{\mathrm{eq}}$ and $z_{+}\sim D_{+}u$. The constant terms must be exponentially small after the Laplace inversion, and thus, we have from Eq. (IV.2.1)

[TABLE]

By the Laplace inversion, we have the following expression in the time domain $t\gg 1/D_{+}u$

[TABLE]

Thus, in the case $1<\alpha_{-}<2$, the relaxation function exhibits the exponential decay at the short time scale [Eq. (53)], whereas it exhibits the power-law decay with an exponential cutoff at the long time scale [Eq. (55)]. As shown in Fig. 2, these theoretical results are consistent with numerical simulations. It is possible to show that Eqs. (52)– (55) are valid also for $\alpha_{\pm}>2$ (it should be assumed that $\alpha_{-}<\alpha_{+}$ and that $\alpha_{-}$ is not an integer). Thus, the power-law dacay $1/t^{1-\alpha_{-}}$ with any values of $\alpha_{-}>1$ is possible in principle. However, for $\alpha_{-}>2$, the decay is relatively fast, and therefore it might be difficult to observe the predicted power-law decay behavior even in numerical simulations.

IV.2.2 Non-equilibrium ensemble

The relaxation function for the typical non-equilibrium ensemble $w^{\mathrm{neq}}(\tau_{1})=p_{+}^{0}\rho^{+}(\tau_{1})+p_{-}^{0}\rho^{-}(\tau_{1})$ with $p_{+}^{0}+p_{-}^{0}=1$ can be derived in the same way as the case of the equilibrium ensemble. By replacing $\hat{\rho}^{\pm,0}(z_{\pm})$ in Eq. (III.3) with $\hat{\rho}^{\pm}(z_{\pm})$, and substituting Eqs. (48) and (49) into the resulting equation, we obtain

[TABLE]

Note that this expression for the typical non-equilibrium ensemble $w^{\mathrm{neq}}$ is valid for any sojourn-time PDF $\rho^{\pm}(\tau)$ with finite means $\mu_{\pm}$.

Now, substituting $\tilde{\rho}^{\pm}(z_{\pm})\simeq\mu_{\pm}z_{\pm}-a_{\pm}z_{\pm}^{\alpha_{\pm}}$, we have a formula for the power-law sojourn-time PDFs as

[TABLE]

For $s\sim D_{+}u$, the first term in Eq. (57) is dominant over the second, and thus the latter can be neglected. Therefore, we have

[TABLE]

or, by the Laplace inversion, we obtain

[TABLE]

Thus, at the short time scale $t\sim 1/D_{+}u$, the exponential relaxation, which is exactly the same as that in the equilibrium ensemble [Eq. (53)], is observed.

However, for $s\ll D_{+}u$, the constant terms can be neglected because they are exponentially small after the Laplace inversion, and thus we have

[TABLE]

or by the Laplace inversion, we obtain a power-law decay again for the long time regime $t\gg 1/D_{+}u$ as

[TABLE]

As shown in Fig. 3, these theoretical results are consistent with numerical simulation. These results [Eqs. (58)– (61)] are valid also for $\alpha_{\pm}>2$ (it should be assumed that $\alpha_{-}<\alpha_{+}$ and that $\alpha_{-}$ is not an integer).

As in the case of the equilibrium ensemble, the relaxation for the non-equilibrium ensemble is given by the exponential decay at the short time scale [Eq. (59)], whereas power-law decay with an exponential cutoff at the long time scale [Eq. (61)]. Thus, the behavior of both ensembles are qualitatively similar except that the power-law decay is slower for the equilibrium ensemble than for the non-equilibrium ensemble. In fact, the power-law index for equilibrium ensemble is greater by one than that for the non-equilibrium ensemble [See Eqs. (55) and (61)]. This difference is due to the fact that the power-law index of the initial sojourn-time PDF for the equilibrium ensemble is smaller by one than that for the non-equilibrium ensemble [See Eq. (87)]. This difference of power-law exponents in the initial sojourn-time PDFs often provides a remarkable initial-ensemble dependence of physical observables such as the correlation function and the diffusion coefficient Akimoto and Aizawa (2007); Akimoto et al. (2018).

IV.3 Power-law distribution: $0<\alpha_{-}<1$

Finally, let us study the case in which $0<\alpha_{-}<1$. In this case, the equilibrium distribution $\rho^{-,\mathrm{eq}}(\tau)$ no longer exists, because the mean value of $\rho^{-}(\tau)$ diverges. Hence, we consider only the non-equilibrium ensemble $w^{\mathrm{neq}}$. On the other hand, $\alpha_{+}$ is assumed to be $1<\alpha_{+}\leq 2$ (For $0<\alpha_{+}<1$, we can have similar results with different scaling exponents. But the behavior for these two cases are qualitatively the same, and thus here we focus only on the case $1<\alpha_{+}\leq 2$). In addition to the assumptions for $\hat{\rho}^{+}(z_{+})$ used in the previous subsection, we further assume that $a_{-}z_{-}^{\alpha_{-}}\ll 1$; this condition is necessary for the expansion given in Eq. (85).

The double Laplace transform $\breve{f}^{D}(u;s|w^{\mathrm{neq}})$ is obtained by replacing $\hat{\rho}^{\pm,0}(z_{\pm})$ in Eq. (III.3) with $\hat{\rho}^{\pm}(z_{\pm})$. Substituting $\hat{\rho}^{\pm}(z_{\pm})\simeq 1$, $\tilde{\rho}^{+}(z_{+})\simeq\mu_{+}z_{+}$ and $\tilde{\rho}^{-}(z_{-})\simeq a_{-}z_{-}^{\alpha_{-}}$ into this equation, we have for $s\sim D_{+}u$

[TABLE]

Strictly speaking, higher order terms should have been taken into account in this calculation, but the final result is unchanged because the higher order terms are canceled out. The Laplace inversion of Eq. (62) gives

[TABLE]

Thus, the relaxation function is given in a stretched-exponential form.

However, for $s\ll D_{+}u$, we have the following formula from the equation obtained by replacing $\hat{\rho}^{\pm,0}(z_{\pm})$ in Eq. (III.3) with $\hat{\rho}^{\pm}(z_{\pm})$

[TABLE]

where we also assumed that $z_{-}^{\alpha_{-}}\ll D_{+}u$ (if this condition is not satisfied, the power-law regime vanishes). The Laplace inversion gives

[TABLE]

Thus, in the case $0<\alpha_{-}<1$, the relaxation function exhibits the stretched-exponential decay at the short time scale, whereas the power-law decay with an exponential cutoff at the long time scale. As shown in Fig. 4, these theoretical results are consistent with numerical simulation.

V Discussion

Recently, Brownian motion with fluctuating diffusivity has been investigated intensively. To characterize such processes, various kinds of correlation functions are used. One of the most frequently used functions might be an ergodicity breaking parameter, which is a fourth order correlation function in terms of the position vector Uneyama et al. (2012, 2015, 2019); Miyaguchi et al. (2016); Miyaguchi (2017); Cherstvy et al. (2013); Cherstvy and Metzler (2016); Metzler et al. (2014). For the case of the LEFD, it is shown that the fluctuating diffusivity cannot be detected with the second order correlation including the MSD Miyaguchi (2017).

In this work, we have studied the relaxation function $\hat{f}^{R}(u;t|w^{0})$. This quantity can be interpreted as a second order position correlation function [Eq. (5)] for the OUFD, and as the self-intermediate scattering function [Eq. (15)] for the LEFD. This correlation function characterizes the fluctuating diffusivity in the sense that it is a Laplace transformation of the integral diffusivity $\int_{0}^{t}dt^{\prime}D(t^{\prime})$.

Under the assumption that the system has two diffusive states, we have derived a general formula [Eq. (III.3)] for the relaxation function $\hat{f}^{R}(u;t|w^{0})$ in terms of the sojourn-time PDFs $\rho^{\pm}(\tau)$. In particular, if the sojourn-time PDF of the slow state is given by a power-law distribution with the index $\alpha_{-}$, the relaxation function at the short time scales shows the exponential relaxation for $1<\alpha_{-}<2$, or the stretched-exponential relaxation for $0<\alpha_{-}<1$. In contrast, at the long time scales, the relaxation function shows a power law decay with an exponential cutoff.

The fact that our model can reproduce the stretched-exponential relaxation function is quite interesting. The stretched-exponential type relaxation behavior has been widely observed in various systems, and generally interpreted as the signature of a heterogeneity. Various models have been proposed to interpret the microscopic origin of the stretched-exponential relaxation Palmer et al. (1984); Odagaki and Hiwatari (1990); Ngai and Rendell (1993); Phillips (1996); Metzler and Klafter (2000). In most of such models, power-law type distributions for some physical quantities such as the waiting time are assumed.

Our model also assumes the power-law distribution for the sojourn time, and thus, in this aspect, it is similar to the conventional models. In particular, it would be worth mentioning a relation of the present model with the quenched-trap model studied in Ref. Odagaki and Hiwatari (1990). Statistical properties of the quenched-trap model with the spatial dimension $d>2$ is asymptotically the same as those of the continuous-time random walks Machta (1985). On the other hand, the LEFD also reduces to the continuous-time random walk in the limit $\mu_{+}\to 0$ with $D_{+}\mu_{+}=\mathrm{const}$ and $D_{-}=0$ Uneyama et al. (2015). Thus, the generation mechanism of the stretched-exponential relaxation (as well as the anomalous diffusion Miyaguchi et al. (2016)) in the present model may be similar to that in the quenched-trap model Odagaki and Hiwatari (1990).

In spite of these similarities between our model and previously studied systems, however, our model might have some advantages in that its dynamics is described by simple equations [Eqs. (1) and (7)]. This is in contrast to some conventional models, where the dynamic equation and the time evolution rule are not expressed in clear and simple equations. We expect that the combination of our model and experimental relaxation data can be employed to interpret molecular-level dynamics of heterogeneous systems. In fact, a model similar to ours was utilized to reproduce results of numerical simulations for glassy systems and supercooled liquids Chaudhuri et al. (2007); Helfferich et al. (2018).

Moreover, the stretched-exponential relaxation is obtained in our model when the equilibrium distribution $\rho^{-,\mathrm{eq}}(\tau)$ does not exist. This means that, the system is not in equilibrium, if it shows the stretched exponential type relaxation. This seems to be consistent with a naive expectation that the dynamic heterogeneity is a precursor of a glass transition which is clearly in non-equilibrium.

As we mentioned, the stretched-exponential relaxation is observed ubiquitously in complex and highly concentrated systems , but its origin is still unclear Phillips (1996). The stretched-exponential relaxation is often expressed as $e^{-(t/\bar{\tau})^{\beta}}$ where $\bar{\tau}$ and $\beta$ are constants; $\beta$ is typically around $0.5$ (experimental values are about $0.4\lesssim\beta\lesssim 0.6$). In our model, this exponent is simply given as $\beta=\alpha_{-}$, and therefore it reflects the tail of the sojourn-time PDF of the slow state at the long time region. This result implies that the power-law sojourn-time PDF for the slow state plays a fundamental role in emergence of the stretched-exponential relaxation. The power-law sojourn-time PDF, which is assumed in our model, is also generally observed in crowded systems Doliwa and Heuer (2003); Helfferich et al. (2014). Thus, our model might bridge a gap between these commonly observed phenomena: the power-law sojourn times and stretched-exponential relaxation.

Although it seems that the fast state plays a minor role, its property affects the characteristic relaxation time; in our model, the characteristic relaxation time $\bar{\tau}$ is given as $\bar{\tau}=[\mu_{+}\Delta Du/a_{-}\Gamma(1+\alpha_{+})]^{1/\alpha_{-}}$. Thus, the characteristic time scale depends on various parameters such as the sojourn-time PDF of the fast state and the two diffusion coefficients.

Acknowledgments

T.M. was supported by Grant-in-Aid (KAKENHI) for Scientific Research C (Grand No. JP18K03417). T.U. was supported by Grant-in-Aid (KAKENHI) for Scientific Research C (Grant No. JP16K05513). T.A. was supported by Grant-in-Aid (KAKENHI) for Scientific Research B (Grand No. JP16KT0021), and Scientific Research C (Grant No. JP18K03468).

Appendix A Equilibrium ensemble

As explained in Sec. II, the initial ensemble is fully specified by the PDF of the first sojourn time $\tau_{1}$ [Recall that the initial distribution of the position $\bm{r}(0)$ is unimportant in investigating the relaxation function $\Phi(t)$]. We derive an explicit expression for the Laplace transform of the equilibrium ensemble $\hat{\rho}^{\pm,\mathrm{eq}}(\sigma):=\int_{0}^{\infty}d\tau e^{-\tau\sigma}\rho^{\pm,\mathrm{eq}}(\tau)$.

The equilibrium ensemble can be derived from a forward-recurrence-time PDF. A forward-recurrence time is a sojourn time until the next transition from some elapsed time $t^{\prime}>0$ Godrèche and Luck (2001). Here, we study the forward-recurrence-time PDF $w_{\pm}(\tau;t^{\prime}|w^{0})$, from which we obtain the equilibrium ensemble by taking a limit $t^{\prime}\to\infty$.

A.1 PDF for forward-recurrence time

The forward-recurrence-time PDF $w_{\pm}(\tau;t^{\prime}|w^{0})d\tau$ is a joint probability that (i) the state is $\pm$ at time $t^{\prime}$, and (ii) the sojourn time from $t^{\prime}$ until the next transition is in an interval $[\tau,\tau+d\tau]$, given that the process starts with $w^{0}(\tau_{1})$ at $t=0$. Here, $\tau$ is called the forward-recurrence time, and it is illustrated in Fig. 5.

To obtain the forward-recurrence-time PDF, we further define the following PDFs for $n=0,1,\dots$:

[TABLE]

where $w_{\pm,n}(\tau;t^{\prime}|w^{0})d\tau$ is a joint probability that (i) the state is $\pm$ at $t=0$, (ii) there are $n$ transitions until time $t^{\prime}$, and (iii) the sojourn time at time $t^{\prime}$ is in an interval $[\tau,\tau+d\tau]$, given that the process starts with $w^{0}(\tau_{1})$ at $t=0$.

The double Laplace transforms $\breve{w}_{\pm,n}(\sigma;s|w^{0}):=\int_{0}^{\infty}dt^{\prime}e^{-t^{\prime}s}\int_{0}^{\infty}d\tau e^{-\tau\sigma}w_{\pm,n}(\tau;t^{\prime}|w^{0})$ with respect to $\tau$ and then $t^{\prime}$ result in ($\tau\leftrightarrow\sigma$ and $t^{\prime}\leftrightarrow s$)

[TABLE]

where $\left\langle\dots\right\rangle_{\pm}$ is a conditional average for a given initial state $\pm$. The ensemble averages in the above equation can be expressed with the Laplace transforms of the initial ensemble $\hat{w}_{\pm}^{0}(s)$ and the sojourn-time PDFs $\hat{\rho}^{\pm}(s)$. In fact, we have

[TABLE]

for $n=1,2,\cdots$, where $\hat{\rho}(s):=\hat{\rho}^{+}(s)\hat{\rho}^{-}(s)$ and $\left\langle\dots\right\rangle_{\pm}$ is a conditional average under the condition that the initial state is $\pm$. Then, using Eqs. (68) and (69) in Eq. (67), we have

[TABLE]

for $n=1,2,\dots$.

Because the Laplace transform of the forward-recurrence-time PDF, $w_{\pm}(\tau;t^{\prime}|w^{0})$, is given by $\breve{w}_{\pm}(\sigma;s|w^{0})=\sum_{n=0}^{\infty}[\breve{w}_{\pm,2n}(\sigma;s|w^{0})+\breve{w}_{\mp,2n+1}(\sigma;s|w^{0})]$, we have

[TABLE]

Equation (A.1) is a general expression of the forward-recurrence-time PDF for the alternating renewal processes.

A.2 Equilibrium ensemble

If both $\mu_{+}$ and $\mu_{-}$ exist, it follows from Eq. (A.1) that

[TABLE]

Since the above limit $t^{\prime}\to\infty$ is equivalent to putting the start time of the process to $-\infty$, the limitting function should be equivalent to the equilibrium ensemble $\hat{w}^{\mathrm{eq}}_{\pm}(\sigma)$:

[TABLE]

where $\hat{\rho}^{\pm,\mathrm{eq}}(\sigma)$ is defined as

[TABLE]

Hence, the equilibrium PDF exists only if $\mu_{\pm}$ is finite as explained in the main text. Note also that the equilibrium ensemble [Eq. (74)] does not depend on the initial ensemble $w^{0}$, as expected.

A.2.1 Exponential distribution

If the sojourn-time PDFs are given by the exponential distributions

[TABLE]

we have its Laplace transforms ($\tau\leftrightarrow\sigma$) as

[TABLE]

In this case, the equilibrium PDF [Eq. (76)] is equivalent to $\rho^{\pm}(\tau)$, i.e., it is easy to check

[TABLE]

Thus, the equilibrium ensemble is given by $p_{+}^{\mathrm{eq}}\rho_{+}(\tau_{1})+p_{-}^{\mathrm{eq}}\rho_{-}(\tau_{1})$.

More generally, the forward-recurrence-time PDF for the exponential sojourn times is calculated as follows. Here, let us assume that the initial ensemble is the typical non-equilibrium ensemble, i.e., $w_{\pm}^{\mathrm{neq}}(\tau_{1})=p_{+}^{0}\rho_{+}(\tau_{1})+p_{-}^{0}\rho_{-}(\tau_{1})$, where $p_{\pm}^{0}$ are the initial fractions with the constraint $p_{+}^{0}+p_{-}^{0}=1$. Then, after some straightforward calculations, the forward-recurrence-time PDF [Eq. (A.1)] is given as

[TABLE]

The double Laplace inversions give

[TABLE]

where $p_{\pm}(t^{\prime})$ are the fractions at time $t^{\prime}$, and are given as

[TABLE]

From Eqs. (81) and (82), the forward-recurrence-time PDF at $t^{\prime}$ is also represented as

[TABLE]

Thus, the forward-recurrence-time PDF with respect to the exponential distributions are also given by the exponential distributions, and only the fractions change with time. This is an outcome of the fact that the exponential distribution is Markovian. In addition, the fractions converge to the equilibrium fractions $p_{\pm}^{\mathrm{eq}}$ exponentially.

A.2.2 Power law distribution

The power law distribution is given by

[TABLE]

where $a_{\pm}$ is a scale factor and $\Gamma(-\alpha_{\pm})$ is the gamma function. Asymptotic forms of the Laplace transforms $\hat{\rho}_{\pm}(\sigma):=\int_{0}^{\infty}d\tau\,e^{-\tau\sigma}\rho_{\pm}(\tau)$ at small $\sigma$ are given by

[TABLE]

where $\mu_{\pm}$ is the mean sojourn time of the state $\pm$.

Now, let us assume that the sojourn-time PDF is given by a power law distribution. For $0<\alpha_{\pm}<1$, the mean sojourn time $\mu_{\pm}$ diverges, and therefore the equilibrium ensemble [Eq. (75)] does not exist. In contrast, for $1<\alpha_{\pm}\leq 2$, the equilibrium ensemble exists and is given by

[TABLE]

where we used Eq. (86). The important point is that the power exponent of $\sigma$ in Eq. (87) is smaller by one than that in the sojourn-time PDF [Eq. (86)]. This means that the mean initial sojourn time of the equilibrium ensemble with $1<\alpha_{\pm}<2$ diverges.

Appendix B Simulation setup

In numerical simulations, we used a simple sojourn-time PDF defined as

[TABLE]

where $\tau_{0}$ is a cutoff time for short trap times; we assume the same cutoff time for both $\rho_{+}(\tau)$ and $\rho_{-}(\tau)$. Then, by comparing Eq. (88) with Eq. (84), it is found that $a_{\pm}$ in Eq. (48) is given by $a_{\pm}=\tau_{0}^{\alpha_{\pm}}|\Gamma(1-\alpha_{\pm})|$. For $\alpha_{\pm}>1$, the mean sojourn times $\mu_{\pm}$ are given as

[TABLE]

To simulate equilibrium processes, we have to generate initial ensembles which follow the first sojourn-time PDFs $\rho_{\pm}^{\mathrm{eq}}(t)$ [Eq. (76)]. This can be achieved with a method presented in Ref. Miyaguchi and Akimoto (2013). First, let $\Theta$ be a random variable following a PDF

[TABLE]

where Eq. (89) is used. Moreover, let $X$ be a random variable, which follows the uniform PDF on the interval $(0,1)$. Then, a random variable $X\Theta$ follows the PDF $\rho^{\pm,\mathrm{eq}}(\tau)$ (See Ref. Miyaguchi and Akimoto (2013) for a derivation).

The Langevin equation in Eq. (1), or its discretized form,

[TABLE]

can be transformed into a dimensionless form with

[TABLE]

The remaining system parameters are $\alpha_{\pm}$, $\tilde{u}$ and the ratio $D_{-}/D_{+}$. In simulations, we set $D_{+}=1.0$ and $\tau_{0}=1$. For the numerical integration of the Langevin equation, the Euler method is employed Kloeden and Platen (2011).