Nonadiabatic escape and stochastic resonance

TL;DR

This paper develops an asymptotic method to analyze particle escape from metastable states under weak periodic forcing, extending to a double well potential to model stochastic resonance in the nonadiabatic limit.

Contribution

It introduces a novel asymptotic approach for solving the Fokker-Planck equation with mixed boundary conditions, applicable to nonadiabatic stochastic resonance scenarios.

Findings

Derived an approximate escape rate for weak periodic forcing

Extended the method to double well potentials with white noise

Formulated a two-state stochastic model in the nonadiabatic limit

Abstract

We analyze the fluctuation-driven escape of particles from a metastable state under the influence of a weak periodic force. We develop an asymptotic method to solve the appropriate Fokker-Planck equation with mixed natural and absorbing boundary conditions. The approach uses two boundary layers flanking an interior region; most of the probability is concentrated within the boundary layer near the metastable point of the potential and particles transit the interior region before exiting the domain through the other boundary layer, which is near the unstable maximal point of the potential. The dominant processes in each region are given by approximate time-dependent solutions matched to construct the approximate composite solution, which gives the rate of escape with weak periodic forcing. Using reflection we extend the method to a double well potential influenced by white noise and weak periodic forcing, and thereby derive a two-state stochastic model--the simplest treatment of stochastic resonance theory--in the nonadiabatic limit

Equations

$$\frac{\partial P}{\partial \tilde{t}}=\frac{\partial \tilde{J}}{\partial \tilde{x}}\text{where}$$

$$\tilde{J}=\left[\frac{d\tilde{U}}{d\tilde{x}}-A\cos (\omega \tilde{t})\right]P+{\sigma }^2\frac{\partial P}{\partial \tilde{x}},$$

$$ϵ=\frac{\sigma }{\sqrt{\Delta U}},r=\frac{A\Delta x}{{\sigma }^2},\text{and}\Omega =\frac{\omega (\Delta x{)}^2}{\Delta U}.$$

$$\frac{\partial P}{\partial t}=\frac{\partial J}{\partial x}=\frac{\partial }{\partial x}\left[\frac{dU}{dx}-{ϵ}^{2}r\cos (\Omega t)\right]P+{ϵ}^2\frac{{\partial }^{2}P}{\partial x^2},$$

$$\frac{\partial P_{B1}}{\partial t}=\frac{\partial }{\partial \eta }(a\eta P_{B1})-ϵr\cos (\Omega t)\frac{\partial P_{B1}}{\partial \eta }+\frac{{\partial }^2{P}_{B1}}{\partial {\eta }^2}.$$

$$P_{B1}(\eta ,t)=n_1\sqrt{\frac{a}{2\pi {ϵ}^2}}\text{exp}\left[-\frac{a}{2{ϵ}^2}(x+1{)}^2\right],$$

$$\frac{\partial P_{B2}}{\partial t}=-\frac{\partial }{\partial \zeta }\left(b\zeta P_{B2}\right)+\frac{{\partial }^2{P}_{B2}}{\partial {\zeta }^2}+O(ϵ).$$

$$P_{B2}\approx D_0\left[1-\sqrt{\frac{b}{2\pi {ϵ}^2}}{\int }_{-\infty }^x\text{exp}\left(-\frac{b{z}^2}{2{ϵ}^2}\right)dz\right]\text{exp}\left(\frac{b{x}^2}{2{ϵ}^2}\right).$$

$$P\sim \text{exp}\left[-{ϵ}^{-2}U+S(x,t)\right],$$

$$\frac{\partial S}{\partial t}+\frac{dU}{dx}\frac{\partial S}{\partial x}=r\frac{dU}{dx}\text{cos}(\Omega t).$$

$$\frac{dx}{dt}=\frac{dU}{dx}\equiv U^{\prime }(x)\text{and}$$

$$\frac{dS}{dx}=r\text{cos}(\Omega t),$$

$$S=S_{min}+r{\int }_{-1}^x\text{cos}\{\Omega \left[t+T(z)-T(x)\right]\}dz,$$

$$T(x)=\int \frac{dx}{U^{\prime }(x)}.$$

$$S$$

