Stochastic resonance for nonequilibrium systems

Stochastic resonance (SR) is a prominent phenomenon in many natural and engineered noisy systems, whereby the response to a periodic forcing is greatly ampliﬁed when the intensity of the noise is tuned to within a speciﬁc range of values. We propose here a general mathematical framework based on large deviation theory and, speciﬁcally, on the theory of quasipotentials, for describing SR in noisy N -dimensional nonequilibrium systems possessing two metastable states and undergoing a periodically modulated forcing. The drift and the volatility ﬁelds of the equations of motion can be fairly general, and the competing attractors of the deterministic dynamics and the edge state living on the basin boundary can, in principle, feature chaotic dynamics. Similarly, the perturbation ﬁeld of the forcing can be fairly general. Our approach is able to recover as special cases the classical results previously presented in the literature for systems obeying detailed balance and allows for expressing the parameters describing SR and the statistics of residence times in the two-state approximation in terms of the unperturbed drift ﬁeld, the volatility ﬁeld, and the perturbation ﬁeld. We clarify which speciﬁc properties of the forcing are relevant for amplifying or suppressing SR in a system and classify forcings according to classes of equivalence. Our results indicate a route for a detailed understanding of SR in rather general systems.


I. INTRODUCTION
Stochastic resonance (SR) is a rather special and somewhat counterintuitive mechanism where noise plays the constructive role of catalyzing the amplification of the response of a system to a weak periodic signal. SR was originally proposed independently by Benzi et al. [1][2][3] and by Nicolis [4,5] as a way to explain the occurrence of periodically spaced ice ages in the Quaternary era despite the presence of extremely weak periodic modulations of the incoming solar radiation due to the Milankovich cycles. Since then, SR has been found and studied in a myriad of natural and engineered systems and has been thoroughly explored through theory, experiments, and numerical simulations. We report examples from laser systems [6], atomic physics [7], nanostructures [8], optics [9], control theory [10], circuits [11], ecology [12], geosciences [13,14], biology [15], physiology [16], neurosciences [17], and psychology [18], among others. Many valuable reviews of the topic are available [19][20][21][22]. Additionally, SR has become the subject of careful mathematical investigations; see, e.g., Refs. [23][24][25][26][27].
The mathematical archetype for SR is the system whose (overdamped) dynamics is described by the following * v.lucarini@reading.ac.uk where V (x) = −ax 2 + bx 4 (a, b > 0) features two stable equilibria at x = ±x 0 = ±(a/2b) 1/2 , ω is the frequency of the periodic forcing, controls the intensity of the periodic forcing, dW is the increment of a Brownian motion, and σ modulates the intensity of the stochastic forcing. Stochastic forcing will lead to trajectories performing transitions between the basin of attractions of the two stable equilibria; we indicate by 2 (1) the basin of attraction of x 0 (−x 0 ). If = 0 and σ > 0, the classical Kramers' theory [28] predicts that, in the weak-noise limit, the average transition rate between the two basins of attraction is The Kramers' formula has been generalized by by Bovier et al. [29] for N-dimensional gradient flows of the form dx(t ) = −∇U (x)dt + σ dW , where x ∈ R N and dW is a vector whose components are N independent increments of a Brownian motion; see also Ref. [30]. The classical result on SR says that, by and large, if we now switch on the periodic forcing by setting > 0, one gets that the periodic component of the expectation value x(t ) is greatly amplified if r 1,σ = r 2,σ ≈ ω/π = 2/T . By tuning the noise in this way, one can obtain a virtually perfect synchronization between the periodic forcing and the transitions between the two basins of attraction for each individual ensemble member; see also Ref. [31]. The problem has been later generalized to the case of more general noise laws [32][33][34][35], while a general treatment of SR in an asymmetric potential with complex stochastic forcing has been presented by Qiao et al. [36]. A useful simplification of the problem is obtained by performing a maximal coarse-graining procedure such that the system is reduced to two discrete states, which corresponding to the stable equilibria in the continuum description [19,37]. The coarse graining leads to concentrating on the interwell hopping and to neglecting the intrawell dynamics. Solid mathematical foundations to this approach can be found in, e.g., Refs. [38,39]. We will refer to these states as x 1 and x 2 , corresponding to the two basin of attractions of x 0 and −x 0 , respectively. The analysis of SR for the two-state model has been presented for the symmetric case by, e.g., Refs. [19,37], while general results have been presented for the asymmetric case and for non-Gaussian stochastic forcing in Refs. [40,41]. We will come back to these results later in the paper.
Most of the results on SR have been derived in the case of one-dimensional systems or, more generally, of Ndimensional gradient flows. The goal of this paper is propose a general formulation of SR able to encompass the case of nonequilibrium systems possessing two metastable states. The classical results valid for system obeying detailed balance are obtained as special cases. We will then consider Ndimensional SDEs with a fairly general class of noise laws and assume that, if the noise is switched off, the deterministic dynamics features two competing asymptotic states. We will rewrite some of the classical results of SR-within the framework of the two-state approximation-in terms of quantities that can be derived from the equations of motion. Our treatment will in principle include the case of stochastically perturbed systems featuring, when noise is removed, two competing chaotic attractors supported on strange sets. The edge state embedded in the basin boundary [42][43][44][45][46] can as well, in principle, be supported on a strange set and feature chaotic dynamics. While some hints at SR for this general case have been proposed in the literature [20], a complete treatment has not been yet presented, to the best of the author's knowledge.
Note that different mechanisms of SR-like phenomena for chaotic systems have been discussed in the literature, where deterministic chaos plays the role of internally generated noise, and no external stochastic forcing is needed. Resonant response to periodic forcing has been found in the case of systems inhabiting preferentially two distinct preferred regions of phase space [47,48], or in the case a parameter periodically fluctuates below and above the value determining the onset of chaotic motions [49].
The paper is structured as follows. In Sec. II we present the mathematical framework we deem useful for studying N-dimensional noisy nonequilibrium systems. In Sec. III we provide a derivation of SR conditions based upon a linear response approach, and we study the resonant amplification of the system's response to a periodic forcing of a very general nature. We draw our conclusion, discuss the limitations of the present work, and present possible future lines of research in Sec. IV. The Appendix contains results pertaining to the statistics of residence times, i.e., the time intervals spent consecutively in each state before a noise-induced transition takes place.

