Noise-induced switching from a symmetry-protected shallow metastable state.

We consider escape from a metastable state of a nonlinear oscillator driven close to triple its eigenfrequency. The oscillator can have three stable states of period-3 vibrations and a zero-amplitude state. Because of the symmetry of period-tripling, the zero-amplitude state remains stable as the driving increases. However, it becomes shallow in the sense that the rate of escape from this state exponentially increases, while the system still lacks detailed balance. We find the escape rate and show how it scales with the parameters of the oscillator and the driving. The results facilitate using nanomechanical, Josephson-junction based, and other mesoscopic vibrational systems for studying, in a well-controlled setting, the rates of rare events in systems lacking detailed balance. They also describe how fluctuations spontaneously break the time-translation symmetry of a driven oscillator.

Introduction

Fluctuation-induced switching from a metastable state underlies a broad range of phenomena and has been attracting much interest in diverse areas, from statistical physics to chemical kinetics, biophysics, and population dynamics. Where fluctuations are weak on average, switching is a rare event, with the rate much smaller than the relaxation rate of the system. For classical and quantum systems in thermal equilibrium, switching has been well understood[1-3]. Here the major mechanisms of switching are thermal activation over the free energy barrier or, for low temperatures, tunneling. The corresponding theory has been standardly used to characterize Josephson junctions[4-6], to study magnetic systems[7-11], and for other applications.

Much progress has been made over the last few decades on the theory of switching in systems far from thermal equilibrium, see refs. [12-18] for a review. However, many problems remain open on the theory side, and much remains to be learned on the experimental side.

The experiments require well characterized nonequilibrium systems that remain stable for a time much longer than the relaxation time. A class of systems that stand out in this respect are resonantly driven mesoscopic vibrational systems where switching occurs between the states of forced vibrations. The examples range from electrons in a Penning trap to cold atoms to nano- and micromechanical systems to Josephson junction based systems, cf.[19-33].

A most detailed comparison of the theory and the experiment can be done for “shallow” metastable states. These are states with a comparatively low barrier for escape. A simple example is a state at the bottom of a shallow potential well. Even a comparatively weak noise can lead to escape from a shallow state with an appreciable rate. This significantly simplifies the experiment. Typically, a stable state becomes “shallow” when one of the parameters of the system approaches the value where the state disappears (a bifurcation point). The rate of switching from a shallow state displays a characteristic scaling with the distance to the bifurcation point in the parameter space. For stable states of forced vibrations such scaling has been found in the classical and quantum regimes[34-36] and has been observed both for the vibrations at the frequency of the driving resonant field[21,23,26,31,32] and for parametrically excited vibrations at half the drive frequency[22,24].

A major feature that underlies the scaling is that, even though the stable states are vibrational, the problem can be mapped on fluctuations of an overdamped particle in a one-dimensional potential well. This is a consequence of the onset of a “soft mode” that controls the motion near the relevant bifurcation points[37]. Moreover, because there is only one slow variable, the corresponding slow motion has detailed balance.

In this paper we consider escape from a shallow metastable state with no detailed balance. We show that such a state exists even in a simple system which has only two dynamical variables. The dynamics is not controlled by soft modes. Rather the emergence of the shallow state is a consequence of the symmetry of the system. The considered model is minimalistic: the scaled equations of motion without noise contain a single parameter.

The physical system we consider is a vibrational mode driven close to triple eigenfrequency. The classical dynamics of such a mode in the absence of noise has been well understood[38]. For sufficiently strong driving, the mode can have three stable period-3 vibrational states, which all have the same amplitude and differ in phase by 2π/3. However, the state with no vibrations (except maybe vibrations at the drive frequency with a very small amplitude) is also stable. It remains stable as the driving amplitude increases. We call it a zero-amplitude state.

The basin of attraction of the zero-amplitude state, i.e., the range in which the mode approaches this state from an initially prepared state, can be seen from Fig. 1(a,b). It is the interior of the region centered at the origin and limited by the solid lines (the separatrices). As we show, the size of this region decreases with the increasing driving. However, it remains nonzero for a finite driving amplitude. Indeed, annihilation of the zero-amplitude state would require that it merges simultaneously with the three saddle points, which correspond to unstable period-3 vibrations and are shown by squares in Fig. 1. The persistence of the zero-amplitude state is therefore a consequence of the symmetry, which is the symmetry of the period-3 vibrational states with respect to incrementing their phase by 2π/3 or, in other words, the symmetry of Fig. 1 with respect to rotation by 2π/3.