$$S_{max}\sim S_{min}+r{\int }_{-1}^{\delta }\text{cos}\{\Omega \left[t+T(z)-T_{max}\right]\}dz,$$

$$S_{min}\sim \text{log}{\left[P_{B1}{e}^{U/{ϵ}^2}\right]}_{\eta \to \infty }\simeq \text{log}\left(n_1\sqrt{\frac{a}{2\pi {ϵ}^2}}\right)-\frac{1}{{ϵ}^2},$$

$$S_{max}\sim \text{log}{\left[P_{B2}{e}^{U/{ϵ}^2}\right]}_{\zeta \to -\infty }\simeq \text{log}{D}_0.$$

$$P\approx$$

$$\text{exp}\left[-\frac{(1+U)}{{ϵ}^2}+r{\int }_{-1}^{\min (x,-ϵ)}\text{cos}\{\Omega \left[t+T(z)-T(x)\right]\}dz\right].$$

$$P\approx$$

$$\text{exp}\left[-\frac{(1+U)}{{ϵ}^2}+r[\min (x,-ϵ)+1]\text{cos}(\Omega t)+r\Omega \sin (\Omega t){\int }_{-1}^{\min (x,-ϵ)}(z+1)T^{\prime }(z)dz\right].$$

$$\frac{d}{dt}{\int }_{-\infty }^{0}Pdx={\int }_{-\infty }^0\frac{\partial J}{\partial x}dx=J{|}_{x=0},$$

$${\int }_{-\infty }^{0}Pdx\simeq ϵ{\int }_{-\infty }^{\infty }{P}_{B1}d\eta =n_1.$$

$$J{|}_{x=0}\simeq {{ϵ}^2\frac{\partial P}{\partial x}|}_{x=0}=-{ϵ}^2{D}_0\sqrt{\frac{b}{2\pi {ϵ}^2}}.$$

$$\frac{d{n}_1}{dt}=-\frac{\sqrt{ab}}{2\pi }\text{exp}\left(-\frac{1}{{ϵ}^2}+r{\int }_{-1}^{\delta }\text{cos}\{\Omega [t+T(z)-T_{max}]\}dz\right)n_1,$$

$$R$$

$$R=\frac{\sqrt{ab}}{2\pi }\text{exp}\left(-\frac{1}{{ϵ}^2}+r\text{cos}\Omega t+r\Omega \text{sin}\Omega t{\int }_{-1}^{\delta }(z+1)T^{\prime }(z)dz\right),$$

$$P_{ad}=n_1\sqrt{\frac{a}{2\pi {ϵ}^2}}\left[1-\sqrt{\frac{b}{2\pi {ϵ}^2}}{\int }_{-\infty }^x{e}^{-b{z}^2/2{ϵ}^2}dz\right]\text{exp}\left[-\frac{(1+U)}{{ϵ}^2}+r(x+1)\text{cos}(\Omega t)\right],$$

$$R_{ad}$$

$$\tilde{U}=U_{max}-\Delta U\left[\frac{(\tilde{x}-x_{max}{)}^2}{2(\Delta x{)}^2}+\frac{(\tilde{x}-x_{max}{)}^3}{3(\Delta x{)}^3}\right],$$

$$U=-3x^2-2x^3,\frac{dU}{dx}=-6x(x+1),\mathrm{and}a=b=6.$$

Full text

I Introduction

The escape of particles from a metastable state under the influence of noise is a classical problem in non-equilibrium statistical mechanics Kramers (1940). Calculating the rate of escape can be approached in a variety of ways (see e.g., Hänggi et al., 1990; van Kampen, 2007; Bressloff, 2017; Forgoston and Moore, 2018, and Refs. therein). An important extension of the original escape problem includes the influence of periodic forcing, with a phase that impacts the escape rate Nicolis (1982); Jung (1993). Of particular relevance here is the problem of stochastic resonance, wherein the combined effect of background noise and weak periodic forcing control the state of the system Benzi et al. (1981, 1982); Gammaitoni et al. (1998). Indeed, because of the compelling consequences of such resonances, there are many methods that have been developed to calculate the escape rate, ranging from eigenfunction expansions Jung (1989) to path-integrals Smelyanskiy et al. (1999); Lehmann et al. (2000). However, the simplest solution used to study the principle characteristics of stochastic resonance appeals to the approximation of the adiabatic limit McNamara and Wiesenfeld (1989).

