Nonadiabatic escape and stochastic resonance

The escape of particles from a metastable state under the influence of noise is a classical problem in non-equilibrium statistical mechanics [1]. Calculating the rate of escape can be approached in a variety of ways [see e.g. 2–5, 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 [6, 7]. 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 [8–10]. 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 [11] to path-integrals [12, 13]. However, the simplest solution used to study the principle characteristics of stochastic resonance appeals to the approximation of the adiabatic limit [14].

Recently, concepts of stochastic resonance have been utilized in numerous fields including sensory biology [e.g. 15, 16], image processing [e.g. 17, 18], signal detection and processing [e.g. 19, 20], and energy harvesting [e.g. 21, 22]. 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 [14]. 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 [23], 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.

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\frac{\partial P}{\partial \tilde{t}} = \frac{\partial \tilde{J}}{\partial \tilde{x}} \qquad {\rm where} \nonumber \\ &\tilde{J}=\left[\frac{{\rm d}\tilde{U}}{{\rm d}\tilde{x}}-A\cos(\omega \tilde{t})\right]P+\sigma^2\frac{\partial P}{\partial \tilde{x}}, \label{eq:FPE} \nonumber \end{align} \tag{ 1 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \epsilon=\frac{\sigma}{\sqrt{\Delta U}}, \qquad r=\frac{A\Delta x}{\sigma^2}, \qquad {\rm and} \qquad \Omega = \frac{\omega(\Delta x)^2}{\Delta U}. \nonumber \end{align} \tag{ 2 }$$

Zoom In Zoom Out Reset image size

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)$, equation (1) becomes

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial P}{\partial t}= \frac{\partial J}{\partial x}= \frac{\partial}{\partial x} \left[\frac{{\rm d}U}{{\rm d}x} - \epsilon^2 r\cos(\Omega t)\right]P+\epsilon^2\frac{\partial^2 P}{\partial x^2} , \label{eqn:main02} \nonumber \end{align} \tag{ 3 }$$

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 $\newcommand{\e}{{\rm e}} \epsilon\ll1$ and $r=O(1)$. The small magnitude of $\newcommand{\e}{{\rm e}} \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{{\rm cos}}(\omega t)$, is weak relative to the thermal noise is embodied by $\newcommand{\e}{{\rm e}} 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.