Figure 1

(a) The phase portrait of the mode that displays period tripling; Q and P are the quadratures (the coordinate and momentum in the rotating frame). Circles and squares show the stable states and the saddle points, respectively. The lines show the separatrices. The plot refers to κ = 0.4 and f = 2 in Eq. (4). (b) The part of the phase portrait inside the dashed square in (a) in the scaled coordinate and momentum q1 = fQ and q2 = fP. For the dynamics is described by Eq. (7). The green line that comes from the stable state shows the most probable path followed in escape for κ = 0.4. (c) The effective Hamiltonian (8) of motion around the shallow state q1 = q2 = 0 in the absence of dissipation. With dissipation, the local maximum at the origin becomes a stable state. (d) The contour plot of the Hamiltonian (c). The squares show the saddle points. These points shift in the presence of dissipation, as seen in (b), but the stable state remains at q1 = q2 = 0.

As the basin of attraction of the zero-amplitude state shrinks with the increasing driving, this state becomes more and more shallow. Thus, for strong driving, the zero-amplitude state is a symmetry-protected shallow metastable state. We are not aware of a simpler shallow metastable state with the dynamics that is not controlled by a soft mode no matter how small the basin of attraction becomes.

Resonant period tripling leads to unusual quantum dynamics of a vibrational mode[39-42]. In the quantum regime, a period tripling was recently observed in a superconducting resonator with several nonlinearly coupled modes[43,44]. There is a significant difference between quantum fluctuations about the shallow zero-amplitude state and the shallow states near bifurcation points. Quantum fluctuations near a bifurcation point are very similar to classical fluctuations, since the dynamics is controlled by a single slow variable[36]. In contrast, the dynamics near the considered shallow zero-amplitude state is described by two non-commuting dynamical variables, and therefore quantum fluctuations qualitatively differ from classical fluctuations. In the present paper we study escape due to classical fluctuations.

The model

Mesoscopic vibrational systems, from nano- and micro-mechanical modes to Josephson-junction based systems, are nonlinear. A major effect of nonlinearity is the dependence of the vibration frequency on the amplitude. To the leading order in the nonlinearity, this dependence is well described by the Duffing model, in which one takes into account the quartic term in the expansion of the potential energy in the mode coordinate q. To resonantly excite period-3 vibrations one has to apply a force at frequency ω close to triple the eigenfrequency ω0. The simplest Hamiltonian that describes the mode dynamics isHere, p is the mode momentum, γ is the nonlinearity parameter, F is the driving amplitude, and we have set the mass of the mode m = 1. The condition of the driving being resonant means that the frequency detuning δω is small,While assuming the detuning small, we will further assume that it is nonzero and that both δω and γ are positive. These conditions are not necessary, the results immediately extend to the cases where δω, γ < 0 and also to δω = 0; the results extend also to the case where the driving is additive and the term in Eq. (1) that describes the driving has the form , cf.[41].

Along with the detuning, we will assume that the nonlinearity and the driving are also small, so that the vibrations are nearly sinusoidal and, equivalently, the change of the mode frequency due to the nonlinearity is small for the typical displacement , i.e., .

The rotating wave approximation

Resonant dynamics of the mode is conveniently described by the amplitude and phase of the vibrations at frequency ω/3, which vary on the time scale that largely exceeds the vibration period 2π/ω. To describe this dynamics we switch to the frame that rotates at frequency ω/3 and introduce the scaled coordinate Q and momentum P in this frame. The corresponding transformation is[note that these are two equations for two real variables Q(t) and P(t)].

The transformed Hamiltonian H contains time-independent terms and terms that oscillate at the frequency 2ω/3 and its overtones. The effect of these terms is small in the considered parameter range and can be disregarded, i.e., the dynamics can be analyzed in the rotating wave approximation (RWA). In the RWA, the Hamiltonian (1) in the variables Q, P becomesThe Hamiltonian g contains a single parameter, the scaled amplitude of the driving field f. In the considered range of small nonlinearity and small detuning δω, parameter f can be large even for weak driving.

The coupling of the mode to a thermal reservoir leads to decay of the vibrations and to noise. For several microscopic mechanisms of the coupling, on the time scale that largely exceeds and the correlation time of the relevant fluctuations of the reservoir, the dynamics of the weakly nonlinear mode is Markovian (Methods). In the simplest case of the coupling linear in the mode coordinate, the dynamical equations of motion in the dimensionless time readHere κ is the dimensionless decay rate of the mode. It is simply related to the friction coefficient Γ used when the mode dissipation is phenomenologically described by a friction force −2Γdq/dt and the friction is weak, ; in this case, . However, Eq. (4) can apply even where in the laboratory frame the dynamics is non-Markovian (Methods).