Recently, concepts of stochastic resonance have been utilized in numerous fields including sensory biology (e.g., Vázquez-Rodríguez et al., 2017; Itzcovich et al., 2017), image processing (e.g., Singh et al., 2016; Gupta and Jha, 2016), signal detection and processing (e.g. Han et al., 2016; Lai and Leng, 2016), and energy harvesting (e.g., Zhang et al., 2016; Kim et al., 2018). The broad impact of stochastic resonance is often viewed as counter-intuitive because rather than background noise obscuring the detection of a weak signal, it leads instead to an enhancement of that signal.

The canonical configuration of stochastic resonance focuses on particles in a double-well potential influenced by white noise and weak periodic forcing. Although there are multiple time-scales involved in the dynamics, the principal interest concerns the time it takes for a particle to transition from one stable point in the potential to the other. Hence, it is common to consider a two-state model using a master equation that describes the time-evolution of the probability density of two discrete states and their exchange rates McNamara and Wiesenfeld (1989). Thus one uses the classical escape rate from one metastable point; when the rate is independent of the slowly varying phase of the periodic forcing, this is called the adiabatic limit. However, considering the extent of the fields in which stochastic resonance plays a role, from climate to engineering to biology Wiesenfeld and Moss (1995), it is of interest to go beyond the adiabatic limit. Here we address the nonadiabatic situation, in the two-state framework, to determine the escape rate when the phase of the periodic forcing does not vary slowly. To achieve this we introduce an asymptotic method to obtain an explicit expression for the escape rate, and in so doing we show how to transform the original double-well potential problem into the two-state model.

II Escape rate under periodic forcing

First, we consider the escape rate from the metastable region of a potential $U(x)$ under the influence of weak periodic forcing $A\text{cos}(\omega t)$ and noise induced fluctuations, as shown in Fig. 1. The Fokker-Planck equation for this situation is

[TABLE]

with boundary conditions $P(\tilde{x}=-\infty,\tilde{t})=P(\tilde{x}=x_{max}+\chi,\tilde{t})=0$, wherein the tilde’s denote dimensional variables. We assume that the magnitude of the noise, $\sigma$, and the amplitude of the periodic forcing, $A$, are both small (in the precise sense outlined below), and that $\chi=O(\Delta x)$. The potential $U$ has the characteristic “diffusivity” scale, $\Delta U=U_{max}-U_{min}$, and length scale, $\Delta x=x_{max}-x_{min}$, which leads to the three parameters

[TABLE]

Using the dimensionless variables $x=(\tilde{x}-x_{max})/\Delta x$, $t=\frac{\Delta U}{(\Delta x)^{2}}\tilde{t}$, $U=(\tilde{U}-U_{max})/\Delta U$ and $J=\tilde{J}/(\omega\Delta x)$, Eq. (1) becomes

[TABLE]

with boundary conditions $P(x=-\infty,t)=P(x=\chi/\Delta x,t)=0$, and the local minimum and maximum are at $x=-1$ and $x=0$ respectively.

The underlying scaling assumptions are that $\epsilon\ll 1$ and $r=O(1)$. The small magnitude of $\epsilon$ is associated with the kinetic energy of a particle near the minimum of the potential being much less than the potential energy, $\Delta U$, necessary to escape from it. The assumption that the external forcing, $A\text{cos}(\omega t)$, is weak relative to the thermal noise is embodied by $r\sim O(\epsilon^{2})$. The dimensionless frequency, $\Omega$, is the ratio of the time-scale for a particle to reach a quasi-stationary state near the potential minimum, $(\Delta x)^{2}/\Delta U$, to the oscillation time-scale $\omega^{-1}$ of the potential. The assumption that $\Omega\ll 1$ implies that the period of the external forcing is much longer than the time required for a particle to reach a quasi-stationary state near the minimum. This assumption facilitates the asymptotic matching procedure near the maximum.

