Viscoelastic subdiffusion: From anomalous to normal

We study viscoelastic subdiffusion in bistable and periodic potentials within the generalized Langevin equation approach. Our results justify the (ultra)slow fluctuating rate view of the corresponding bistable non-Markovian dynamics which displays bursting and anticorrelation of the residence times in two potential wells. The transition kinetics is asymptotically stretched exponential when the potential barrier $V_0$ several times exceeds thermal energy $k_{B}T$ $\left[V_0\sim \left(2–10\right)k_{B}T\right]$ and it cannot be described by the non-Markovian rate theory (NMRT). The well-known NMRT result approximates, however, ever better with the increasing barrier height, the most probable logarithm of the residence times. Moreover, the rate description is gradually restored when the barrier height exceeds a fuzzy borderline which depends on the power-law exponent of free subdiffusion $\alpha$. Such a potential-free subdiffusion is ergodic. Surprisingly, in periodic potentials it is not sensitive to the barrier height in the long time asymptotic limit. However, the transient to this asymptotic regime is extremally slow and it does profoundly depend on the barrier height. The time scale of such subdiffusion can exceed the mean residence time in a potential well or in a finite spatial domain by many orders of magnitude. All these features are in sharp contrast with an alternative subdiffusion mechanism involving jumps among traps with the divergent mean residence time in these traps.

Figure 1

(Color online) Frictional power-law memory kernel (in units of ${\eta }_{\alpha }$) and its two different approximations versus time $t$ (arbitrary units) for $\alpha =0.5$ and ${\nu }_0={10}^3$. Notice that the approximation with $b=2$ practically coincides with the exact kernel in this plot and the choice $b=10$ is also a very good one in spite of logarithmic oscillations which are barely seen. The inset magnifies a part of plot.

Figure 2

(Color online) (a) Stochastic simulations of Eqs. (5, 6, 7) with embedding dimension $D=18$ vs exact solution of GLE [Eq. (3)] for free subdiffusion. Coordinate is scaled in some arbitrary units $L$ and time is scaled in units of ${\tau }_{\beta }={\left(L^2{\eta }_{\alpha }/k_{B}T\right)}^{1/\alpha }$; $\alpha =0.5$ and the damping parameter $r_{\beta }={\tau }_{\beta }/{\tau }_{\text{bal}}=10$ with ${\tau }_{\text{bal}}=L/v_T$ and $v_T=\sqrt{k_{B}T/m}$. Particles are initially localized at $x=0$ with Maxwellian distributed velocities. The number of particles $n={10}^4$ is used in simulations (stochastic Heun algorithm, time step $\Delta t={10}^{-4}$) for the ensemble averaging. The time averaging for a single trajectory is done using Eq. (A5). The agreement proves ergodicity. The inset compares two numerically exact solutions of GLE, one with the strict power-law kernel and one with its 16-exponential approximation. (b) Relative numerical error of our simulations (versus the exact solution) is presented. The relative deviation of the solution with the approximate memory kernel from the solution with exact power-law kernel is also presented. Notice the occurrence of logarithmic oscillations within less than 1% error margin, which, however, does not play any essential role.

Figure 3

(Color online) A sample trajectory of bistable transitions for $\beta V_0=6$. Time is given in units of ${\tau }_r^{\left(b\right)}={\left[{\eta }_{\alpha }{x}_0^2/\left(4V_0\right)\right]}^{1/\alpha }$ and coordinate in units of $x_0$. The residence time distributions are extracted by using the two thresholds (red lines). Stochastic Heun algorithm with the time step $\Delta t={10}^{-3}$ and Mersenne Twister pseudorandom number generator, combined with the Box-Muller algorithm, are used. Markovian embedding with $N=16$ is used as described in the text. $\alpha =0.5$, $r=10$, $b=10$, and ${\nu }_0=1000$.

Figure 4

(Color online) Two-dimensional distribution of the logarithmic-transformed residence times in the potential wells for $\beta V_0=2$. The black symbol corresponds to the non-Markovian rate theory result (for the studied 16-exponential expansion of the power-law memory kernel). The orientation and the structure of the plateau of maximal probabilities manifest anticorrelation of the residence times in two potential wells. This feature is completely beyond the non-Markovian rate theory. Statistical software R [62] is used for data analysis and to produce the plot.

Figure 5