The terms and in Eq. (4) are independent δ-correlated Gaussian noises,For thermal noise .

The phase portrait of the mode in the rotating frame in the absence of noise is shown in Fig. 1(a). For , Eq. (4) has three stable stationary solutions, which correspond to three stable period-3 vibrational states with nonzero scaled squared amplitude in the laboratory frame, three saddle points, which correspond to unstable vibrational states with scaled squared amplitude , and a stable state with Q = P = 0,The stable states and the saddle points with the squared amplitudes (6) are located at the vertices of equilateral triangles, cf. Fig. 1(b). In a standard fashion, the separatrices go through the saddle points and divide the phase plane (Q, P) into the basins of attraction of different stable states. Importantly, the basins of attraction of the states with nonzero vibrational amplitude have a common boundary only with the zero-amplitude state, but not with each other.

With the increasing amplitude of the driving field f, the stable states with nonzero amplitude (6) move away from the origin. In contrast, the saddle points move toward the origin, for large f. Therefore, the basin of attraction of the zero-amplitude state shrinks, but at the same time, this state remains stable for any finite f.

Results

Dynamics near the zero-amplitude state for large driving amplitudes

For strong driving, the dynamics of the system in the vicinity of the zero-amplitude state can be conveniently studied by switching to the variables q1 = fQ and q2 = fP. They can be thought of as components of a vector q ≡ (q1, q2). For , from Eq. (4) the equations of motion for q1,2 in the range have the formHere is the Levi-Civita tensor, and ; is the gradient,and ξ1,2(τ) are two independent white Gaussian noises,

Function is an analog of a Hamiltonian, if one thinks of q1 as a coordinate and q2 as a momentum. It determines the conservative motion of the strongly modulated mode in the absence of dissipation and noise. It has no parameters. Its structure is therefore universal. It is shown in Fig. 1(c,d). As seen in this figure, has a local maximum at q = 0 and three saddle points. The local maximum corresponds to the zero-amplitude state, in the presence of dissipation.

The phase portrait of the system (7) in the presence of dissipation but with no noise is shown in Fig. 1(b). Dissipation shifts the saddle points of , but they remain saddle points. We use the same color coding for the saddle points of in Fig. 1(d) and the saddle points of the full dynamics (7) in Fig. 1(a). The separatrices in Fig. 1(b) are the boundaries of the basin of attraction to the zero-amplitude state.

Overall, in the absence of noise the dynamics described by Eq. (7) depends on only one parameter, the scaled decay rate κ. Varying κ leads only to a quantitative change of the phase portrait, the structure remains intact, as seen from Fig. 1(b). The simple topology of the phase portrait indicates the universality of the dynamics. Quite remarkably, this universality is not related to the onset of a soft mode. There is no slowing down.

Another important feature of the dynamics (7) is that the effective noise intensity increases with the increasing driving amplitude F ∝ f. Such increase is essentially a consequence of the rescaling of the dynamical variables. As the basin of attraction of the zero-amplitude state on the original (Q, P) phase plane shrinks with the increasing f, the noise becomes effectively stronger. This is what makes the zero-amplitude state “shallow”.

By linearizing equations of motion about q = 0 we see that, for weak noise, the mean-square displacement about the zero-amplitude state isThe linearization makes sense only if , whereis the squared distance of the saddle points from the origin, with the account taken of the dissipation. The positions of the three saddle points can be written as with . The latter equation gives three values of ϕ that correspond to different saddle points and differ by 2π/3.

With the increasing , the mean-square displacement increases. Once it approaches , fluctuations about the stable state may no longer be assumed small. Physically, it means that the system placed initially near the zero-amplitude state quickly escapes from the basin of attraction of this state and switches to one of the states of period-3 vibrations. For smaller f the escape rate is smaller, but still it may be not exceedingly small even for a weak noise, justifying the term “shallow metastable state” as applied to the zero-amplitude state.

Escape rate from the shallow state with no detailed balance