II.1 Potential Minimum – Boundary Layer 1 (BL1): $|x+1|\sim O(\epsilon)$

Near $x=-1$, the potential $U(x)$ can be approximated as $U\simeq-1+\frac{1}{2}a(x+1)^{2}$, where $a\equiv|U^{\prime\prime}_{min}|$. Thus, we rewrite Eq. (3) in terms of a state variable $P_{B1}(\eta,t)$ that depends on the stretched coordinate $\eta=(x+1)/\epsilon$:

[TABLE]

The leading-order solution, written in terms of the original position variable, is

[TABLE]

where $n_{1}$, which will be determined as part of the matching procedure, is a slowly-varying function of time satisfying $\frac{1}{n_{1}}\frac{dn_{1}}{dt}\ll 1$. Implicit in this solution is therefore that the probability density reaches a quasi-steady state around the potential minimum, which arises because the weak noise in system drives only a small leakage of probability across the barrier at the maximum.

II.2 Potential Maximum – Boundary Layer 2 (BL2): $|x|\sim O(\epsilon)$

Near the maximum ($x=x_{max}$) we let $x=\epsilon\zeta$ and the potential $U$ can be approximated as $U\simeq-\frac{1}{2}b\epsilon^{2}\zeta^{2}$, where $b\equiv|U^{\prime\prime}_{max}|$. Hence, Eq. (3) is rewritten in terms of the state variable $P_{B2}(\zeta,t)$ as

[TABLE]

Equation (6) must be solved subject to $P_{B2}(\zeta=\infty,t)=0$ and that the solution match to that for the interior region between the extrema of the potential, outlined presently. The match, however, implies that the solution for the interior delivers a probability flux to the boundary layer around the maximum that varies periodically in time with frequency $\Omega$. This precludes a straightforward solution of Eq. (6) if $\Omega=O(1)$.

Instead, we avoid solving the boundary-layer problem for BL2 for general $\Omega$, and adopt the convenient approximation that the oscillation frequency is small, $\Omega\ll 1$. This allows us to neglect the left-hand side of Eq. (6) and write the quasi-stationary approximation,

[TABLE]

II.3 Interior Region

Within the interior region between the two extrema in the potential, the rapid exponential decline of the probability $P$ suggests that we adopt a WKBJ-type ansatz, viz.,

[TABLE]

from which the leading-order Fokker-Planck equation is

[TABLE]

The characteristic curves of Eq. (9) are given by

[TABLE]

which begin at $x=-1$ when $t\to-\infty$ where $S=S_{min}$, and converge to $x=0$ as $t\to\infty$ where $S=S_{max}$. The solution is

[TABLE]

where

[TABLE]

Note that, in the limit that $\Omega\ll 1$, an integration by parts furnishes the simpler approximation,

[TABLE]

in which we set $\sin\{\Omega[t+T(z)-T(x)]\}\approx\sin\Omega t$, assuming that $\Omega t$ may be $O(1)$ but $\Omega T(x)\ll 1$. This approximation exposes an issue with general solution in Eq. (11): the transit time function $T(x)$ in Eq. (12) diverges logarithmically for $x\to x_{min}$ or $x\to x_{max}$. This does not present a problem at the potential minimum in view of the integration limits in either Eq. (11) or (13), but it does obscure the limit to the maximum. In fact, for $|x|=O(\delta)$ with $1\gg\delta\gg\epsilon$, the interior solution should actually be matched to the boundary-layer solution, leaving $T\sim-b^{-1}\log\delta+O(1)$. Thus we write Eq. (12) as

[TABLE]

where $T_{max}\sim t_{0}-b^{-1}\log\delta$ and $t_{0}$ is an (undetermined) order-one constant time shift that should, in principle, be fixed by a matching argument.

II.4 Asymptotic Matching & Uniform approximation

Asymptotic matching near the local minimum leads to

[TABLE]

and near the local maximum it is required that

[TABLE]

The preceding results suggest an approximation that is valid throughout the two boundary layers and the interior region:

[TABLE]

In the small frequency ($\Omega\ll 1$) approximation, Eq. (17) becomes

[TABLE]

