Quantum fluctuations in modulated nonlinear oscillators

With a modulated oscillator, we study several effects of quantum fluctuations far from thermal equilibrium. One of them is quantum heating, where quantum fluctuations lead to a finite-width distribution of a resonantly modulated oscillator over its quasienergy (Floquet) states. We also analyze large rare fluctuations responsible for the tail of the quasienergy distribution and switching between the states of forced vibrations. We find an observable characteristic of these fluctuations, the most probable paths followed by the quasienergy in rare events, and in particular in switching. We also explore the discontinuous change of the most probable switching path where the detailed balance condition is broken. For oscillators modulated by a nonresonant field, we compare different mechanisms of the field-induced cooling and heating of the oscillator.


Introduction
The last few years have seen an upsurge in the interest in the dynamics of modulated nonlinear oscillators [1]. There have emerged several new areas of research where this dynamics plays a central role, such as nanomechanics, cavity optomechanics, and circuit quantum electrodynamics. The vibrational systems of the new generation are mesoscopic. On the one hand, they can be individually accessed, similar to macroscopic systems, and are well-characterized. On the other hand, since they are small, they experience comparatively strong fluctuations of thermal and quantum origin. This makes their dynamics interesting on its own and also enables using modulated oscillators to address a number of fundamental problems of physics far from thermal equilibrium.
Many nontrivial aspects of the oscillator dynamics are related to the nonlinearity. Essentially all currently studied mesoscopic vibrational systems display nonlinearity. For weak damping, even small nonlinearity becomes important. It makes the frequencies of transitions between adjacent oscillator energy levels different. Where several levels are occupied, the dynamics strongly depends on the interrelation between the width of the ensued frequency comb and the oscillator decay rate. An important consequence of the nonlinearity is that, when an oscillator is resonantly modulated, it can have coexisting states of forced vibrations, i.e., display bistability [2].
One of the general physics problems addressed with modulated nonlinear oscillators is fluctuation-induced switching in systems that lack detailed balance, see [3,4,5,6,7,8,9,10,11,12,13,14] for the classical and [15,16,17,18,19,20,21,22,23,24,25,26] for the quantum regime. A remarkable property of the switching rate in the quantum regime is fragility. The rate W sw calculated for T = 0, where the system has detailed balance [27], is exponentially different from the rate calculated for T > 0, where the detailed balance is broken [15,18]. Recently the effect of fragility of the rates of rare events was also found in the problem of population dynamics [28]. There, too, a small change of the control parameter (infinitesimal, in the semiclassical limit) leads to an exponentially strong rate change. The nature of the dynamics and the sources of fluctuations in a quantum oscillator and in population dynamics are totally different, and it is important to understand how it happens that they display common singular features.
An important source of quantum fluctuations is the coupling of the oscillator to a thermal bath. It leads to oscillator relaxation via emission of excitations in the bath accompanied by transitions between the oscillator energy levels. The transitions lead to relaxation only on average, in fact they happen at random, giving rise to a peculiar quantum noise. For a resonantly modulated oscillator, the noise causes diffusion over the oscillator quantum states in the external field, which are the quasienergy (Floquet) states. As a result, even where the bath temperature is T = 0, the distribution over the states has a finite width, the effect of quantum heating [29].
We discuss quantum heating for a resonantly modulated oscillator and compare the predictions with the recent experiment [30] where the effect was observed. The spectral manifestation of quantum heating is considered, with the focus on the influence of dissipation on the oscillator spectral characteristic of interest for sideband spectroscopy, the technique which was nicely implemented in the experiment [30] using a microwave cavity with an embedded qubit.
We also study switching between the stable states of forced vibrations of an oscillator modulated close to its eigenfrequency. As quantum heating, switching occurs because of the quantum-noise induced diffusion over the oscillator states. It reminds switching of a classical Brownian particle over the potential barrier due to diffusion over energy [31] and therefore is called quantum activation. Generally, the rate of quantum activation largely exceeds the rate of switching via quantum tunneling. We develop an approach to calculating the rate of quantum activation, which naturally connects to the conventional formulation of the rare events theory in chemical and biological reaction systems and in population dynamics [32,33]. This approach provides a new insight into the fragility of the switching rate of the oscillator.
The dynamics of a periodically modulated harmonic quantum oscillator coupled to a thermal reservoir is one of exactly solvable problems of physical kinetics [34,35,36].
However, the solution disregards the fact that nonresonant modulation can open a new channel for oscillator relaxation, where a transition between the oscillator energy levels is accompanied by an emission (absorption) of an excitation in the medium, while the energy deficit is compensated by the modulation. Alternatively, the role of the medium can be played by another mode with the relaxation rate higher than that of the oscillator [37]. This mechanism underlies the cooling studied in cavity optomechanics. We provide a brief comment in order to unify various mechanisms of the change of the quantum distribution of the oscillator by nonresonant modulation.