The theory of escape from a metastable state of a white-noise driven system is well-established[13]. The underlying idea is that, when the noise is weak on average, escape occurs as a result of a rare fluctuation. In this fluctuation, a large outburst of noise drives the system over the boundary of the basin of attraction of the initially occupied stable state. Once the boundary is crossed, noise is no longer needed, the system moves to another state “on its own”. An outburst of noise is a certain realization of a random force. The outbursts needed for switching are exponentially unlikely for a weak Gaussian noise. In addition, the probabilities of different appropriate outbursts are exponentially different. The rate of escape is determined by the most probable of them, i.e., by the least improbable appropriate evolution of the random force in time. Through the equations of motion, such force leads to the corresponding trajectory of the system[45]. This trajectory is often called[6,46] the most probable escape path (MPEP).

We will consider escape from the zero-amplitude state using Eq. (7). The conventional analysis refers to systems with no symmetry. In contrast, our system has a three-fold symmetry. The attraction basin of the zero-amplitude state is bound by three separatrices, which can be obtained from each other by a rotation by 2π/3 on the (q1, q2) plane, see Fig. 1(b). The probability to cross any of them in escape is the same, and therefore the total escape rate is three times the escape rate for crossing one of them.

For weak noise, it is most probable to cross a separatrix near the saddle point[34,47,48]. For the system (7), the rate of escape over one of the separatrices W0 for small noise intensity has the formwhere C is a constant that smoothly (nonexponentially) depends on the parameters. Equation (9) reminds the Kramers formula[1] for the rate of activated escape from a potential well. However, the dynamics (7) does not correspond to Brownian motion in a potential well, and R is not a height of a potential barrier.

The effective activation energy R for the considered fluctuating dissipative system is given by the action of an auxiliary conservative system[13],Here is the Lagrangian of an auxiliary system; this is a Hamiltonian system, with no relaxation and no noise. Its appropriate Hamiltonian trajectory qopt(τ) gives the MPEP, i.e., the trajectory which the initial dissipative noisy system is most likely to follow in escape. This trajectory starts at q = 0 for τ → −∞ and goes to one of the saddle points for τ → ∞ (Methods).

The effective activation energy in the limiting cases

Equation (10) allows one to calculate the effective activation energy R for any κ by numerically solving the variational equations of motion, or the corresponding Hamiltonian equations, with the appropriate boundary conditions (Methods). Analytically, explicit values of R can be calculated in the limiting cases of large and small κ. We start with the case . For κ = 0 the noise-free trajectories of the system (7) are closed loops with , as seen in Fig. 1(d). For small κ the noise-free trajectories become tight spirals that spiral toward q = 0, with slowly increasing toward . The MPEP qopt(τ) is also a tight spiral, but it spirals from q = 0 toward the saddle points. To the lowest order in κ, the value of R is determined by the action (10) accumulated on this trajectory until reaches its value at the saddle point[34,49,50].

Importantly, the multiplicity of the saddle points makes no difference, to the leading order in κ. Indeed, for κ → 0 all saddle points have the same , cf. Fig. 1(c,d). Therefore the variational equations for the optimal path can be solved in the same way as in systems with a single saddle point. The resulting general expression for R is[34,49,50]where is the interior of the orbit with a given . Such orbits are shown in Fig. 1(d). From Eq. (7), the orbit is a circle for small and becomes a triangle for , with the two sides given by and the 3rd side given by q1 = 1/2.

From Eq. (7), , and therefore

This expression shows that the effective activation energy R linearly increases with the scaled decay rate κ where κ is small.

We now consider the case of fast decay, . Still we assume that , so that the dynamics in the considered part of the phase plane is described by Eq. (7). In this case it is convenient to change variables in Eqs. (7) and (10) by setting . Then , whereas, to the leading order in κ, the components of the vector k become with , and . Therefore we can write the expression for the activation energy as

The parameter κ has been scaled out of Eq. (13). The Lagrangian does not contain any parameters. Therefore R′ is a number. This number can be found by solving the variational problem (13) numerically for the extreme trajectories that start at q′ = 0 at τ′ → −∞ and for τ → ∞ go to one of the saddle points, which are located at (cf. Methods). The numerical solution gives . Therefore,

The asymptotic results (12) and (14) explicitly describe the dependence of the effective activation energy of escape on the single parameter of the nonlinear dynamics, κ, in the limiting cases and .

Numerical results