(Color online) Survival probability in one potential well versus time (in units of ${\tau }_r^{\left(b\right)}={\left[{\eta }_{\alpha }{x}_0^2/\left(4V_0\right)\right]}^{1/\alpha }$) for $\beta V_0=2$. The asymptote is stretched exponential with $\gamma \approx 0.65$. Fit with Eq. (12) is done using $d=0.81$ and $R=2.67\times {10}^{-2}$. These values were derived from the numerical data using the first two moments of the numerical distribution and Eqs. (13, 14).

Figure 6

(Color online) The low-frequency part of the power spectrum of fluctuations for $\beta V_0=2$. The $1/f^{\gamma }$ feature for smallest frequencies is defined chiefly by the stretched-exponential asymptotics of the survival probability (cf. Fig. 3). The correlogram method from the singular spectrum analysis-multitaper method (SSA-MTM) toolkit for spectral analysis [63] is used to produce the spectrum.

Figure 7

(Color online) Two-dimensional distribution of the logarithmic-transformed residence times in the potential wells for $\beta V_0=6$. The non-Markovian rate theory result (black symbol) yields the most probable value of $\text{ln}\left({\tau }_{1,2}\right)$. Kinetics is, however, clearly nonexponential even asymptotically.

Figure 8

(Color online) Survival probability in one potential well versus time (in units of ${\tau }_r^{\left(b\right)}$) for $\beta V_0=6$. Asymptotics is stretched exponential with $\gamma \approx 0.76$. Initial kinetics is stretched exponential with $\gamma \approx \alpha$. Fit with Eq. (12) is done using $d=0.634$ and $R=4.48\times {10}^{-4}$. These values were derived from the numerical data using the first two moments of the distribution and Eqs. (13, 14).

Figure 9

(Color online) Survival probability in one potential well versus time (in units of ${\tau }_r^{\left(b\right)}$) for $\beta V_0=10$. Using a fit to the survival probability, the whole distribution is reasonably well described by a stretched exponential with $\gamma \approx 0.90$. The maximum likelihood fit to the numerical data yields but a slightly different value $\gamma \approx 0.85$ (cf. in Fig. 10), indicating that the whole distribution is not properly stretched exponential. It is rather described by Eq. (12) with $d=0.335$ and $R=8.26\times {10}^{-6}$. These values were derived from the numerical data using the first two moments of the distribution and Eqs. (13, 14).

Figure 10

(Color online) Power exponent of a (single) stretched-exponential fit to the overall residence time distribution versus the barrier height for two different subdiffusion exponents $\alpha$. The maximum likelihood approach is used to derive $\gamma$ from data. Notice that this $\gamma$ does not correspond to the asymptotic $\gamma$ in the previous figures and text but rather presents an average of $\gamma \left(t\right)$ changing slowly in time. This subtle difference diminishes when $\gamma$ approaches one. For $\alpha =0.75$, $C_{0.75}=1.885$; other parameters in the kernel approximation are the same as for $\alpha =0.5$. For $V_0$ exceeding some borderline value $V_c\left(\alpha ,T\right)$, $\gamma$ tends to one and kinetics becomes gradually single exponential.

Figure 11

(Color online) Diffusion in washboard potentials. The distance is measured in units of period $L$ and the time in units of ${\tau }_{\beta }={\left(\beta L^2{\eta }_{\alpha }\right)}^{1/\alpha }$. $\alpha =0.5$ and the damping parameter $r_{\beta }={\tau }_{\beta }/{\tau }_{\text{bal}}=10$ with ${\tau }_{\text{bal}}=L/v_T$ and $v_T=\sqrt{k_{B}T/m}$. Particles are initially localized at $x=0$ with Maxwellian distributed velocities. The number of particles $n={10}^4$ is used in simulations (stochastic Heun algorithm, time step $\Delta t={10}^{-4}$) for the ensemble averaging. Initially, the diffusional broadening is always ballistic due to inertia effects. In potentials, this regime is followed by transient rattling oscillations due to a combination of the cage effect and the influence of the potential force. On the space scale $\delta x>L$, diffusion becomes ultraslow $⟨\delta x^2\left(t\right)⟩\propto t^{{\alpha }_{\text{eff}}\left(t\right)}$ with a slowly changing ${\alpha }_{\text{eff}}\left(t\right)$ which depends on the potential height. Asymptotically, ${\alpha }_{\text{eff}}\left(t\to \infty \right)\to 0.5$. This universal asymptotics is almost reached for $\beta V_0=2$.