II. MATHEMATICAL FRAMEWORK
We consider an SDE in Itô form written as where x, F ∈ R N , dW j is the increment of an N-dimensional Brownian motion, C i j (x) = s ik (x)s jk (x) is the noise covariance matrix with s i j (x) ∈ R N×N , and σ 0. We assume thaṫ has multiple steady states, so that the phase space is partitioned between the basins of attraction B j of the attractors j and the boundaries ∂B l , l = 1, . . . , L separating such basins. Orbits initialized on the basin boundaries ∂B l , l = 1, . . . , L are attracted towards invariant saddles. Such saddles l , l = 1, . . . , L are called edge states [42][43][44][45] and can feature chaotic dynamics [46]. In this latter case we refer to the edge states as Melancholia states [46,50,51]. In the absence of noise, the asymptotic state is uniquely determined by the initial condition, while noise allows trajectories to hop across boundaries between the various basins of attraction.

A. Computing the quasipotential
In the case of elliptic (and possibly hypoelliptic) diffusion processes, the Freidlin-Wentzell [52] theory and modifications thereof [53][54][55] show that in the weak-noise limit σ → 0 the (unique) invariant measure can be written as a large deviation law: where the rate function (x) is referred to as quasipotential, and we have neglected the preexponential term. Specifically, the symbol ∼ in Eq. (4) implies that (x) = −2 lim σ →0 σ 2 log σ (x). The function (x) can be obtained as follows. We solve the stationary Fokker-Planck equation corresponding to Eq. (3): where J is the current density, we consider the weak noise limit, and we use as an ansatz the expression given in Eq. (4). In this paper we adopt the Einstein convention on repeated indices. We obtain the following Hamilton-Jacobi equation [56]: The previous equation allows one to express in terms of the drift and volatility fields. The quasipotential can also be computed by solving the variational problem associated with the Freidlin-Wentzell action [57]. Finally, alternative routes for computing have been proposed [58,59]. 1 The explicit computation of is, in general, far from trivial, yet of great interest in many applications; see, e.g., Ref. [61] for the case of biological systems. Brackston et al. [62] have recently proposed an algorithm for estimating in the case the governing equations are polynomial, and it involves solving an optimization over the coefficients of a polynomial function. Alternatively, Tang et al. [63] proposed a variational method for estimating in the the populations corresponding to each determininistic attractor without resorting to computing the invariant measure.
Following Refs. [53,54], we now describe the dynamical meaning of . Indeed, solving the previous Hamilton-Jacobi equation corresponds to the fact that it is possible to write the drift vector field as the sum of two vector fields obeying the orthogonality condition R i (x)∂ i (x) = 0. In the case Eq. (3) describes a thermodynamical system near equilibrium, R defines the time reversible dynamics, while F − R defines the irreversible, dissipative dynamics [64]. One finds that As a result, just as in the case of gradient flow, has the role of a Lyapunov function whose decrease describes the convergence of an orbit to the attractor. Specifically, (x) has local minima at the deterministic attractors j , j = 1, . . . , J and has a saddle behavior at the edge states l , l = 1, . . . , L. If an attractor (or an edge state) is chaotic, has constant value over its support, which can then be a strange set [53,54]. The standard case of gradient flow and noise correlation matrix proportional to the identity [obtained by setting

B. Noise-induced escape from the attractor
The probability that an orbit with initial condition in B j does not escape from it over a time p decays as wherer j,σ is the escape rate and where j = ( l ) − ( j ) is the lowest pseudo-potential barrier height [55], i.e., has the lowest value in l compared to all the other edge states neighboring j . In general, one may need to add a correcting prefactor to P j,σ (p) [55]. Equation (9) defines the residence-time distribution for basin of attraction B j . See the Appendix for further discussion of this key statistical property.
Note thatr j,σ given in Eq. (9) does not contain the preexponential factor, as opposed to Eq. (2). Bouchet and Reygner [65] provided an expression for such a preexponential factor for general nonequilibrium diffusion processes under the assumption that attractors and edge states are simple points, thus generalizing the results by Bovier et al. [29]. As we also aim at treating a more general setting for the geometry of attractors and edge states, below we pay the price of having to take the preexponential factors as phenomenological parameters one can find from experiments or numerical simulations. We also remark that, in the zero-noise limit, the transition paths outside a basin of attraction follow the instantons. Instantons are defined as solutions of that connect a point in j to a point in l . Instantonic trajectories have a reversed component of the gradient contribution to the vector field compared to regular-relaxationtrajectories.

III. STOCHASTIC RESONANCE
Let us now assume that, generalizing Eq. (1), we perturb the Eq. (3) as follows: As a result of the perturbation, the rate function (x) will depend on the parameter . Assuming is small, in the spirit of linear response theory, we can expand indicates higher order terms. Substituting this expansion in Eq. (6) and collecting the first-order terms, we obtain (12) Solving the previous linear equation with respect to (x) allows us to derive the first-order correction to the rate function, so that → + , with the ensuing modifications in, e.g., Eqs. (4) and (9). In the latter, taking again a linear approximation, the quasipotential difference is evaluated by considering the unperturbed attractor and edge state. Bouchet et al. [57] showed in great generality that the correction term can indeed be found, and they proposed an algorithmic procedure to compute the perturbative terms at all orders in .
In the weak-noise limit, 2 the escape rate from the basin of attraction B j is where we have assumed α j /r j 1. We have explicit expressions for the rate in terms of the quasipotential of the system. The escape times implied by the rates in Eq. (13)-(15) are very long compared to the dynamical timescales of the system within basins of attraction. In our case, the prefactors A j , j = 1, 2 should be estimated from numerical experiments performed with different values of the noise 2 We do not take the limit σ → 0 because this leads to concentrating the measure over the deterministic attractor featuring the lower value of the quasipotential, leading to the disappearance of bistability; see Ref. [51] for a discussion of an associated first-order phase transition in a climate model. strength σ . Nonetheless, unless A 1 /A 2 is very different from 1 (which amounts to having radically different properties of the quasipotential near the two attractors), the results below depend weakly (compared to σ ) on A 1 /A 2 .
We now treat the case of a time-dependent variant of the perturbed evolution given in Eq. (11), where we consider → cos(ωt ). This generalizes the standard onedimensional case described in Eq. (1), where = x. If the period of the oscillation T = 2π/ω is much longer than the internal timescales of the system within each attractor, we obtain that, using the quasiadiabatic approximation [19], the escape rates of the perturbed system can be written as We then perform a coarse graining and consider the two-state system corresponding to the two unperturbed attractors 1 and 2 . The master equation for the population of state 1, n 1 (t ), isṅ where n 1 (t ) + n 2 (t ) = 1. The construction of a meaningful master equation relies on the presence of clear timescale separation between the relaxation motions near each attractor and those across the edge state, which depends critically on the presence of weak noise.
Within the two-state approximation, the time-dependent expectation value of given where O j is the expectation value of O in the measure supported on the unperturbed attractor j . In the usual case described by Eq. (1), one typically chooses O = x.
The value of R indicates whether we are in SR conditions or not, because R measures the strength of the periodic motion of the ensemble mean of the trajectories. R tends to zero as σ → 0 and σ → ∞ (keeping in mind that the latter limit goes against our weak noise assumption), and one expects that a maximum for R is achieved for intermediate values of σ . Such a maximum defines conditions of SR. As discussed in Refs. [36,40], the resonance is, ceteris paribus, weakened by the presence of strong asymmetries in the system. Taking the standard symmetric case where A 1 = A 2 and 1 = 2 , one gets, by maximizing R, the following transcendental equation for σ defining the SR condition: The resonance condition we find agrees, obviously, with the result presented in Ref. [19]; the main improvement we get in our result is that we can relate all parameters in the previous equation to the unperturbed equations of motion via . We will comment below on the relevance of the specific functional form of the perturbation field G. Note that, by definition, 1 − 2 = 2 − 1 . A second measure of SR 3 is obtained by studying under which conditions the periodic forcing and the noise interact constructively to create in the power spectrum of a general observable a strong spectral feature at the frequency ω of the periodic forcing. We study the t-averaged correlation function for a general observable O: and, in particular of its symmetrized Fourier transform } is the Fourier transform of C O (τ ) and ν is the angular frequency. 3 Different measures of the quality of the SR based on the synchronization between the periodic forcing and occurrence of transitions between the two basins of attraction have been proposed in Ref. [24].
In order to find the correct expression of S s O (ν) one needs to consider transient behavior as well, as opposed to the case of the estimate of R above. Following the careful calculations in Ref. [40], one finds that where the first two terms refer to the singular components of the spectrum and the second two describe the continuum component. Specifically, one has while the rather convoluted (smooth) function is not reported here for reasons that become apparent below. Note that the zero-frequency component can be removed by redefining Following Ref. [19], one defines the (linear) spectral amplification SNR as follows: which leads to the final result: As we well know, in the small limit, SNR does not depend on ω. Additionally, the parameter SNR is clearly maximized in the symmetric case A 1 = A 2 , 1 = 2 , where we obtain SNR = π/2| 1 − 2 | 2 A 1 exp(−2 1 /σ 2 )/σ 4 . In the symmetric case, SNR is maximized if σ 2 = 1 , regardless of the perturbation field, which generalizes what is given in Ref. [19].
The expression of R in Eq. (21) and of SNR in Eq. (30) indicate that, in the weak-perturbation and weak-noise limit, the choice of the perturbation field G impacts the strength of the signal (both in conditions of SR or not) exclusively through the factor | 1 − 2 |. Clearly, perturbation fields G's differing in the transversal [(with respect to the gradient structure; see Eq. (7)] component have the same effect in terms of SR.
In particular, we find that SR is entirely suppressed if 1 = 2 . In other terms, if the change in the quasipotential is the same in 1 and 2 , there is no SR phenomenon at all. Again, this condition can be realized for a very large class of nontrivial G's. In this case, adiabatically we see an periodic increase and decrease of the both r 1,σ and r 2,σ . This amounts to a slow modulation in the overall interwell timescale of the system and has no differential effects on the transition 1 → 2 and 2 → 1.