We have studied escape from the zero-amplitude state numerically using three approaches. First, we performed numerical simulations of the full stochastic equations of motion (4) for several values of the scaled driving force f. The results for f = 5 are shown in Fig. 2 by crosses. Then we performed simulations of the scaled stochastic equations (7) that refer to large f and do not contain f other than in the noise intensity . The results are shown in Fig. 2 by full circles. Then we solved the noise-free variational problem (10). The numerical integration of the stochastic equations was done following the standard routine, cf.[51]. To find the escape rate W0, for each parameter value we assembled 3000 trajectories that went from the vicinity of the zero-amplitude state to the area well behind the basin of attraction of this state. We made sure that the result was independent of the chosen boundary of this area.

Figure 2

The logarithm of the rate of escape from the zero-amplitude state log W0 as a function of the scaled noise intensity as obtained by numerically simulating the Langevin dynamics. The full circles show the results obtained for the dynamics described by Eq. (7), where the driving field amplitude f was scaled out. The crosses show the results obtained by simulating the full Eq. (4) for f = 5; respectively, the noise intensity in these simulations was . The results show excellent scaling in the both cases. The straight dashed lines are guide for the eye. The data with brown, yellow, cyan, black, magenta, green, red, and blue full circles refer to κ = 4.5, 3.5, 3, 2.5, 2, 1.5, 0.5, 0.25.

The results of the numerical simulations shown in Fig. 2 clearly demonstrate the activation dependence of the escape rate on the noise intensity for all values of the scaled decay rate we explored: logW0 is linear in . Already for f = 5, the results obtained from the full Eq. (4), which explicitly contain f and refer to the noise intensity , were extremely close to those obtained from the scaled Eq. (7) that refer to the limit.

In Fig. 3 we compare the results of the simulations with the results obtained by solving the variational problems (10) and (13) (Methods). The effective activation energy of escape increases with the increasing decay rate κ. The results show that the effective activation energy R is well-described by the asymptotic small-κ expression (12) in the range . In the range the asymptotic large-κ expression (14) works reasonably well.

Figure 3

The dependence of the effective activation energy R of escape from the shallow metastable zero-amplitude state on the scaled decay rate of the mode κ. The results refer to strong driving, The solid lines show the solution of the full variational problem (10). The red and green dashed lines show the asymptotic results (12) and (13) for small and large κ, respectively. The full circles show the values of R obtained from numerical simulations of the stochastic equation of motion (7).

It is also seen from Fig. 3 that the values of the activation energy obtained by simulations and from the analytical theory are in excellent agreement. This agreement holds for all values of κ we explored.

Discussion

The results of this paper bear on two problems. One is the possibility to facilitate the comparison of the theory and the experiment on escape from a metastable state of a generic thermally nonequilibrium system in the presence of weak noise. The second is the spontaneous breaking of the time-translation symmetry in vibrational systems where the symmetry-preserving state is dynamically stable. We have shown that both problems can be addressed by studying period tripling in mesoscopic vibrational systems, including nano-mechanical and Josephson-junction based systems in particular. Where such systems are driven close to triple the eigenfrequency, they can have three stable period-3 states and also a stable “zero-amplitude” state where they do not oscillate or oscillate with a small amplitude at the drive frequency.

The stability of the zero-amplitude state is symmetry-protected, and because of the nature of the symmetry, the dynamics in the vicinity of this state has no detailed balance. In this sense, this dynamics is generic for a nonequilibrium system. At the same time, escape from the zero-amplitude state into one of the period-3 states, and thus the breaking of the symmetry of the driving with respect to time translation by 2π/ω, can have a comparatively large probability for strong driving even where the noise is weak. This means that the zero-amplitude state becomes shallow while remaining dynamically stable. We have shown that the escape rate W0 displays universal scaling with the parameters of the vibrational mode as well as the amplitude F and the frequency ω of the drive. For thermal noise Eq. (9) takes the formwhere κ is the scaled decay rate. The dependence of the effective activation energy of escape R on the single parameter κ has been found analytically in the limiting cases and numerically in the general case. Quite remarkably, the simple explicit results in these limiting cases well describe the escape rate almost in the entire range of the parameters of the system.

The possibility to exponentially strongly increase the escape rate by varying the parameters of the driving is important for studying escape in the experiment. The explicit scaling of the escape rate with the parameters of the drive opens a way of a detailed quantitative comparison of the theory of escape in systems lacking detailed balance with the experiment.

Methods

Markovian dynamics of a driven oscillator in slow time