Near $x=\pm 1$, this approximation reduces to the two boundary layer solutions in Eqs. (5) and (7), with $D_{0}$ given by the relevant approximation of Eq. (16), whereas in the interior it reduces to the solution implied by Eq. (13). Note that the limit of the last integral in (18) introduces the approximation $\delta\approx\epsilon$, thereby furnishing a solution that depends on only a single small parameter and avoids any exercise in matching.

II.5 Exit Rate

Now, we construct the exit rate by a suitable integration of the Fokker-Planck equation (3) as follows. Note that

[TABLE]

and because the probability is principally concentrated near the minimum, $x=-1$, we have

[TABLE]

The flux at the origin is given by

[TABLE]

Hence, by combining Eqs. (16) and (17) we obtain

[TABLE]

thereby giving the escape rate $R\equiv\frac{1}{n_{1}}\frac{dn_{1}}{dt}$ as

[TABLE]

In the small frequency ($\Omega\ll 1$) approximation, by substituting Eq. (18) into $J|_{x=0}=\epsilon^{2}\partial P/\partial x|_{x=0}$ we have

[TABLE]

II.6 Adiabatic Limit

In the adiabatic limit, we discard the terms containing factors of $\Omega$ in Eq. (18), to arrive at

[TABLE]

with the associated escape rate

[TABLE]

II.7 The cubic potential

For the cubic potential

[TABLE]

we have

[TABLE]

Hence we have

[TABLE]

[TABLE]

and

[TABLE]

where

[TABLE]

if we again set $\delta=-\epsilon$, and

[TABLE]

captures the suppression of the periodic adiabatic variation of the escape rate by non-adiabatic effects.

In Fig. 2(a) we compare a numerical solution of the Fokker-Planck Equation, (3), with the non-adiabatic (Eqs. 17 and 18) and adiabatic (Eq. 25) analytical solutions. Our numerical method for Eq. (3) is based on the implicit finite difference scheme introduced by Chang and Cooper Chang and Cooper (1970). Both of the non-adiabatic analytical solutions match the numerical solution at the percentage level of accuracy, save for the transition region from the interior to the boundary layer near the maximum ($x=0$). However, the adiabatic solution differs substantially from both of the others, as is particularly evident when $\Omega t=\pi/2$ where the non-adiabatic contribution, $\text{sin}(\Omega t)$, is maximal.

In Fig. 2 (b) we show the time evolution of the escape rates for the four solutions, defined as $J_{x=0}/\int_{-\infty}^{0}Pdx$ for the numerical solution. To calculate the analytical solutions, we use the approximation $n_{1}=1$, which is accurate to machine precision for the parameter settings and times used in the figure. (Likewise, for the double-well potential below, we use the approximation $n_{1}=n_{2}=0.5$ in comparing the asymptotic predictions with numerics in Fig. 6.) The non-adiabatic analytical solutions compare well with the numerical solution, whereas there is a pronounced deviation of the adiabatic solution in both phase and amplitude of the maximum escape rates.

Finally, in Fig. 3 we bring out the deviations of the various approximations as a function of frequency. Note in particular the substantial differences in phase and amplitude between Figs. 3(a) and (c). In particular, while the non-adiabatic analytic solutions compare well with the numerical solution, there is a pronounced deviation of the adiabatic solution in both phase and amplitude of the maximum escape rates.

III Double-well potential and Stochastic Resonance

We now treat Brownian particles in a double-well potential under the influence of weak periodic forcing, which is the original configuration of stochastic resonance Benzi et al. (1981, 1982). By reflection of Fig. 1 we extend the approach described above to construct the approximate solutions in the five regions shown in Fig. 4. The potential $\tilde{U}$ is scaled as before, so that $U=(\tilde{U}-U_{max})/\Delta U$, where $\Delta U$ is now a measure of the height of the barrier, and we define $\Delta x$ as half the distance between the two minima. As the potential may not be symmetrical, this translates to a scaled potential that vanishes at $x=0$ and takes the values $U_{1}$ and $U_{2}$ at the two minima $x_{1}$ and $x_{2}$, respectively. We replace the absorbing boundary near the local maximum with the usual boundary condition, $P(\pm\infty,t)=0$, insuring that the probability is conserved throughout the entire domain as particles move between the two minima.