IV. CONCLUSIONS
Our analysis in this paper bridges the investigation of the statistical properties of general nonequilibrium systems with that of SR. Our findings give as a special case the classical results valid for systems obeying a detailed balance, as in the case of N-dimensional gradient flows forced by standard additive noise.
Our approach revolves around the computation of the quasipotential of the unperturbed dynamics. The quasipotential is proportional to the rate function defining the large deviation law describing the invariant measure of the system and controls the rate of escape of trajectories from a given basin of attraction through a relevant edge state. Additionally, the quasipotential is associated with the decomposition of the drift term into the sum of a gradient component (defined by the quasipotential) and a transversal component, which are mutually orthogonal. Finally, acts as a Lyapunov function for the dynamics and reaches local minima at the attractor(s) of the system. In the case of multistable system, has a saddle behavior at the edge states embedded in the boundaries of the basins of attraction.
If we consider a system possessing two metastable states undergoing a periodic modulation to the drift term, the previous framework allows us to derive some of the main classical results of SR within the two-state approximation with the great advantage that the parameters contained in the SR conditions can be derived from the drift and volatility terms of rather general equations of motion. Indeed, we need to compute the correction to the quasipotential , the former being associated with the spatial pattern of the periodically driven perturbation field. We are able to deal with the fairly general case where the the attractors and/or the edge state feature chaotic dynamics.
We have clarified that SR is an intrinsic property of the unperturbed system, which is catalyzed by the presence of the periodically oscillating perturbation field. The details of the perturbation to the drift term impact the SR intensity through a simple factor depending on . Additionally, we have explained that one can define classes of equivalence of perturbations in terms of their SR properties, with each class constituted by elements differing only by the transversal component with respect to the gradient structure determined by . Similarly, one finds that (see the Appendix) the statistical properties of the residence times can be fully described in terms of and .
Our results can in principle be extended to the case of multiple metastable states connected through a potentially nontrivial network of channels defined by the edge states between them. Addressing the challenges of a complex topology of transition paths would probably benefit from taking advantage of the sophisticated supersymmetric techniques developed by Tȃnase-Nicola and Kurchan [66] in the context of noisy gradient flows.
A relevant improvement to our results would come from the possibility of using a general formula for the preexponential factor of the escape rates valid also for the case of nontrivial attracting and saddle sets. The lack of such a formula makes our finding somewhat phenomenological, yet hopefully useful. Clearly, the results above would greatly benefit from a more rigorous level of mathematical formulation, which is currently beyond the abilities of the author.
The findings above pave the way for the rigorous study of SR-related phenomena in many systems, especially taking into account that the assumption of a (one-dimensional) gradient structure for the drift component of the flow is far from verified (or meaningful, in fact) in general; see discussion in Ref. [67]. We foresee applications in the many areas where SR has proved to be a valuable and useful concept. The fact that several algorithms for computing the quasipotential have been recently presented in the literature makes possible a systematic exploitation of the results discussed here.
As for the inclination of the author, further investigation will focus on the study of SR in geophysical systems, where multistability is often encountered and of great relevance [68,69]. In the spirit of the first investigations of SR, the author will attempt a careful numerical and analytical examination of SR in the multistable climate model presented in Refs. [50,51], where noise permits transitions between the competing snowball and warm climate states. Furthermore, following Refs. [13,14], the author aims at studying using the formalism presented here the SR acting between the two competing on and off states of the deep ocean circulation, featuring a strong and virtually absent overturning circulation in the Atlantic Ocean, respectively [70,71].