Hamiltonian in the rotating frame
A major type of the internal oscillator nonlinearity of interest for the effects we will discuss is the Duffing nonlinearity, where the potential energy has a term quartic in the oscillator displacement q; in quantum optics, it corresponds to the Kerr nonlinearity. The simplest types of resonant modulation that lead to the bistability of the oscillator are additive modulation at frequency ω F close to the oscillator eigenfrequency ω 0 and parametric modulation (modulation of the oscillator frequency) at frequency ≈ 2ω 0 [2]. The analysis of quantum fluctuations in these two systems has much in common [38], and the method that we will develop here applies to the both types of systems. For concreteness, we will consider here additive modulation. The oscillator Hamiltonian is where q and p are the oscillator coordinate and momentum, the mass is set equal to one, γ is the anharmonicity parameter, and A is the modulation amplitude. We assume that the the modulation is resonant and not too strong, so that A periodically modulated oscillator is described by the Floquet, or quasienergy states Ψ ε (t). They provide a solution of the Schrödinger equation where t F = 2π/ω F . This expression defines quasienergy ε. To find quasienergies and to describe the oscillator dynamics it is convenient to change to the rotating frame. This is done by the standard canonical transformation U(t) = exp −ia † a ω F t , where a † and a are the raising and lowering operators of the oscillator. We introduce slowly varying dimensionless coordinate Q and momentum P in the rotating frame, U † (t)qU(t) = C(Q cos ϕ + P sin ϕ), U † (t)pU(t) = −Cω F (Q sin ϕ − P cos ϕ).
Here, ϕ = ω F t and the scaling factor is C = |8ω F (ω F − ω 0 )/3γ| 1/2 . The commutation relation between P and Q has the form Parameter λ ∝ plays the role of the Planck constant in the quantum dynamics in the rotating frame. It is determined by the oscillator nonlinearity, λ ∝ γ. For concreteness we assume that γ, δω > 0; the oscillator displays bistability for γδω > 0. In the range (2) the oscillator dynamics can be studied in the rotating wave approximation (RWA). The RWA Hamiltonian is where β is the scaled modulation intensity and E sl = γC 4 is the characteristic energy of motion in the rotating frame. This motion is slow on the time scale ω −1 F . Note that E sl ≪ ω 2 0 q 2 ∼ ω 2 0 C 2 . Operatorĝ = g(Q, P ) is the dimensionless Hamiltonian in the rotating frame. In the RWA, the Schrödinger equation for the RWA wave function ψ(Q) in dimensionless slow time τ reads Operatorĝ is independent of time and has a discrete spectrum,ĝ|n = g n |n . The eigenvalues g n give the quasienergies in the RWA, ε n = (3E sl /8)g n (we are using an extended ε-axis rather than limiting ε to the analog of the first Brillouin zone 0 ≤ ε < ω F ). Function g(Q, P ) has the shape of a tilted Mexican hat and is shown in Fig. 1(a); the quasienergy levels are shown in Fig. 1(b) . In contrast to the Hamiltonian H 0 ,ĝ is not a sum of the kinetic and potential energies. As seen from Fig. 1, the eigenstates localized near the local maximum of g(Q, P ) correspond to semiclassical orbits on the surface of the "inner dome" of g(Q, P ); these states become stronger localized as g n increases toward the local maximum of g(Q, P ). The quasienergy level spacing ∝ λE sl is small compared to the distance between the oscillator energy levels in the absence of modulation, |ε n − ε n+1 | ∼ λE sl ≪ ω 0 .

Master equation for linear coupling to the bath
The analysis of the oscillator dynamics is often done assuming that the oscillator is coupled to a thermal bath in such a way that the coupling energy is linear in the oscillator coordinate q and thus in the oscillator ladder operators a, a † [34]. In this case the coupling energy H i and the typical relaxation rate Γ are of the form where h b depends on the bath variables only and . . . b denotes thermal averaging over the bath states. Relaxation (6) corresponds to transitions between neighboring energy levels of the oscillator in the lab frame, with energy transferred to bath excitations, see Fig. 1(c). The renormalization of the oscillator parameters due to the coupling is assumed to have been incorporated. For a smooth density of states of the bath, resonant modulation does not change the decay rate parameter, Γ(ω F ) ≈ Γ(ω 0 ). However, it excites the oscillator, as sketched in Fig. 1(c). In a stationary state of forced vibrations (in the lab frame) the energy provided by the modulation is balanced by the relaxation.
To the second order in the interaction (6), the master equation for the oscillator density matrix ρ in dimensionless time τ readṡ Here, the term ∝ [ρ,ĝ] describes dissipation-free motion, cf. (5). Operatorκρ describes dissipation and has the same form as in the absence of oscillator modulation [39,40]; κ is the dimensionless decay rate, andn is the oscillator Planck number, In the classical limit λ → 0 the oscillator described by (7) can have one or two stable states of forced vibrations. Their positions in the rotating frame (Q a , P a ) are given by the stable stationary solutions of the classical equations of motion of the oscillatoṙ Equations (9) are, essentially, the mean-field equations for the moments Tr(Qρ), Tr(P ρ) for λ → 0. For small damping Q a and P a are close to the extrema of g(Q, P ).

Balance equation
We will concentrate on the oscillator dynamics in the case where the oscillator is strongly underdamped and its motion is semiclassical, In this case the number of quasienergy states localized about the extrema of g(Q, P ) is large, ∝ 1/λ [the scaled quasienergies of such states g n lie between the value of g at the corresponding extremum and the saddle point value g S of g(Q, P ) in Fig. 1]. Also, the spacing between the levels is large compared to their width, |g n − g n±1 | ≫ λκ. Where the latter condition is met, the off-diagonal matrix elements of the density matrix on the quasienergy states ρ nm ≡ n|ρ|m (n = m) are small. To the lowest order in κ the oscillator dynamics can be described by the balance equation for the populations ρ nn of quasienergy states. From (7) where a nm ≡ n|a|m . We disregard tunneling when defining functions |n ≡ ψ n (Q), i.e., we use the wave functions localized about the extrema of g(Q, P ); the effect of tunneling is exponentially small for λ ≪ 1. We count the localized states off from the corresponding extremum, i.e., for a given extremum the state with n = 0 has g n closest to g(Q, P ) at the extremum. An important feature of the rates of interstate transitions W mn is that, even for T = 0 (and thusn = 0), there are transitions both toward and away from the extrema of g(Q, P ). This is because the wave functions |n are linear combinations of the wave functions of the oscillator Fock states, see Fig. 1(c). Therefore, even though relaxation corresponds to transitions down in the oscillator energy in Fig. 1(c), the transitions up and down the quasienergy have nonzero rates. One can show that, for the both extrema of g(Q, P ), the rates of transitions toward an extremum are larger than away from it. Therefore, depending on where the system was prepared initially, it would most likely move to one or the other extremum of g(Q, P ). This is why the extrema correspond to the stable states of forced vibrations of the modulated oscillator in the classical limit of a large number of localized states.
For small effective Planck constant λ ≪ 1, the rates W mn can be calculated in an explicit form by finding the matrix elements a mn in the WKB approximation [15,41]. The problem is then related to that of classical conservative motion with Hamiltonian g(Q, P ) and with equations of motion of the formQ = ∂ P g,Ṗ = −∂ Q g. A significant simplification comes from the fact that the classical trajectories Q(τ ; g) are described by the Jacobi elliptic functions. As a result, Q(τ ; g) is double-periodic on the complexτ plane, with real period τ (1) p (g) and complex period τ [42]. It can be calculated along an appropriately chosen closed contour on the complex τ -plane and is determined by the pole of a(τ ; g n ).
In particular, for the states localized about the local maximum of g(Q, P ) we obtain Here, τ * ≡ τ * (g n ) and τ is the pole of Q(τ ; g) closest to the real axis; ν(g) is the dimensionless frequency of vibrations in the rotating frame with quasienergy g. To the leading order in λ, we have W n n+k = W n−k n for n, n ± k ≫ 1.