The asymptotic solution is

[TABLE]

where $a_{\{1,2\}}=U^{\prime\prime}(x_{\{1,2\}})$, and the “constants” of integration, $n_{1}$, $n_{2}$, $S_{1}$, $S_{2}$, $D_{0}$ and $D_{1}$ must be connected by matching the five solutions together. In particular, we find

[TABLE]

[TABLE]

if we again make the approximation that the match can be accomplished at $x=\pm\epsilon$. Again, more compactly we have

[TABLE]

In the small frequency ($\Omega\ll 1$) approximation, Eq. (35) becomes

[TABLE]

Global conservation of probability implies that $\int_{-\infty}^{\infty}P\simeq n_{1}+n_{2}=1$, where $\frac{dn_{1}}{dt}=J|_{x=x_{max}}$, which leads to

[TABLE]

or equivalently

[TABLE]

where

[TABLE]

are the escape rates from $x_{\{1,2\}}$ through $x_{max}$, with

[TABLE]

Equation (38) has the same form as the two-state Master equation used in other approaches to stochastic resonance theory (e.g., McNamara and Wiesenfeld, 1989). The factors $\Upsilon_{\{1,2\}}(\Omega)$ again represent the suppression of the adiabatic variation of the escape rates by non-adiabatic effects; $\Theta_{\{1,2\}}(\Omega)$ represent additional phase shifts. The functions $\Upsilon_{\{1,2\}}(\Omega)$ are equal for a symmetrical base potential, and are illustrated in figure 5 for the quartic potential with $U=\frac{1}{4}x^{4}-\frac{1}{2}x^{2}$. In the theory of stochastic resonance presented in McNamara and Wiesenfeld (1989), the adiabatic variation of the escape rates (given by the replacements $\Upsilon_{\{1,2\}}(\Omega)\to 1$ and $\Theta_{\{1,2\}}(\Omega)\to 0$ in (39)) is responsible for the characteristic improvement in the signal-to-noise ratio. Thus, the suppression factors $\Upsilon_{\{1,2\}}(\Omega)$ determine the destructive effects of non-adiabaticity in a generalization of that analysis. In view of our small-frequency approximation, the exponential decline of $\Upsilon_{1}=\Upsilon_{2}$ for large $\Omega$ seen in figure 5 is again only approximate.

In Fig. 6 we compare the numerical and analytic solutions from Eqs. (35) and (36), for the double well potential $U=\frac{1}{4}x^{4}-\frac{1}{2}x^{2}$. Because our approach is a reflection of the escape calculation described in the previous section, we expect the match to be good, save for the same transition region, here between the two interior regions and the central maximum.

IV Conclusion

We have developed an asymptotic method of calculating the probability density function and the associated escape rate of Brownian particles from a metastable state under weak periodic forcing. The approach uses boundary layers near the two extremes, where the potential $U(x)$ is approximately quadratic and the time-dependent linear Fokker-Planck equations can be solved. In the interior layer separating these, an advection-dominated solution is constructed and the three approximate solutions are matched. Because the evolution of the total probability is equal to the probability flux at the absorbing boundary, we can integrate Fokker-Planck Equation over the complete domain and determine the escape rate in the non-adiabatic limit. Finally, by reflection we extended this asymptotic approach to the problem of a double-well potential with weak periodic forcing to find a solution to the problem of stochastic resonance in the non-adiabatic case. In particular, the ease with which Eq. 39 can be used, and its limits understood through Fig. 5, provide substantial applicability. Given the ubiquity of stochastic resonance, this result is likely of the broadest relevance. Additionally, the approach we take here is complimentary to other general approaches, which focus on the universality of fast-slow systems in stochastic resonance and two state systems Bergland and Gentz (2005, 2006); Lim et al. (2019).

Acknowledgements.

WM and JSW acknowledge the support of Swedish Research Council grant no. 638-2013-9243. WM acknowledges a Herchel-Smith postdoctoral fellowship and JSW a Royal Society Wolfson Research Merit Award for support.