APPENDIX: RESIDENCE-TIME DISTRIBUTION
The residence time for the system in a state is the random variable describing the time interval spent consecutively in such a state before a noise-induced transition occurs and the system reaches another state. The statistics of the residence times is one of the essential properties of stochastically perturbed multistable systems, and the typical exponential residence-time distribution for the system in absence of periodic forcing has been introduced in Eq. (9). When a periodic forcing is applied, the simple exponential law is fundamentally altered, in such a way that escape from a state is periodically enhanced and suppressed.
In the classic setting of Eq. (1), one finds peaks in the residence time in each state (modulated by the exponential decay associated with the time-averaged transition rate) at times p = (2n + 1)T /2, as a result of the alternating increase and decrease of the potential barrier for either state. Key is the fact that the time-dependent anomaly of the potential is at all times in opposite phase in the two states. If the system is in state 1 at time 0, after having performed a transition from state 2 when the potential barrier of state 2 is at a minimum, it has to wait half a period before the potential barrier of state 1 reaches a minimum. This explains the first peak. If the transition does not take place, the system has to wait a full period before reaching again favorable conditions for the jump. The intensity of such an effect depends dramatically on /σ 2 , up to reaching an almost perfect synchronization of the transitions with the phase of the periodic forcing [19]. A very accurate analysis of the residence-time statistics has been given by Choi et al. [72] and Berglund and Gentz [73] in the case of symmetric potential described by Eq. (1).
In this Appendix, we build on previous results by Löfstedt and Coppersmith [74], who studied residence-time statistics within the two-state approximation and considered asymmetries in the two states and in their response to the periodic perturbation. We rewrite the escape rate from the state j = 1, 2 as where r j,σ, =1/T T +t t dt 1 r j,σ, (t 1 ) and T +t t dt 1 δ {r j,σ, (t 1 )} = 0, with T = 2π/ω. We now consider Eq. (13) and substitute → cos(ωt ). By expanding in power series the exponential function, one obtains where F (x) = ∞ j=0 1 2 2 j 1 ( j!) 2 x 2 j . The function F increases rapidly with x, with F (0) = 1 [see the linear expansion of r j,σ, (t ) in Eq. (16)], while the first nontrivial term is given by 1/4x 2 , Therefore, as a result of the periodic modulation of the argument of the exponent, the average rate of escape within a period is larger than the unperturbed one, with the first correction being quadratic in . We then derive Following Löfstedt and Coppersmith [74], we have that the residence-time distribution P j,σ, (p), j = 1, 2 can be written as P j,σ, (p) = N j,σ, exp(−r j,σ, p)G j,σ, (p), where G j,σ, (p) is periodic of period T , i.e., G j,σ, (p + T ) = G j,σ, (p) and N j,σ, is a normalization constant such that ∞ 0 dτ P j,σ, (τ ) = 1. The function G j,σ, (p) is G j,σ, (p) = T 0 dt 1 r j,σ, (t 1 + p) where the function Y j,σ, is defined as follows: Y j+1,σ, (s) = r j,σ, (s) 1 − exp(r j,σ, T ) where Y j+1,σ, (s) = Y 1,σ, (s) if j = 2. The previous equation provides a link between the escape rates from the two states and can be solved recursively starting, e.g., from the initial ansatz Y 1,σ, (s) = const. A closed expression for Y 1,σ, and Y 2,σ, can be found by assuming small and taking a lin-ear approximation. If = 0, one recovers the result shown in Eq. (9), where G j,σ, (s)| =0 = r j,σ , Y j+1,σ, (s)| =0 = 1/T , and N j,σ, | =0 = 1, j = 1, 2. Using Eqs. (A1)-(A4), one can define the residence-time distribution in terms of the drift and volatility field of the system, according to the procedure described in Sec. II, and exclusively through the definition of the quasipotentials and ; the transversal components of the fields do not matter. This is a key advantage of the formulation of the problem proposed in this paper.
Building on Ref. [74], one understands that, depending on the intensity of the perturbation quasipotential and of its asymmetry in the two states 1 and 2, and on the asymmetry of the unperturbed quasipotential between the two states, the intensity of the periodic modulation contained in the function G j,σ, (p) can vary enormously. In particular, it is reasonable to guess that the properties of the periodic modulation of G j,σ, (p) are strongly influenced by the presence of a difference in the perturbation of the quasipotential 1 and 2 at all times (the classical setting corresponds to 1 = − 2 ). One then expects that the presence of characteristic peaks at times p = (2n + 1)T /2 mentioned above depends critically on the quantity | 1 − 2 |/σ 2 . Nonetheless, the investigation of these properties and the analysis of the conditions conducive to SR from this point of view are beyond the scope of this paper.