Effective temperature of vibrations about a stable state
Equation (12) has to be modified for states very close to the extrema of g(Q, P ). Near these extrema the classical motion of the oscillator in the rotating frame is harmonic vibrations. One can introduce raising and lowering operators b and b † for these vibrations (via squeezing transformation) and expand g(Q, P ) near an extremum as [(Q a , P = 0) is the position of the considered extremum; it is given by equation The derivatives of g in (13) are evaluated at (Q a , P = 0). Parameter φ * is given by equation tanh ϕ * = (|∂ 2 Q g| 1/2 − |∂ 2 P g| 1/2 )/(|∂ 2 Q g| 1/2 + |∂ 2 P g| 1/2 ). From (11) and (13), near an extremum of g we have W m+1 m = 2κ(m + 1)(n e + 1), W m m+1 = 2κ(m + 1)n e , n e =n cosh 2 φ * + (n + 1) sinh 2 φ * , whereas W m m+k = 0 for |k| > 1. Equation (14) is a familiar expression for the transition rates between the states of a harmonic oscillator coupled to a thermal bath, withn e being the Planck number of the excitations of this fictitious bath at the frequency of vibrations in the rotating frame ν 0 δω. From (14), the stationary distribution of the modulated oscillator over its quasienergy states near an extremum of g(Q, P ) is of the Boltzmann type, with effective temperature T e = λν 0 / ln[(n e + 1)/n e ] [18,38]. In agreement with the qualitative picture discussed above, this temperature is nonzero even where the temperature of the true bath is T = 0. This is the effect of quantum heating due to quantum fluctuations in a nonequilibrium system.
Quantum heating of a resonantly modulated oscillator was recently observed in an elegant experiment [30] using a mode of a microwave cavity with an embedded Josephson junction [43]. The occupation of the excited quasienergy states was revealed using a two-level system (a transmon qubit) as a probe. As seen from Fig. 2, the results of the experiment are in a qualitative agreement with the above theory. The agreement improves for larger scaled field intensity β (4), where the ratio κ/ν 0 is smaller. It is in the range of small κ/ν 0 that the quantum temperature is a good  Figure 2. (a) The effective Planck numbern e of the vibrations about the largeamplitude state of the modulated oscillator, which corresponds to the minimum of g(Q, P ) in the small-damping limit; β is the scaled driving field intensity (4). Squares: experimental data [30]. Solid line: equation (14) forn = 0 (see also [38]). Triangles: the estimate of the experimentally measured parameter discussed in Appendix A. (b) The scaled power spectra of the oscillator occupation numbern = a † a for vibrations about the large-amplitude stable state for |δω|/κ = 3.9 (the value used in [30]). The black and red curves correspond to β = 0.17 and 0.8. The triangles in (a) are determined from the ratio r Φ of the heights of the lower and higher peaks of Φ nn as r characteristic of the distribution over quasienergy states, as the lifetime of these states largely exceeds the reciprocal level spacing (scaled by ). Still, even for not too small κ/ν 0 , the technique developed in [30] makes it possible to reveal the broadening of the stationary distribution of the modulated oscillator due to quantum fluctuations far from equilibrium. A characteristic of this effect is discussed in Appendix A and the corresponding results are shown in Fig. 2.

Switching between the stable states
The effect of diffusion over quasienergy states due to quantum fluctuations is not limited to the quantum heating described above. Along with small fluctuations, which lead to comparatively small deviations of quasienergy from its value at an extremum of g(Q, P ), there occasionally occur large fluctuations. They push the oscillator far away from the initially occupied extremum. It is clear that, if as a result of such fluctuation, the oscillator goes "over the quasienergy barrier" to states localized about the other extremum, with probability ≈ 1 it will then approach this other extremum. Such transition corresponds to switching between the stable states of forced vibrations via the quantum activation mechanism. As seen from Fig. 1(a), with an accuracy to a factor ∼ 1/2 the switching rate W sw is determined by the probability to reach the saddle-point value g S of g(Q.P ). The switching rate is small, as switching requires that the oscillator makes many interlevel transitions with rates W mn smaller than the rates of transitions in the opposite direction, W nm . Therefore, before the oscillator switches, there is formed a quasistationary distribution over its states localized about the initially occupied extremum of g(Q, P ). This is similar to what happens in thermally activated switching over a high barrier [31]. However, in contrast to systems in thermal equilibrium, a modulated oscillator generally does not have detailed balance. Its statistical distribution has a simple Boltzmann form with temperature T e only for small damping and only close to the extrema of g(Q, P ). Therefore the standard technique developed for finding the switching rate in quantum equilibrium systems [44,45,46,47] does not apply. Also, even for T → 0 an oscillator modulated close to its eigenfrequency generally does not switch via tunneling (see [16,19,21,48,49,50,51,52] for the theory of tunneling switching for additive and parametric modulation). Switching via quantum activation is exponentially more probable.

Relation to chemical kinetics and population dynamics
For small scaled decay rate κ the switching rate W sw can be obtained from the balance equation (11). An approach to solving this equation was discussed earlier [15,18].
Here we provide a formulation that gives an insight into how the oscillator actually moves in switching and also makes a direct connection with the technique developed in chemical kinetics and population dynamics. The balance equation is broadly used in these areas. It describes chemical or biochemical reactions in stirred reactors (no spatial nonuniformity). The reactions can be thought of as resulting from molecular collisions in which molecules change, and if the collision duration is small compared to the reciprocal collision rate the kinetics is described by a Markov equation [53] ρ(X, τ ) = Here, X = (X 1 , X 2 , ...) is the vector that gives the numbers of molecules X i of different types i, and ρ is the probability for the system to be in a state with given X; W (X, r) is the rate of a reaction in which the number of molecules changes from X to X + r.
Typically, X i are large, X i ∝ N ≫ 1, where N is the total number of molecules. In contrast, the change of the number of molecules in an elementary collision is |r| ∼ 1, because it is unlikely that many molecules would collide at a time. Equation (15) is also often used in population dynamics, including epidemic models, cf. [54]. In this case the components of X give populations of different species. Since the number of molecules (population) is large, N ≫ 1, fluctuations are small on average. Disregarding fluctuations corresponds to the mean-field approximation. In this approximation one can multiply (15) by X and sum over X while assuming that the width of the distribution ρ(X) is small. This gives the equation of motion for the scaled mean number of molecules (population) Stable solutions of (16) give the stable states of chemical (population) systems. There may be also unstable stationary or periodic states. In population dynamics, an unstable stationary solution of (16) can be the state where one of the species goes extinct.
Equation (15) describes diffusion in the space of variables x. Along with small (∝ N −1/2 ) fluctuations around the stable states, this diffusion leads to rare large deviations (∼ O(N −1 ) in x-space) and to switching between the stable states. There is an obvious similarity between diffusion over the number of molecules and diffusion over quasienergy states of a modulated oscillator, but there are also some subtle differences, which we discuss below. There is also an obvious difference, with profound consequences: in the case of an oscillator the transition rates W mn (11), (12) are not limited to |m − n| ∼ 1.

The eikonal approximation
The role of the large number of molecules (population) in a modulated oscillator is played by the reciprocal effective Planck constant λ −1 , which determines the number of states localized about the extrema of g(Q, P ), cf. Fig. 1. For λ ≪ 1 it is convenient to switch from the state number n to the classical mechanical action I for the Hamiltonian orbits Q(τ ; g), P (τ ; g), which are described by equationsQ = ∂ P g(Q, P ),Ṗ = −∂ Q g(Q, P ), where ν(g) is the vibration frequency for given g [2πI gives the area of the cross-section of the surface g(Q, P ) in Fig. 1(a) by plane g =const]. One can show that, in spite of the nonstandard form of g(Q, P ), the semiclassical quantization condition has the familiar form I n ≡ I(g n ) = λ(n + 1/2). In the semiclassical approximation the rates of transitions between quasienergy states W mn become functions of the quasicontinuous variable I and can be written as W mn = W (I m , n − m). The dependence of W on I is smooth, as seen from (11) This equation shows how the oscillator is most likely to evolve. Using that the matrix element a mn in the expression (11) for the rate W mn is the (n − m)th Fourier component of function (2λ) −1/2 [Q(τ ; g m )+iP (τ ; g m )], one can show by invoking the Stokes' theorem that the time evolution of I is extremely simple, where I S is the value of I(g) for g approaching the saddle-point value g S from the side of the extremum of g(Q, P ) of interest; the values of I S are different on the opposite sides of g S . Equation (19) coincides with the result for the evolution of g for a classical modulated oscillator [3]. We note that the semiclassical approximation breaks down very close to the saddle point (in particular, the relation W (I m , n − m) ≈ W (I n , n − m) clearly ceases to apply), but the width of the corresponding range of I goes to zero as λ → 0. We now consider the quasistationary distribution ρ n about the initially occupied stable state. It is formed on times (κ|δω|) −1 ≪ t ≪ W −1 sw . To find ρ n far from the stable state we use the eikonal approximation [15,18], but in the form similar to that used in chemical kinetics and population dynamics [33]. We set and assume that |∂ I R| ≪ λ −1 . Then ρ n+r ≈ ρ n exp[−r∂ I R] for |r| ≪ λ −1 and, to the leading order in λ, the balance equation (11) becomes Equation (21)

Optimal switching trajectory
An advantageous feature of the formulation (21) is that it provides an insight into how the quantum oscillator evolves in large fluctuations that lead to occupation of quasienergy states far from the initially occupied extremum of g(Q, P ). Even though the diffusion over quasienergy states is a random process and different sequences of interstate transitions can bring the system to the given quasienergy state, the probabilities of such sequences are strongly different. Of physical interest is the most probable sequence, known as the optimal fluctuation. For classical fluctuating systems it has been understood theoretically and shown in experiment and simulations [33,55,56,57,58] that the evolution of the system in the optimal fluctuation, i.e., the optimal fluctuational trajectory is given by the classical trajectory of the auxiliary Hamiltonian system, which in the present case is described by equatioṅ The concept of the optimal fluctuational trajectory can be extended to the quantum oscillator. Such trajectory for the action variable I is well-defined, since any information of the oscillator phase is automatically erased and the range of the I values largely exceeds the quantum uncertainty in I, which is ∝ λ. Therefore the optimal fluctuational trajectory I(t) can be measured in the experiment in the same way as it is done in classical systems.
In (22) we have set R(0) = 0 and thus ignored the normalization factor (an analog of the reciprocal partition function) in the expression (20) for ρ n . This factor leads to a correction ∝ λ to R(0). Since, to logarithmic accuracy, the switching rate is determined by the probability of approaching the saddle-point value I S of I, we have Parameter R A plays the role of the effective activation energy for switching via quantum activation, with the effective Planck constant λ replacing the temperature in the conventional expression for thermally activated switching. Optimal fluctuations away from the extremum of g(Q, P ) are described by optimal trajectories that emanate from I = 0, which is reflected in (22). The value of the momentum p I ≡ ∂ I R for I → 0 on the trajectory can be found by noticing that the distribution over quasienergy near the extremum of g(Q, P ) is of the form of the Boltzmann distribution with effective temperature T e , and thus R ∝ I/T e ; from (14) p I = ln[(n e + 1)/n e ] for I → 0. Then from (21),İ = 2κI on the optimal trajectory for I → 0. As expected, the system moves along the optimal fluctuational trajectory away from the stable state of fluctuation-free dynamics.
The facts that p I = 0 at the starting point of the optimal trajectory and that the state I = 0 lies on the boundary of the available values of I are connected with each other and present a distinctive feature of the oscillator dynamics. In chemical kinetics and population dynamics usually stable states lie in the middle of the space of dynamical variables X. The probability distribution has a Gaussian maximum at such X, and then the momentum on the optimal trajectory is equal to zero [28,33]. The states (I = 0, p I = 0) and (I = 0, p I = ln[(n e + 1)/n e ]) are stationary points of the Hamiltonian H(I, p I ). From (14),İ =ṗ I = 0 at these points. The motion of the system near these points is exponential in time and is shown in Fig. 3. Figure 3(a) shows the mean-field (fluctuation-free) trajectory I(t) and the optimal trajectory I(t) obtained numerically from equations (19) and (22), respectively. An interesting feature of the considered model of the modulated quantum oscillator is that it satisfies the detailed balance condition for T = 0 and thusn = [exp( ω 0 /k B T )−1] −1 = 0 [27]. This is seen from the explicit expression for the rates (11) and (12), as forn = 0 they meet the familiar detailed balance condition W n n+k /W n+k n = exp(−k/ξ n ) [the explicit form of ξ n ≡ ξ(I n ) follows from (12)]. Therefore p I = 1/ξ(I), and one can show from (22) thatİ = 2κI. As a consequence, the optimal fluctuational trajectory I(t) is the time-reversed mean-field trajectory I(t). This is a generic feature of classical systems with detailed balance, see [55]. Our results show that the symmetry also holds in quantum systems provided the notion of a trajectory is well-defined.
Of special interest is the vicinity of the saddle-point value of the action variable I S , see Fig. 3. In a dramatic distinction from chemical kinetics, there is no slowing down of I(t) near I S . The quantity I S is a boundary value of I for states localized about a given extremum of g(Q, P ) in Fig. 1. Functionsİ andİ are discontinuous there. This is an artifact of the balance equation approximation, which applies in the weak Figure 3. (a) The mean-field (fluctuation free) and optimal fluctuational trajectories of the action variables. Because the system has detailed balance forn = 0, the optimal trajectory in this case is the time-reversed mean-field trajectory. The data refer to the trajectories for the local maximum of g(Q, P ) in Fig. 1 for β = 0.035. The shape of the trajectory changes discontinuously wheren becomes nonzero; the trajectories for n = 0 andn → 0 coincide only for I < In. The limitn → 0 is taken with the constraint n ≫ λ damping limit where the dimensionless frequency ν(g) ≫ κ. For g → g S the frequency ν(g) → 0, and the approximation breaks down. With account taken of decay, in the region of bistability the oscillator has a "true" unstable stationary state in the neglect of fluctuations. Both the mean-field trajectory and the optimal trajectory in phase space are moving away/approaching this state exponentially in time, cf. [15], but the region of I where it happens is very narrow for small κ.

Fragility in the problem of large rare fluctuations
A striking feature of optimal fluctuational trajectories obvious from Fig. 3 is that these trajectories have different shapes depending on whether the oscillator Planck number is n = 0 orn > 0. The discontinuous with respect ton change of the trajectories and the associated change of the logarithm of the distribution R(I) and of the activation energy for switching R A show the fragility of the detailed-balance solution forn = 0 [15,18]. It has been found that the fragility also emerges in a very different type of problem, the problem of population dynamics described by equation (15) [28]. In particular, the well-known result for the rate of disease extinction in the presence of detailed balance [59,60,61,62] can change discontinuously with the varying elementary rates W (X, r) as the detailed balance is broken.
We now show that the condition for the onset of fragility proposed in [28] applies also to the modulated oscillator, even though the divergence it reveals shows up in a different fashion. The condition relies on the expression for the switching exponent. Similar to how it was done for the oscillator, this exponent can be found by seeking the solution of the master equation (15) in the eikonal form ρ(X) = exp[−NR(x)]. To the leading order in 1/N, the problem is then mapped onto Hamiltonian dynamics of an auxiliary system with mechanical actionR(x). From (15), the Hamiltonian of the auxiliary system is [as before, we use that W (X − r, r) ≈ W (X, r)]. If the system is initially near a stable state x a [a stable solution of (16)],R(x) is determined by the Hamiltonian trajectories that emanate from x a . From (24),R(x) = x xa pdx. The rate of switching from x a (or extinction, in the extinction problem) is W sw ∝ exp(−NR A ). Similar to the quantum oscillator,R Here x S is the saddle point of the deterministic dynamics (16); it can be shown that it is the Hamiltonian trajectory that goes to the saddle point that provides the switching or extinction exponentR A , cf. [3,33,63]. Both x a and x S are stationary points of the Hamiltonian H, and the integral over time in (25) goes from −∞ to ∞. This is a significant distinction from the modulated oscillator problem; there equations (22) and (23) for the activation exponent can be written as where we set the instant where I(t) reaches I S on the optimal trajectory to be t = 0. A small change of the reaction rates W (X, r) → W (X, r) + ǫW (1) (X, r) (ǫ ≪ 1) leads to the linear in ǫ change of the Hamiltonian, H → H + ǫH (1) , as seen from (24). The action is then also changed. To the first order in ǫ,R A →R A +ǫR (1) A . The correction term is given by a simple expression familiar from the Hamiltonian mechanics [2], where the integral is calculated along the unperturbed trajectory x(t), p(t). In the extinction problem the integral (26) can diverge at the upper limit, t → ∞. This is because in this problem p(t) remains finite for t → ∞, and therefore if w (1) (x S , r) is nonzero, H (1) = 0 for t → ∞. The divergence indicates the breakdown of the perturbation theory; in the particular example studied in [28], for ǫ → 0 the change of R A was ∼R A . For the modulated oscillator, the role of the small parameter ǫ is played by the Planck numbern. If w (0) (I, r) is the transition rate forn = 0, then from (11) the thermally induced term in the transition rate has the formnw (1) (I, r) = n[w (0) (I, r) + w (0) (I, −r)]. Where the perturbation theory applies, the correction to the effective activation energy of switching reads As we saw, in contrast to reaction/population systems, p I = 0 for t → −∞. However, w (1) (I, r) ∝ I ∝ exp(2κt) for t → −∞, therefore (27) does not diverge for t → −∞.
There is also no accumulation of perturbation for large t, as the integral goes to t = 0. Therefore the cause of the fragility should be different from that in population dynamics/reaction systems. As mentioned earlier, in contrast to reaction systems, for the oscillator the values of r in (27) can be large. Then the correction R (1) A can diverge because of the divergence of the sum over r. This happens if on the optimal trajectory w (1) (I, r) decays with r slower than exp(rp I ). From (12), w (1) (I, r) decays with r exponentially; in particular, w (1) (I, r) ∝ exp[−2rν(g)τ * (g)] for r ≫ 1. The region of the values of p I where r w (1) (I, r) exp(rp I ) remains finite is shown in Fig. 3(b). As seen from this figure, the value of p I on then = 0-trajectory can be too large for the sum over r to converge. Then the perturbation theory becomes inapplicable. The trajectory followed in switching changes discontinuously wheren changes fromn = 0 ton > 0. The probability distribution also changes discontinuously. We note that |p I | ∼ 1 ≪ λ −1 on the optimal fluctuational trajectory, which justifies the approximation (21) that underlies the above analysis. It is clear that the optimal fluctuational trajectory I(t) corresponds to the optimal fluctuational trajectory of the quasienergy g(t), since the I and g variables are related by ∂ g I = ν −1 (g).
The instanton approximation relies on the assumption that the mean square fluctuations provide the smallest scale in the problem, similar to the wavelength in the WKB approximation [42]. If the system is perturbed and the perturbation is small, it can be incorporated into the prefactor of the rate of rare large fluctuations. If the perturbation is still small but exceeds the small parameter of the theory, it can be incorporated into the instanton Hamiltonian and leads to a correction to the exponent of the rare event rates. This correction is generically linear in the perturbation. However, this is apparently not a universal behavior, as the unperturbed solution can be fragile with respect to a perturbation. So far the fragility has been found in cases where the perturbation breaks the time-reversal symmetry.
An important problem is the crossover between the instanton solutions without and in the presence of the perturbation. For a modulated quantum oscillator it was recently addressed in [41] (but the most probable fluctuational trajectories were not studied in this paper). The analysis [41] shows that the very instanton approximation breaks down by thermal fluctuations, function ∂ I R is not smooth forn > 0, rather it displays a kink. The threshold for the onset of this behavior is exponentially low inn, with | lnn| λ −1 . It corresponds to the regime where the rate of transitions between oscillator states induced by absorption of thermal excitations, which is ∝n, becomes comparable with the switching rate W sw calculated forn = 0. The region where the instanton approximation is inapplicable extends ton λ 3/2 . This is why we indicate that the optimal trajectories in Fig. 3 forn → 0 correspond to vanishingly smalln compared to the small parameter of the theory λ, yetn λ 3/2 .

Nonresonant modulation: a brief summary
Much attention has attracted recently the possibility of cooling mesoscopic oscillators, and the whole new area, the cavity optomechanics, has emerged, see [64] for a recent review. The cooling is performed by nonresonant modulation, with frequency ω F significantly different from the oscillator frequency ω 0 . The very idea of cooling different types of quantum systems by a high-frequency field goes back to the mid-70s [37,65,66], about the same time when the laser cooling of atomic motion was proposed [67,68]. The change of the distribution can be understood from Fig. 4 [40,65]. It refers to a system coupled by the modulating field to another system, which can be a thermal bath or a mode with a relaxation time much shorter than that of the system of interest, so that it serves effectively as a narrow-band thermal reservoir. The modulation provides a new channel of relaxation for the relatively slowly relaxing system of interest. Figure 4 indicates possible transitions between the states of the system accompanied by energy exchange with the thermal reservoir. For example in (a), a transition of the system from the excited to the ground state is accompanied by a transition of the reservoir to the excited state with energy ω b = (ω 0 + ω F ), with the energy deficit compensated by the modulation. On the other hand, a transition of the system from the ground to the excited state requires absorbing an excitation in the thermal reservoir, which is possible only when such excitation is present in the first place. The ratio of the state populations of the system is determined by the ratio of the rates of transitions up and down in energy, and thus by the population of the excited states of the thermal reservoir with energy ω b . If the corresponding process is the leading relaxation process, the effective temperature of the system becomes T * = (ω 0 /ω b )T . It means there occurs effective cooling for ω 0 ≪ ω b . Similarly, for ω 0 ≫ ω b the modulation leads to heating of the system, see Fig. 4(b). In the case sketched in Fig. 4 (c), the induced transitions from the ground to the excited state of the system are more probable then from the excited to the ground state, which leads to a negative effective temperature for strong modulation.
In the case of an oscillator, the system has many levels, but the above picture still applies. The goal of this section is to outline and compare different microscopic mechanisms of the coupling of the oscillator to the modulation and the bath. The unexpected feature is that the distribution of the oscillator over its Fock states can be of the Boltzmann form with an effective temperature determined by the strength and frequency of the modulation [37]. However, this is the case only provided the major mechanism of oscillator relaxation in the absence of modulation is the conventional mechanism (6), which in a phenomenological classical description of oscillator dynamics corresponds to a friction force proportional to the oscillator velocity.
A simple model of the modulation-induced dissipation is where the external field parametrically modulates the coupling of the oscillator to a thermal bath. The coupling Hamiltonian is Figure 4. Modulation-induced relaxation processes leading to cooling (a), heating (b), and population inversion (c); ω 0 , ω F , and ω b are the frequency of the system (the oscillator, in the present case), the modulation frequency, and the frequency of the mode (or a thermal bath excitation) to which the oscillator is coupled by the modulation, respectively; the relaxation time of the mode is much shorter than that of the oscillator. Strong modulation imposes on the oscillator the probability distribution of the fast-decaying mode in (a) and (b) and leads to population inversion in (c). If, in the absence of modulation, oscillator relaxation is described by the standard model (6), (7), the distribution over the Fock states in the presence of modulation is of the form of the Boltzmann distribution [37]; in (c) the distribution over low-energy Fock states is described by negative temperature and the oscillator vibrates close to its eigenfrequency.
Here h depends on the variables of a thermal bath, or it can be the coordinate of a comparatively quickly decaying mode coupled to a thermal bath. The interaction (28) has the same structure as the interaction (6), except that it can lead to decay processes with the energy transfer (ω 0 ± ω F ), cf. Fig. 4(a) and (b). Therefore the structure of the master equation for the oscillator should not change, but the decay parameters and the Planck numbers of excitations created in decay should change appropriately. The interaction can also lead to decay processes with energy transfer ω F − ω 0 , for the appropriate modulation frequencies. In this case absorption of bath excitations is accompanied by oscillator transitions down in energy. Respectively, in the master equation (7) in the expression for the rates of transitions due to excitation absorption one has to formally replacen(ω 0 ) →n(ω 0 − ω F ) = −n(ω F − ω 0 ) − 1, which means that the friction coefficient becomes negative.
The above qualitative arguments can be confirmed by a formal analysis similar to that in [37]. It shows that in the RWA the master equation for the oscillator with account taken of the modulation-induced relaxation processes has the form (7) with the relaxation parameter Γ and the Planck numbern replaced by Γ F = Γ + Γ + + Γ − − Γ inv andn F , (29) where Here, Γ ± give the rates of transitions at frequencies ω 0 ± ω F , which correspond to the processes sketched in Fig. 4(a) and (b); Γ inv gives the rate of processes sketched in Fig. 4(c), where excitation of the oscillator is accompanied by excitation of the thermal bath. If these processes dominate, they lead to vibrations of the oscillator at frequency ≈ ω 0 , with amplitude determined by other mechanisms of losses [40]. Parameters Γ − and Γ inv in (30) refer to the cases where ω 0 > ω F and ω 0 < ω F , respectively. From (29) and (30), the probability distribution of the oscillator is characterized by effective temperature A similar behavior occurs if the modulation is performed by an additive force A cos ω F t, but the interaction with the narrow-band thermal reservoir is nonlinear in the oscillator coordinate, H b . This case was considered in [37]. It reduces to the above formulation if one makes a canonical transformation is replaced with h (2) b . The analysis of cooling of a vibrating mirror in an optical cavity can be also often mapped onto the analysis for the interaction (28). A quantum theory in this case was developed in [69,70]. It considers an oscillator (the mirror) coupled to a cavity mode driven by external radiation. If the radiation is classical, in the appropriately scaled variables the coupling and modulation are described by Hamiltonians H where q and q b are the coordinates of the mirror and the mode. In cavity optomechanics one usually writes H ; the following discussion immediately extends to this form of the interaction.
In the absence of coupling to the mirror, the cavity mode is a linear system, hence q b (t) = q b 0 (t) + [χ b (ω F ) exp(−iω F t) + c.c.]A/2, where q b 0 (t) is the mode coordinate in the absence of modulation and χ b (ω) is the susceptibility of the mode [70]. The coupling H (m) i in the interaction representation then has a cross-term ∝ q 0 (t)q b 0 (t) exp(±iω F t).
Since the cavity mode serves as a thermal bath for the mirror, this term is fully analogous to H

Conclusions
It follows from the results of this paper that a modulated nonlinear oscillator displays a number of quantum fluctuation phenomena that have no analog in systems in thermal equilibrium. Oscillator relaxation is accompanied by a nonequilibrium quantum noise. It leads to a finite-width distribution of the oscillator over its quasienergy states even for the bath temperature T → 0. For resonant modulation, the distribution is Boltzmannlike near the maximum. We have discussed the recent experiment that confirmed this prediction and the effect of oscillator damping on the outcome of a sidebandspectroscopy based measurement of the distribution.
The quantum noise also leads to large rare events that form the far tail of the distribution of the oscillator over quasienergy and to switching between the coexisting states of forced vibrations. We have developed an approach to the analysis of the distribution tail and the switching rate, which makes a direct connection with the analysis of the corresponding problems in chemical and biological systems and in population dynamics. We show that, in a large deviation, the quasienergy of an underdamped oscillator most likely follow a well-defined real trajectory in real time. This trajectory is accessible to measurement. For T = 0, where the oscillator has detailed balance, the most probable fluctuational trajectory is the time-reversed trajectory of the fluctuation-free (mean-field) relaxation of the oscillator to the stable state. Thermal fluctuations break the detailed balance condition and, even where the thermal Planck numbern is small compared to the effective Planck constant, lead to ann-independent change of the most probable fluctuational trajectory. We show that the criterion of the fragility, i.e., of a discontinuous change of the optimal fluctuation trajectory with the varying parameter can be formulated in a general form, that applies both to reaction systems with classical fluctuations and to the modulated quantum oscillator.
An interesting effect of nonresonant modulation of the oscillator is that it can significantly change the oscillator distribution over the Fock states, leading to heating, cooling, or excitation of self-sustained vibrations depending on the modulation frequency and the coupling to the thermal bath or a mode with a relaxation time shorter than that of the oscillator. We show that different coupling and modulation mechanisms can be described in a similar way and derive explicit expressions for the effective decay rate and temperature of a modulated oscillator.

Acknowledgments
We are grateful to P. Bertet and V. N. Smelyanskiy for the discussion. This research was supported in part by the ARO, grant W911NF-12-1-0235, and the Dynamics Enabled Frequency Sources program of DARPA. VP also acknowledges support from the ERC Starting Grant OPTOMECH. (we provide the definition of the relevant correlator for arbitrary operators K, L). The major contribution to this power spectrum comes from small-amplitude fluctuations about the stable states of forced vibrations; in the limit of weak damping these states correspond to the extrema of function g(Q, P ). We will disregard fluctuation-induced transitions between the stable states and calculate the correlator (A.1) for each of these states separately; the averaging ... then means averaging over small-amplitude fluctuations about the corresponding state. However, we will not limit ourselves to small damping; moreover, we will assume that the scaled damping rate κ ≫ λ, so that the spectra do not display the fine structure related to the nonequidistance of the levels g n near the stable states [29,38]. Operatorn(t) smoothly depends on time: it does not have fast-oscillating factors ∝ exp(±iω 0 t). However, for small κ,n(t) contains terms which oscillate at frequencies ∼ ν 0 δω and ∼ 2ν 0 δω with ν 0 being the dimensionless frequency of vibrations about the considered stable state (13). To see this, we first note that classical motion about the stable state (Q a , P a ) is described by linearized equations (9) for δQ = Q − Q a , δP = P − P a . For small κ this motion is decayed vibrations [3] with dimensionless frequency (the imaginary part of the eigenvalues of the equations for δQ, δP ) ν a = κ 2 + 3r 4 a − 4r 2 a + 1 1/2 , r 2 a = Q 2 a + P 2 a . (A.2) From (13) and (A.2), ν a → ν 0 for κ → 0. We now write the operators of the oscillator as a = a a + δa, a a = (2λ) −1/2 (Q a + iP a ), δa = (2λ) −1/2 (δQ + iδP ),n ≈ a † a a a + a † a δa + δa † a a + δa † δa. This immediately shows that, indeed,n oscillates at dimensional frequencies ν a δω and, with smaller amplitude (quadratic in δQ, δP ), 2ν a δω. From (A.1) and (A.3), the leading term in the power spectrum ofn is Re n,n ω ≈ (r 2 a /2λ)|δω| −1 Φ nn (ω), Φ nn (ω) = |δω| Re [ δa, δa † ω + δa † , δa ω ] (A.4) Functions δa, δa † ω and δa † , δa ω were found earlier [71] (see also [38] where the current notations were adopted). Using their explicit form, one can show that, for small decay rate, κ ≪ ν a , function Φ nn (ω) has two Lorentzian peaks. They are located at ≈ ±ν a |δω| and have halfwidth κ|δω|.
Examples of the spectra Φ nn (ω) are shown in Fig. 2(b). For the chosen parameters the spectra have two well-resolved peaks. A straightforward but somewhat tedious calculation shows that the ratio of the heights of the peaks of Φ nn (ω) approaches n e /(n e + 1) for κ/ν a → 0, wheren e is the effective Planck number for vibrations about the stable state (14). Therefore by measuring this ratio one can reveal and quantitatively characterize the effect of quantum heating. However, for the parameters in Fig. 2(b), even though the peaks are well resolved, the ratio of their heights is different from n e /(n e + 1). This ratio as a function of the scaled intensity of the modulating field β is shown by triangles in Fig. 2(a).
In the experiment [30] the occupation of excited quasienergy states of the oscillator was detected by attaching the oscillator to a two-level system (qubit). There was applied an extra field ∝ F q cos ω q t at frequency ω q close to the transition frequency of the qubit ω ge , and then the resulting population of the excited qubit state was measured. The frequencies were chosen in such a way that |ω 0 − ω q | ≫ |ω q − ω ge |, therefore transitions between the qubit states accompanied by transitions between the Fock states of the oscillator have negligible probability. However, one can think that the oscillator modulates the coupling of the qubit to the field F q . Then phenomenologically one can write the effective Hamiltonian of the driven qubit as where σ z , σ ± = σ x ± iσ y are Pauli matrices.