2.1. Potential minimum—boundary layer 1 (BL1): ${\newcommand{\e}{{\rm e}} |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 $\newcommand{\e}{{\rm e}} a \equiv |U''_{\min}|$. Thus, we rewrite equation (3) in terms of a state variable $\newcommand{\e}{{\rm e}} P_{B1}(\eta,t)$ that depends on the stretched coordinate $\newcommand{\e}{{\rm e}} \eta = (x+1)/\epsilon$:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial P_{B1}}{\partial t}=\frac{\partial}{\partial\eta}(a\eta P_{B1}) - \epsilon r \cos(\Omega t)\frac{\partial P_{B1}}{\partial\eta} +\frac{\partial^2P_{B1}}{\partial\eta^2}. \nonumber \end{align} \tag{ 4 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{B1}(\eta,t)=n_1\sqrt\frac{{a}}{{2\pi\epsilon^2}}{{\rm exp}} \left[-\frac{a}{2\epsilon^2}(x+1)^2\right], \label{eq:Pmin} \nonumber \end{align} \tag{ 5 }$$

where n₁, which will be determined as part of the matching procedure, is a slowly-varying function of time satisfying $\frac{1}{n_1}\frac{{\rm d}n_1}{{\rm d}t} \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.

2.2. Potential maximum—boundary Layer 2 (BL2): ${\newcommand{\e}{{\rm e}} |x| \sim O(\epsilon)}$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \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(\epsilon). \label{eq:max} \nonumber \end{align} \tag{ 6 }$$

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, as 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 equation (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\ll1$. This allows us to neglect the left-hand side of equation (6) and write the quasi-stationary approximation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{B2} \approx D_0\left[1-\sqrt{\frac{b}{2\pi\epsilon^2}} \int_{-\infty}^{x} {{\rm exp}}\left(-\frac{bz^2}{2\epsilon^2}\right){\rm d}z\right] {{\rm exp}} \left(\frac{b x^2}{2\epsilon^2}\right). \label{eq:pstarb2} \nonumber \end{align} \tag{ 7 }$$

2.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.,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P \sim {{\rm exp}}\left[-\epsilon^{-2}U+S(x,t)\right], \nonumber \end{align} \tag{ 8 }$$

from which the leading-order Fokker–Planck equation is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{\partial S}{\partial t} + \frac{{\rm d}U}{{\rm d}x}\frac{\partial S}{\partial x} = r \frac{{\rm d}U}{{\rm d}x}{{\rm cos}}(\Omega t). \label{eq:S} \nonumber \end{align} \tag{ 9 }$$

The characteristic curves of equation (9) are given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &\frac{{\rm d}x}{{\rm d}t} = \frac{{\rm d}U}{{\rm d}x} \equiv U'(x) \qquad {\rm and}\nonumber \\ &\frac{{\rm d}S}{{\rm d}x}= r {{\rm cos}}(\Omega t), \nonumber \end{align} \tag{ 10 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S=S_{\min}+ r \int_{-1}^{x}{{\rm cos}}\{\Omega\left[t+T(z)-T(x)\right]\}{\rm d}z, \label{eq:ss} \nonumber \end{align} \tag{ 11 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T(x)=\int\frac{{\rm d}x}{U'(x)}. \label{eq:T} \nonumber \end{align} \tag{ 12 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S \simeq S_{\min}+r(x+1){{\rm cos}}(\Omega t) + r \Omega \sin(\Omega t) \int_{-1}^{x}(z+1)T'(z){\rm d}z, \label{eq:sx} \nonumber \end{align} \tag{ 13 }$$

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)\ll1$. This approximation exposes an issue with general solution in equation (11): the transit time function $T(x)$ in equation (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 equation (11) or (13), but it does obscure the limit to the maximum. In fact, for $|x|=O(\delta)$ with $\newcommand{\e}{{\rm e}} 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 equation (11) as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_{\max} \sim S_{\min}+ r \int_{-1}^{\delta}{{\rm cos}}\{\Omega\left[t+T(z)-T_{\max}\right]\}{\rm d}z, \label{eq:ssa} \nonumber \end{align} \tag{ 14 }$$

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

2.4. Asymptotic matching & uniform approximation

Asymptotic matching near the local minimum leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_{\min} \sim {{\rm log}}\left[P_{B1} {\rm e}^{U/\epsilon^2} \right]_{\eta\to\infty} \simeq {{\rm log}}\left(n_1\sqrt{\frac{a}{2\pi\epsilon^2}}\right) - \frac{1}{\epsilon^2}, \label{eq:s} \nonumber \end{align} \tag{ 15 }$$

and near the local maximum it is required that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_{\max} \sim {{\rm log}}\left[P_{B2} {\rm e}^{U/\epsilon^2} \right]_{\zeta\to-\infty} \simeq {{\rm log}}D_0 . \label{eq:logD} \nonumber \end{align} \tag{ 16 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P \approx & n_1\sqrt{\frac{a}{2\pi\epsilon^2}} \left[1-\sqrt{\frac{b}{2\pi\epsilon^2}} \int_{-\infty}^{x} {\rm e}^{-b z^2 / 2 \epsilon^2} {\rm d}z\right] \nonumber \\ &\times{{\rm exp}}\left[-\frac{(1+U)}{\epsilon^2} + r \int_{-1}^{{\rm min}(x,-\epsilon)}{{\rm cos}}\{\Omega\left[t+T(z)-T(x)\right]\}{\rm d}z\right]. \label{uni_ori} \nonumber \end{align} \tag{ 17 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P \approx & n_1\sqrt{\frac{a}{2\pi\epsilon^2}} \left[1-\sqrt{\frac{b}{2\pi\epsilon^2}} \int_{-\infty}^{x} {\rm e}^{-b z^2 / 2 \epsilon^2} {\rm d}z\right] \nonumber \\ &\times{{\rm exp}}\left[-\frac{(1+U)}{\epsilon^2} + r [{\rm min}(x,-\epsilon)+1]{{\rm cos}}(\Omega t) + r \Omega \sin (\Omega t) \int_{-1}^{{\rm min}(x,-\epsilon)} (z+1) T'(z){\rm d}z\right]. \label{unicorn} \nonumber \end{align} \tag{ 18 }$$

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

2.5. Exit rate

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}}{{\rm d}t}\int_{-\infty}^{0}P{\rm d}x = \int_{-\infty}^{0}\frac{\partial J}{\partial x} {\rm d}x = J|_{x=0}, \label{eq:flux} \nonumber \end{align} \tag{ 19 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{-\infty}^{0}P{\rm d}x \simeq \epsilon \int_{-\infty}^{\infty}P_{B1}{\rm d}\eta = n_1. \nonumber \end{align} \tag{ 20 }$$

The flux at the origin is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle J|_{x=0} \simeq \left. \epsilon^2\frac{\partial P}{\partial x} \right|_{x=0} =-\epsilon^2D_0\sqrt{\frac{b}{2\pi\epsilon^2}}. \nonumber \end{align} \tag{ 21 }$$

Hence, by combining equations (16) and (17) we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}n_1}{{\rm d}t} = -\frac{\sqrt{ab}}{2\pi} {{\rm exp}}\left(-\frac{1}{\epsilon^2}+ r \int_{-1}^{\delta}{{\rm cos}}\{\Omega[t+T(z)-T_{\max}]\}{\rm d}z\right)n_1, \label{eq:n1_evolution} \nonumber \end{align} \tag{ 22 }$$

thereby giving the escape rate $\newcommand{\e}{{\rm e}} R \equiv \frac{1}{n_1}\frac{{\rm d}n_1}{{\rm d}t}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R = \frac{\sqrt{ab}}{2\pi} {{\rm exp}}\left(-\frac{1}{\epsilon^2} + r \int_{-1}^{\delta}{{\rm cos}}\{\Omega[t+T(z)-T_{\max}]\}{\rm d}z\right) . \label{eq:Rapprox} \nonumber \end{align} \tag{ 23 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R = \frac{\sqrt{ab}}{2\pi} {{\rm exp}}\left(-\frac{1}{\epsilon^2}+r{{\rm cos}}\Omega t+r\Omega{{\rm sin}}\Omega t\int_{-1}^{\delta}(z+1)T'(z){\rm d}z\right). \label{eq:Rsmallomega} \nonumber \end{align} \tag{ 24 }$$

2.6. Adiabatic limit

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P_{ad} = n_1\sqrt{\frac{a}{2\pi\epsilon^2}} \left[1-\sqrt{\frac{b}{2\pi\epsilon^2}} \int_{-\infty}^{x} {\rm e}^{-b z^2 / 2 \epsilon^2} {\rm d}z\right] {{\rm exp}}\left[-\frac{(1+U)}{\epsilon^2} + r(x+1){{\rm cos}}(\Omega t) \right], \label{unicorna} \nonumber \end{align} \tag{ 25 }$$

with the associated escape rate

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R_{ad} = \frac{\sqrt{ab}}{2\pi}{{\rm exp}}\left[r \cos(\Omega t) - \frac{1}{\epsilon^2}\right]. \label{eq:Rad} \nonumber \end{align} \tag{ 26 }$$

2.7. The cubic potential

For the cubic potential

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \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], \nonumber \end{align} \tag{ 27 }$$

we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle U=-3x^2-2x^3, \qquad \frac{{\rm d}U}{{\rm d}x}=-6x(x+1), \qquad {\rm and} \qquad a=b=6. \nonumber \end{align} \tag{ 28 }$$

Hence we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle T(x) = \frac{1}{6}{{\rm log}}\left(\frac{1+x}{-x}\right), \nonumber \end{align} \tag{ 29 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S=S_{\min}+ \frac{3r}{2} \int_{-\infty}^{T(x)}\frac{{{\rm cos}}\{\Omega[\tau+t-T(x)]\}} {{{\rm cosh}}^2(3\tau)}{\rm d}\tau, \nonumber \end{align} \tag{ 30 }$$

and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R = \frac{3}{\pi}\exp\left[-\frac{1}{\epsilon^2} + r \Upsilon(\Omega) \cos\theta\right], \label{eq:nonadR} \nonumber \end{align} \tag{ 31 }$$

where

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \theta \equiv \Omega(t-T_{\max}) \simeq \Omega \left(t+ \frac{1}{6}\log\epsilon \right), \nonumber \end{align*}$$

so that if we again set $\newcommand{\e}{{\rm e}} \delta=-\epsilon$, for which

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Upsilon(\Omega) = \frac{\frac{1}{6}\Omega\pi}{\sinh(\frac{1}{6}\Omega\pi)} \label{Uppy} \nonumber \end{align} \tag{ 32 }$$

then we can capture the suppression of the periodic adiabatic variation of the escape rate by non-adiabatic effects.

In figure 2(a) we compare a numerical solution of the Fokker–Planck equation, (3), with the non-adiabatic (equations (17) and (18)) and adiabatic (equation (25)) analytical solutions. Our numerical method for equation (3) is based on the implicit finite difference scheme introduced by Chang and Cooper [24]. 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, ${{\rm sin}}(\Omega t)$, is maximal.

Zoom In Zoom Out Reset image size

In figure 2(b) we show the time evolution of the escape rates for the four solutions, defined as $J_{x=0}/\int_{-\infty}^0 P {\rm d}x$ for the numerical solution. To calculate the analytical solutions, we use the approximation n₁  =  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 figure 6.) In both amplitude and phase, the non-adiabatic solution in (31) compares more favourably with the numerical solution than the adiabatic result in (24).

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

Zoom In Zoom Out Reset image size

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 [8, 9]. By reflection of figure 1 we extend the approach described above to construct the approximate solutions in the five regions shown in figure 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₁ and U₂ at the two minima x₁ and x₂, 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.

Zoom In Zoom Out Reset image size

The asymptotic solution is

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle P = \left\{\begin{array}{@{}ll@{}} n_1 \sqrt{\frac{a_1}{2\pi\epsilon^2}} \exp \left[- \frac{a_1}{2\epsilon^2}(x-x_1)^2\right], & {{\rm BL}_{1}} \nonumber \\ \exp\left(S_1-\frac{U}{\epsilon^2} +r\int_{-1}^{x}\cos[\Omega(t+T(z)-T(x))]{\rm d}z\right), & {{\rm Interior} \ 1} \nonumber \\ \left[D_0+D_1\sqrt{\frac{b}{2\pi\epsilon^2}}\int_{-\infty}^{x} \exp\left(-\frac{bz^2}{2\epsilon^2}\right){\rm d}z \right] \exp \left(\frac{bx^2}{2\epsilon^2}\right), & {{\rm BL}_{\max}} \nonumber \\ \exp\left(S_2-\frac{U}{\epsilon^2} +r\int_{1}^{x}\cos[\Omega(t+T(z)-T(x))]{\rm d}z\right), & {{\rm Interior} \ 2} \nonumber \\ n_2 \sqrt{\frac{a_2}{2\pi\epsilon^2}} \exp \left[- \frac{a_2}{2\epsilon^2}(x-x_2)^2\right], & {{\rm BL}_{2}} \end{array}\right. \nonumber \end{align*}$$

where $a_{\{1,2\}}=U''(x_{\{1,2\}})$, and the 'constants' of integration, n₁, n₂, S₁, S₂, D₀ and D₁ must be connected by matching the five solutions together. In particular, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle S_{\{1,2\}} = \log\left(n_{\{1,2\}}\sqrt{\frac{a_{\{1,2\}}}{2\pi\epsilon^2}}\right) +\frac{U_{\{1,2\}}}{\epsilon^2}, \nonumber \end{align} \tag{ 33 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle &D_0={{\rm exp}}\left(S_1+r \int_{-1}^{-\epsilon}\cos[\Omega(t+T(x)-T(-\epsilon))]{\rm d}x \right) \qquad {{\rm and}} \nonumber \\ &D_0+D_1={{\rm exp}}\left(S_2+r \int_{1}^{\epsilon}\cos[\Omega(t+T(x)-T(\epsilon))]{\rm d}x\right), \nonumber \end{align} \tag{ 34 }$$

if we again make the approximation that the match can be accomplished at $\newcommand{\e}{{\rm e}} x=\pm\epsilon$. More compactly we have

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P &\simeq n_1\sqrt\frac{{a_1}}{{2\pi\epsilon^2}}{{\rm exp}}\left(-\frac{U}{\epsilon^2} +r\int_{-1}^{{{\rm min}}(x,-\epsilon)}{{\rm cos}}[\Omega(t+T(x')-T({{\rm min}}(x,-\epsilon)))]{\rm d}x'\right) \nonumber \\ &+\left[n_2\frac{\sqrt{a_2 b}}{2\pi\epsilon^2}{{\rm exp}}\left(-\frac{U}{\epsilon^2} +r\int_{1}^{{{\rm max}}(x,\epsilon)}{{\rm cos}}[\Omega(t+T(x')-T({{\rm max}}(x,\epsilon)))]{\rm d}x'\right) \right. \nonumber \\ &\left. -n_1\frac{\sqrt{a_1 b}}{2\pi\epsilon^2}{{\rm exp}}\left(-\frac{U}{\epsilon^2} +r\int_{-1}^{-\epsilon}{{\rm cos}}[\Omega(t+T(x')-T(-\epsilon))]{\rm d}x'\right) \right] \int_{-\infty}^{x}{{\rm exp}}\left(-\frac{b}{2\epsilon^2}x'^2\right){\rm d}x'. \label{uni_ori_double} \nonumber \end{align} \tag{ 35 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle P &\simeq n_1\sqrt\frac{{a_1}}{{2\pi\epsilon^2}}{{\rm exp}}\left(-\frac{U+1}{\epsilon^2} +r(x+1){{\rm cos}}(\Omega t)+r\Omega{{\rm sin}}(\Omega t)\int_{-1}^{{{\rm min}}(x,-\epsilon)} (z+1)T'(z){\rm d}z\right) \nonumber \\ &+\left[n_2\frac{\sqrt{a_2 b}}{2\pi\epsilon^2}{{\rm exp}}\left(-\frac{U+1}{\epsilon^2} +r(x-1){{\rm cos}}(\Omega t) + r\Omega{{\rm sin}}(\Omega t)\int_{1}^{{{\rm max}}(x,\epsilon)}(z-1)T'(z){\rm d}z\right) \right. \nonumber \\ &\left. -n_1\frac{\sqrt{a_1 b}}{2\pi\epsilon^2}{{\rm exp}}\left(-\frac{U+1}{\epsilon^2} +r(1-\epsilon){{\rm cos}}(\Omega t)+r\Omega{{\rm sin}}(\Omega t)\int_{-1}^{-\epsilon}(z+1)T'(z){\rm d}z\right)\right] \int_{-\infty}^{x}{{\rm exp}}\left(-\frac{b}{2\epsilon^2}x'^2\right){\rm d}x'. \label{uni_ori_smallomega} \nonumber \end{align} \tag{ 36 }$$

Global conservation of probability implies that $\int_{-\infty}^{\infty}P \simeq n_1+n_2 = 1$, where $\frac{{\rm d}n_1}{{\rm d}t}=J|_{x=x_{\max}}$, which leads to

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}n_1}{{\rm d}t} = \epsilon^2D_1\sqrt{\frac{b}{2\pi\epsilon^2}}, \nonumber \end{align} \tag{ 37 }$$

or equivalently

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm d}n_1}{{\rm d}t}=-R_1n_1+R_2n_2=R_2-(R_1+R_2)n_1, \label{SM} \nonumber \end{align} \tag{ 38 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle R_{\{1,2\}} = \frac{\sqrt{a_{\{1,2\}}b}}{2\pi} \exp\left(-\frac{U_{\{1,2\}}}{\epsilon^2}\right) \exp\left(r \Upsilon_{\{1,2\}} \cos [\Omega(t-T(\{-\epsilon,\epsilon\})) + \Theta_{\{1,2\}}] \right) \label{escape2} \nonumber \end{align} \tag{ 39 }$$

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

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{\{-1,1\}}^{\{-\epsilon,\epsilon\}} \cos[\Omega(t+T(x)-T(\{-\epsilon,\epsilon\}))] {\rm d}x \equiv \Upsilon_{\{1,2\}}(\Omega) \cos [\Omega(t-T(\{-\epsilon,\epsilon\})) + \Theta_{\{1,2\}}(\Omega)]. \nonumber \end{align} \tag{ 40 }$$

Equation (38) has the same form as the two-state Master equation used in other approaches to stochastic resonance theory [e.g. 14]. 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 [14], the adiabatic variation of the escape rates (given by the replacements $\Upsilon_{\{1,2\}}(\Omega)\to1$ and $\Theta_{\{1,2\}}(\Omega)\to0$ 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.

Zoom In Zoom Out Reset image size

In figure 6 we compare the numerical and analytic solutions from equations (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.

Zoom In Zoom Out Reset image size

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 the 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 equation (39) can be used, and its limits understood through figure 5, provide substantial applicability. Given the ubiquity of stochastic resonance, this result is likely of the broadest relevance. Finally, 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 [25–27].

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. Part of this work was performed at the Geophysical Fluid Dynamics Summer Study Program at the Woods Hole Oceanographic Institution, which is supported by the National Science Foundation and the Office of Naval Research under No. OCE-1332750.