Coupling of a dynamical system to a thermal reservoir leads to relaxation and fluctuations in the system. Generally, the ensuing dynamics is non-Markovian. Phenomenologically, it is often described by adding to the Hamiltonian equations of motion a retarded friction force and a Gaussian noise with the corresponding correlation function[52,53]. A simple microscopic model that underlies such description is the coupling to a thermal bath that consists of harmonic oscillators, with the coupling being bilinear in the coordinate of the considered dynamical system and the coordinates of the oscillators[54-56]. If the system can be thought of as a particle in a potential well, and if the correlation time of the thermal reservoir is short and the decay time of the system is long compared to the period of intrawell vibrations, the long-time dynamics can be mapped on an effective Markov process with energy-dependent decay and diffusion coefficients, see[57,58] and references therein.

The dynamics of the system considered in the present paper, i.e., a weakly nonlinear oscillator weakly coupled to a thermal reservoir, can be described differently and for a fairly general form of the coupling. This is because the oscillator is a “narrow-band filter”, that is, it is susceptible to fluctuations in the thermal reservoir in a narrow frequency range centered at the oscillator eigenfrequency. For a harmonic oscillator bilinearly coupled to a bath of harmonic oscillators, the rigorous analysis by Bogolyubov[54] showed that, in slow time compared to the vibration period 2π/ω0, the oscillator dynamics is Markovian provided the density of states of the bath weighted with the interaction is smooth near the oscillator eigenfrequency.

For the considered oscillator, the nonlinear part of the vibration energy is small compared to the harmonic part,Because of the nonlinearity, the vibration frequency depends on the vibration amplitude, but for small nonlinearity it remains close to ω0. The frequency ω/3 is also close to ω0. Therefore the physical argument that the oscillator is sensitive to fluctuations in the thermal reservoir at frequency ≈ω0 applies to the driven nonlinear oscillator. The oscillator dynamics in the slow time can be analyzed[34,59] using the asymptotic technique, which is to some extent similar to the conventional ladder approximation in the Green’s function theory of electronic systems[60]. This technique applies in the quantum and classical limits, and in the classical limit leads to Eq. (4) if the interaction with the thermal reservoir is of the form H = qhb, where hb depends only on the dynamical variables of the thermal reservoir. In particular, if the reservoir is a set of harmonic oscillators, hb can have linear and nonlinear terms in the coordinates of these oscillators. For the considered oscillator with mass m the decay rate is[the eigenfrequency ω0 is also renormalized by a term ]. The approximation applies for , where Δω includes the frequency detuning ω0 − ω/3 and the nonlinear frequency shift .

In addition to being coupled to a thermal reservoir, mesoscopic vibrational systems are often subjected to an external noise[29], which in many cases is Gaussian. An oscillator is particularly perceptive to the noise at frequencies close to ω0. If the power spectrum of the noise is smooth around ω0, in the slow time the noise acts on the oscillator as a white noise. To allow for the presence of this noise, the noise intensity D in Eq. (4) is chosen as a parameter, which can be larger than the thermal noise intensity.

Calculating the most probable escape path

A simple numerical procedure of finding the MPEP is based on solving the Hamiltonian equations of motion for the auxiliary system with the Lagrangian given by Eq. (10). The Hamiltonian of this system and the corresponding equations areThe MPEP corresponds to the trajectory that starts for τ → −∞ at the stable state of the original system and arrives for τ → ∞ to one of the saddle points (a real time instanton), cf.[13,34,46,50,61]. On this trajectory . The activation energy of escape R is given by the dynamical action accumulated by moving along this trajectory,with S(τ → −∞) → 0.

The Hamiltonian Eq. (17) can be solved by first linearizing them near the stationary point q = p = 0; this gives , , and explicitly determines the matrix ; we find . One can then use the shooting method by choosing a small q(0), with p(0), S(0) given by the above expression, and integrating Eqs. (17) and (18) so that the trajectory approaches one of the three saddle points. An example of the MPEP obtained this way is shown in Fig. 1(b). The results were used to obtain the solid line in Fig. 3.

The green dashed line in Fig. 3 was obtained by the same procedure, with the Hamiltonian of the form of , where the vector .

We note that our analysis refers to a system with weak damping, which becomes multistable only due to resonant periodic driving. The analysis is done in the rotating wave approximation, and therefore the Hamiltonian for the escape trajectories (17) is independent of time. This is qualitatively different from the problem of escape of a Brownian particle from a periodically modulated metastable potential well, which has attracted much attention in the literature[62,63] and where, for an overdamped particle, both the exponent and the prefactor in the escape rate have been found[64-67]. The Hamiltonian for escape trajectories in this case explicitly depends on time. Recent numerical results on escape from a periodically modulated potential well can be found in ref. [68].