Nonlocal random walk over Floquet states of a dissipative nonlinear oscillator

TL;DR

This paper investigates how a dissipative nonlinear oscillator driven near triple its eigenfrequency exhibits nonlocal transitions between Floquet states, leading to a random walk in state space and unique Floquet dynamics features.

Contribution

It introduces a model of nonlocal Floquet state transitions in a driven oscillator, revealing new insights into multi-state Floquet dynamics and quantum switching behavior.

Findings

Transition rates for Floquet states are derived.

The stationary distribution involves remote state transitions.

Switching rates scale with system parameters.

Abstract

We study transitions between the Floquet states of a periodically driven oscillator caused by the coupling of the oscillator to a thermal reservoir. The analysis refers to the oscillator that is driven close to triple its eigenfrequency and displays resonant period tripling. The interstate transitions result in a random ``walk'' over the states. We find the transition rates and show that the walk is nonlocal in the state space: the stationary distribution over the states is formed by the transitions between remote states. This is to be contrasted with systems in thermal equilibrium, where the distribution is usually formed by transitions between nearby states. The analysis of period tripling allows us to explore the features of the multi-state Floquet dynamics including those missing in the previously explored models of driven oscillators such as the absence of detailed balance for low temperatures. We use the results to study switching between the period-3 states of the oscillator due to quantum fluctuations and find the scaling of the switching rates with the parameters.

Equations

$$H$$

$$|\delta \omega |,|\gamma |⟨q^2⟩\ll {\omega }_0,\delta \omega =\frac{1}{3}{\omega }_F-{\omega }_0.$$

$$U^{†}(t)qU(t)=C(Q\cos \varphi +P\sin \varphi ),$$

$$U^{†}(t)pU(t)=-C\frac{{\omega }_F}{3}(Q\sin \varphi -P\cos \varphi ),\varphi ={\omega }_{F}t/3,$$

$$[Q,P]=i\lambda ,\lambda =\frac{27\gamma \hslash }{8{\omega }_F^2|\delta \omega |}.$$

$$g(Q,P)=\frac{1}{4}[Q^2+P^2-\mathrm{sgn}(\delta \omega ){]}^2-\frac{1}{3}f(Q^3-3PQP)$$

$$f=F_0/(8{\omega }_F\gamma |\delta \omega |{)}^{1/2}.$$

$$i\lambda \frac{\partial }{\partial \tau }|\varphi ⟩=g(Q,-i\lambda {\partial }_Q)|\varphi ⟩.$$

$${\hat{N}}_3{\varphi }_n^{(k)}=\exp (-2\pi ik/3){\varphi }_n^{(k)}.$$

$${\epsilon }_n^{(k)}=\frac{\hslash |\delta \omega |}{\lambda }{g}_n^{(k)}+\frac{2\pi k}{3}\hslash {\omega }_F.$$

$$g$$

$$\displaystyle-\frac{f}{3}\sqrt{2\lambda}(a^{\dagger 3}+a^{3})\Bigr{]}+\frac{1}{4}(1+\lambda^{2}).$$

$${\psi }_{\nu }(Q;n)=C_n^{(\nu )}({\partial }_{P}g{)}^{-1/2}\exp [i{S}^{(\nu )}(Q,g_n)/\lambda ],$$

$${\partial }_Q{S}^{(\nu )}=P^{(\nu )}(Q,g_n)$$

$${}_{\nu }⟨n|g|n{⟩}_{\nu }=g_n,{}_{\nu }⟨n|g|n{⟩}_{\nu \pm 1}=J_{\nu \pm }(g_n).$$

$${\varphi }_n^{(k)}(Q)\approx 3^{-1/2}\sum _{\nu }{\psi }_{\nu }(Q;n)e^{-2\pi i\nu k/3},$$

$$g_n^{(k)}-g_n\approx 2\mathrm{Re}[J_{0+}(g_n)\exp (-2\pi ik/3)].$$

$${\psi }_{\nu ;\mathrm{W}}(Q;n)=3^{-1/2}\sum _k{\varphi }_n^{(k)}(Q)\exp (2\pi i\nu k/3).$$

$$\frac{d}{d\tau }\rho =$$

$$\mathcal{D}[a]\rho =$$

$$+\bar{n}(a{a}^{†}\rho -2a^{†}\rho a+\rho a{a}^{†}).$$

$$\frac{d}{d\tau }{\rho }_n=-\sum _{n^{\prime }}{W}_{n{n}^{\prime }}{\rho }_n+\sum _{n^{\prime }}{W}_{n^{\prime }n}{\rho }_{n^{\prime }},$$

$$W_{n{n}^{\prime }}=2\kappa \left[(\bar{n}+1){\left|⟨n^{\prime }|a|n⟩\right|}^2+\bar{n}{\left|⟨n|a|n^{\prime }⟩\right|}^2\right]$$

$$⟨n+m|\hat{a}|n⟩\equiv a_m(g_n).$$

$$a_m(g)=\frac{\omega (g)}{2\pi }{\int }_{-\pi /\omega (g)}^{\pi /\omega (g)}d\tau \exp [-im\omega (g)\tau ]a(\tau ;g).$$

$$\frac{d}{d\tau }Q={\partial }_{P}g,\frac{d}{d\tau }P=-{\partial }_{Q}g.$$

$$W_{nn+m}=2\kappa [(\bar{n}+1)|a_m(g_n){|}^2+\bar{n}|a_{-m}(g_n){|}^2].$$

$$|a_{m<0}(g_n)|\approx$$

$$\times \exp [-|m|\omega (g_n){\tau }_{\infty }(g_n)].$$

$$|a_{m>0}(g_n)|\approx$$

$$\times \exp [-m\omega (g_n){\tau }_{\infty }(g_n)].$$

$$\omega (g){\tau }_{\infty }\approx -\frac{1}{2}\ln \Delta g,\Delta g\equiv \frac{g-g_{\min }}{g_s-g_{\min }}\ll 1,$$

$$\frac{W_{n{n}^{\prime }}}{W_{n^{\prime }n}}=\frac{W_{n{n}^{\prime \prime }}}{W_{n^{\prime }{n}^{\prime \prime }}}\frac{W_{n^{\prime \prime }{n}^{\prime }}}{W_{n^{\prime \prime }n}}.$$

$$W_{n+mn}>W_{nn+m},m>0,$$

$$g\approx g_{\min }+\frac{1}{2}{g}_{PP}{P}^2+\frac{1}{2}{g}_{QQ}(Q-Q_0{)}^2,$$

Full text

I Introduction

The distribution over the eigenstates of a quantum system weakly coupled to a thermal reservoir is formed by the coupling-induced interlevel transitions. In thermal equilibrium, the transitions are balanced, the rates of transitions $m\to n$ and $n\to m$ weighted with the state populations are equal for any states $m,n$ Einstein (1917). This detailed balance condition does not hold away from thermal equilibrium. When the distribution over the states of a nonequilibrium system is stationary, it means that the overall probability to come to a state from all other states is equal to the overall probability to leave this state. In this case, if the distribution over the states is steep, the transitions to a given state from remote states may become important. This is despite the fact that the rates of such transitions, which are determined by the overlapping of the corresponding wave functions, are exponentially small.

One can think of the interstate transitions due to the coupling to a reservoir as a random walk over the states. Where the transitions between remote states are important, the walk becomes nonlocal, in the state space.

An important group of nonequilibrium system where the locality of the walk over the states can be explored is provided by periodically driven systems. The interest in such system has surged in recent years in various areas of physics and in various contexts, cf. Peano and Thorwart (2006); Martin et al. (2017); Weinberg et al. (2017); Desbuquois et al. (2017); Schmidt (2018); Seetharam et al. (2018); Lohse et al. (2018) and references therein. Periodically driven systems have well-defined quantum states, the Floquet or quasienergy states  Shirley (1965); Zel’dovich (1967); Ritus (1967), which form a complete set. It has been appreciated Kohn (2001) that, when a periodically driven system is brought in contact with a thermal reservoir, there is no detailed balance, and the resulting stationary distribution over the Floquet states is complicated.

A remarkable and important exception to the lack of detailed balance in Floquet systems is provided by a nonlinear oscillator driven close to its eigenfrequency or parametrically modulated close to twice the eigenfrequency Drummond and Walls (1980); Kryuchkyan and Kheruntsyan (1996). For the basic model of such an oscillator, that includes the quartic nonlinearity of the potential and a simple dissipation mechanism, the oscillator has detailed balance, if the temperature of the thermal reservoir is zero. The detailed balance breaks down for nonzero temperature Dykman and Krivoglaz (1979); Dykman and Smelyanskii (1988); Marthaler and Dykman (2006). However, the very occurrence of it for $T=0$ is fascinating. Since it happens for two major models, one may ask whether this is a generic feature of periodically driven dissipative oscillators.

In this paper, we study the rates of transitions between the Floquet states and the stationary distribution over these states for a nonlinear oscillator driven close to triple its eigenfrequency $\omega_{0}$. The quantum coherent behavior of such an oscillator displays peculiar features such as the nontrivial geometric phase and the oscillations of the Floquet wave functions in the classically inaccessible region Guo et al. (2013a); Zhang et al. (2017). An observation of period tripling in the quantum regime with a coupled-modes superconducting cavity has been already reported Svensson et al. (2017). Here we explore the dissipative dynamics of this system and the associated quantum fluctuations. It turns out that the analytical tools developed for the previously explored oscillator models Dykman and Smelyanskii (1988); Marthaler and Dykman (2006); Peano and Dykman (2014) do not apply in the present case. Moreover, even for $T=0$ the oscillator does not have detailed balance. This suggests that the detailed balance for resonant and parametric driving could be an artifact.

Another important feature of the oscillators driven close to the eigenfrequency or modulated close to twice the eigenfrequency is that their stationary distribution over the Floquet states is formed locally, i.e., by transitions between a few close states. This is the case both at zero temperature and above a certain very low temperature set by $\hbar$ weighted with a combination of the nonlinearity and driving parameters Guo et al. (2013b); Peano and Dykman (2014). This again poses the question of whether the locality is a generic property of the $T=0$- stationary distribution of a driven oscillator. Answering this question for period tripling is one of the central topics of the present paper.

Classically, an oscillator that displays period tripling has four stable states: three states of period-three vibrations, with the same amplitude but the phases that differ by $2\pi/3$, see Fig. 1 (a), and also the quiet state with no vibrations excited. Quantum fluctuations lead to switching between these states. For the driving at frequency close to $\omega_{0}$ and $2\omega_{0}$ it was found that, unless the dissipation rate is exponentially small, the switching occurs via transitions over the quasienergy barrier even for $T=0$. The mechanism was called quantum activation Marthaler and Dykman (2006). The period-3 vibrations are different in several respects. One has to find how the oscillator switches between the period-3 states and between these states and the quiet state. In particular, the oscillator can possibly switch directly between the period-3 states or first go to the quiet state and then switch from this state.

In Sec. II below, for the oscillator that displays resonant period tripling, we provide the Hamiltonian in the rotating wave approximation (RWA) and discuss the Floquet eigenstates and the Wannier-type states, which are relevant where the wells of the RWA Hamiltonian have many bound states. In Sec. III we develop the WKB theory of the rates of dissipation-induced transitions between the intrawell states and demonstrate the breaking of the detailed balance even for $T=0$. In Sec. IV we study the distribution over the states using an eikonal-type approximation. In this approximation the distribution is formed by transitions between nearby states. Sec. V is one of the central parts of the paper. Here we analyze the breaking of the local approximation for $T=0$. We also show that, already for small Planck numbers of the oscillator, the locality is restored. The results are compared with numerical analysis based on calculating the transitions rates using the Wannier-type wave functions. In Sec. VI we study switching between different wells of the RWA Hamiltonian in the regime of weak damping and also in the vicinity of the bifurcation points where the period-3 states of the oscillator emerge. Sec. VII contains conclusions. The Appendices describe technical details of the calculations.

II The model and the rotating wave approximation

Period tripling generically occurs in strongly driven nonlinear systems. However, for weakly damped oscillators it may emerge already for a relatively weak drive provided the frequency of the driving force $\omega_{F}$ is close to triple the oscillator eigenfrequency $\omega_{0}$ Nayfeh and Mook (2004). Several aspects of the quantum formulation of the problem in the absence of dissipation were discussed previously Zhang et al. (2017). To set the scene, here we briefly summarize and generalize the formulation.

A simple Hamiltonian of the oscillator that displays resonant period tripling is

[TABLE]

where $p$ and $q$ are the momentum and coordinate of the oscillator. The conditions of resonant driving and comparatively weak nonlinearity are

[TABLE]

Without loss of generality, we set $\gamma>0$ and $F_{0}>0$. The analysis of Zhang et al. (2017) referred to the case $\delta\omega>0$. An alternative scaling that allows one to consider the case $\delta\omega=0$ is considered in Appendix A.

To study the oscillator dynamics near its eigenfrequency, we switch to the rotating frame at one third of the driving frequency via the unitary transformation $U(t)=\exp(-i\omega_{F}ta^{\dagger}a/3)$, with $a=[(\omega_{F}q/3)+ip](3/2\hbar\omega_{F})^{1/2}$. We then introduce slowly varying in time coordinate $Q$ and momentum $P$ of the oscillator in the rotating frame,

[TABLE]

where $C=(8\omega_{F}|\delta\omega|/9\gamma)^{1/2}$. The operators $Q,P$ satisfy the commutation relation

[TABLE]

Parameter $\lambda$ plays the role of the Planck constant in the rotating frame. It is dimensionless, and in what follows we use the WKB approximation, which applies for $\lambda\ll 1$. In terms of the raising and lowering operators, $Q=(\lambda/2)^{1/2}(a+a^{\dagger})$ and $P=-i(\lambda/2)^{1/2}(a-a^{\dagger})$.

In the rotating wave approximation (RWA), the motion in the rotating frame is described by the canonically transformed Hamiltonian $U^{\dagger}HU-i\hbar U^{\dagger}\dot{U}=(\hbar|\delta\omega|/\lambda)g(Q,-i\lambda\partial_{Q})$, where

[TABLE]

Here, $g(Q,P)$ as a function of the operators $Q,P$ is the scaled dimensionless Hamiltonian of the oscillator in the rotating frame calculated in the RWA. It depends on a single parameter, the scaled field amplitude $f$, and does not depend on time.

The Hamiltonian $g(Q,P)$ as a function of the classical coordinate and momentum is plotted in Fig. 1 (b). In the case shown in the figure, $g(Q,P)$ has three local minima at a nonzero $Q^{2}+P^{2}$. For the frequency detuning $\delta\omega>0$ the onset of the minima has no threshold in the scaled field amplitude $f$, whereas for $\delta\omega<0$ the minima emerge for $f>2$, see Appendix A. In the classical limit and in the presence of weak dissipation, these minima correspond to the three period-3 vibrational states in the laboratory frame, cf. Fig. 1 (a).

The Hamiltonian $g(Q,P)$ also has a local maximum or minimum at $Q=P=0$ depending on whether $\delta\omega>0$ or $\delta\omega<0$. It corresponds to a classically stable quiet state with no vibrations. The dissipative classical dynamics is outlined in Appendix B.

II.1 The quasienergy spectrum

The Schrödinger equation in the rotating frame in terms of the scaled RWA Hamiltonian $g$ reads

[TABLE]

We have introduced here the dimensionless time $\tau=|\delta\omega|t$. The Hamiltonian $g$ is invariant with respect to a rotation in the $Q,P$ phase plane by $2\pi/3$. Such rotation is described by the operator $\hat{N}_{3}=\exp(-2\pi ia^{+}a/3)$, which commutes with $g$. Therefore the eigenfunctions $\phi_{n}^{(k)}$ of $g$ are also eigenfunctions of $\hat{N}_{3}$,

[TABLE]

The superscript $k=0,1,2$ here determines the eigenvalue of the rotation operator $\hat{N}_{3}$. The subscript $n$ is the other number that enumerates the eigenstates of the Hamiltonian $g(Q,-i\lambda\partial_{Q})$; this Hamiltonian has an infinite countable set of eigenfunctions and eigenvalues.

With the account taken of the interrelation between the oscillator Hamiltonian in the laboratory frame (1) and $g$, the eigenvalues $g_{n}^{(k)}$ of $g$ determine the quasienergies $\varepsilon_{n}^{(k)}$ of the driven oscillator in the RWA Zhang et al. (2017),

[TABLE]

As seen from this expression, the eigenvalues $g_{n}^{(k)}$ give the the scaled quaisenergies calculated in the reduced Brillouin zone $0\leq\varepsilon<\hbar\omega_{F}/3$. In what follows we call these eigenvalues scaled RWA energies, or just RWA energies, for brevity.

To calculate the eigenvalues and eigenfunctions of the scaled RWA Hamiltonian $g$, it is convenient to rewrite $g$ in terms of the ladder operators $a$ and $a^{\dagger}$,

[TABLE]

Several low-lying eigenvalues of $g$ as functions of the scaled driving amplitude $f$ are shown in Fig. 1(d). It is clear from Eq. (6) and can be also directly seen from Eq. (II.1) that the eigenvalues $g_{n}^{(k)}$ form groups with different values of $k=0,1,2$. The eigenvalues with different $k$ can cross with the varying $f$, as they correspond to different eigenvalues of $\hat{N}_{3}$. The eigenvalues with the same $k$ anticross.

An important qualitative feature of the evolution of the spectrum is that, with the increasing $f$, the eigenvalues with different $k$ cluster into triples. This can be understood from the form of the function $g(Q,P)$ in Fig. 1 (b) and (c). For nonzero $f$ and $\delta\omega>0$, this function has three minima separated by three saddle points. The clustering occurs once the values of $g_{n}^{(k)}$ become smaller than the saddle-point value $g_{s}$ of $g(Q,P)$, which is shown as the red dashed line in Fig. 1. The states with $g_{n}^{(k)}<g_{s}$ correspond to linear combinations of the intrawell states of $g(Q,P)$ that have the same RWA energy in the neglect of interwell tunneling. The tunneling between the three wells leads to the level splitting within a triple of the eigenstates. It sharply falls off as the states go deeper into the wells, see Sec. II.2. As the driving strength increases, the wells of $g(Q,P)$ become deeper whereas the saddle points get closer to the local maximum of $g(Q,P)$ at $Q=P=0$. Thus, the number of intrawell states increases, whereas the number of the states localized near $Q=P=0$ decreases.

II.2 Intrawell states

Of central interest for this paper is the regime of a relatively strong drive where the three wells of the Hamiltonian function $g(Q,P)$ are well-separated in phase space compared to the typical quantum scale of the phase-space area $\lambda$. One can then think of the states localized mostly inside the wells. The number of such states is $\propto 1/\lambda$. We will be interested in the case where this number is large.

We will denote the intrawell states by $|n\rangle_{\nu}\equiv|\psi_{\nu}(Q;n)\rangle$, where $\nu=0,1,2$ enumerates the wells. In what follows we choose the $\nu=0$ well to have a minimum on the $P=0$-axis, see Fig. 1. The wells with $\nu=1,2$ are symmetrically located with respect to the $P=0$ axis, with their minima at the same $Q$. We use $({}Q_{\nu},{}P_{\nu})$ to denote the positions of the minima of the wells. These points are vortices of an equilateral triangle on the $(Q,P)$ plane, see Fig. 1. We denote the minimal value of $g$ as $g_{\min}\equiv g({}Q_{\nu},{}P_{\nu})$.

The wave functions $\psi_{\nu}(Q;n)$ are approximate eigenstates of the Hamiltonian $g$. Near the bottom of a well $\nu$, functions $\psi_{\nu}(Q;n)$ have the form of the eigenfunctions of the harmonic oscillator in their central part. The Hamiltonian of this oscillator is obtained by expanding $g(Q,P)$ about $({}Q_{\nu},{}P_{\nu})$ to the second order in $Q-{}Q_{\nu},P-{}P_{\nu}$. The lowest state was described earlier Zhang et al. (2017). The large-$n$ states can be obtained in the WKB approximation,

[TABLE]

Here, $P^{(\nu)}(Q,g)$ is the solution of the equation $g(Q,P)=g$ that corresponds to the Hamiltonian orbit circling the $\nu$th minimum of $g(Q,P)$ with the classical RWA energy $g$. Parameter $C^{(\nu)}_{n}$ is the normalization constant. For large $n$, one finds $g_{n}$ by calculating $S_{\nu}(Q,g)$ and applying the standard Bohr-Sommerfeld quantization condition Landau and Lifshitz (1977). By symmetry the values of $g_{n}$ are the same for all wells $\nu=0,1,2$.

The intrawell wave functions of different wells overlap, which is the tunneling effect. However, for small $\lambda$ the overlapping is exponentially small. As a result we have, to the leading order in the overlapping,

[TABLE]

The hopping integral $J_{\nu\pm}(g_{n})$ for the lowest intrawell states was found in Ref. Zhang et al., 2017; for higher states it is discussed in Sec. VI.1. The notion of intrawell states is meaningful if $|J_{\nu\pm}(g_{n})|\ll|g_{n}-g_{n\pm 1}|$. Since the symmetry operator $\hat{N}_{3}$ can be thought of as rotating all wells together, $\nu\to\nu+1,\forall\nu$, we see that $J_{\nu\pm}$ is independent of $\nu$ and $J_{\nu+}\equiv J_{0+}=J_{\nu-}^{*}\equiv J_{0-}^{*}$.

To the lowest order in the overlap integrals, the eigenfunctions $\phi_{\nu}^{(k)}$ and eigenvalues $g_{n}^{(k)}$ ($k=0,1,2$) of the RWA Hamiltonian $g$ can be found in the same way as the Bloch wave functions and the energy levels in the tight-binding approximation in solid-state physics Ziman (1979),

[TABLE]

From Eqs. (II.2) - (II.2), the approximately calculated Bloch-type wave functions $\phi_{n}^{(k)}(Q)$ have maxima inside the wells of $g(Q,P)$ and are small outside the wells. It is clear from qualitative arguments that this is also true for the wave functions $\phi_{n}^{(k)}(Q)$ calculated exactly, for example, by numerically diagonalizing the Hamiltonian (II.1). It is important that, since the operator $g(Q,-i\lambda\partial_{Q})$ is real, its eigenfunctions can be also chosen to be real.

Using the exact Bloch-type wave functions $\phi_{n}^{(k)}(Q)$, one can construct Wannier-type wave functions

[TABLE]

Since the exact wave functions $\phi_{n}^{(k)}(Q)$ are orthogonal and normalized, the Wannier functions (11) are also orthogonal and normalized, in contrast to $\psi_{\nu}(Q;n)$. However, one can make functions $\psi_{\nu}(Q;n)$ and $\psi_{\nu;{\rm W}}(Q;n)$ with the same $\nu$ and $n$ very close to each other inside the $\nu$th well, simultaneously for all $\nu$. This requires choosing the exact functions $\phi_{n}^{(k)}(Q)$ to be real and properly choosing the signs of all $\phi_{n}^{(k)}(Q)$ with a given $n$. As Eq. (II.2) suggests, inside the well $\nu=0$, where only $\psi_{0}(Q;n)$ is large (in the absolute value), if we want it to be close to $\psi_{0;{\rm W}}(Q;n)$, all exact wave functions $\phi_{n}^{(k)}(Q)$ should be close to each other. This determines the choice of the otherwise arbitrary signs of the numerically calculated $\phi_{n}^{(k)}(Q)$ with different $k$.

III Dissipation-induced intrawell transitions

Coupling of the oscillator to a thermal reservoir leads to dissipation and to transitions between the intrawell states $|n\rangle_{\nu}\equiv|\psi_{\nu}(Q;n\rangle$ of the Hamiltonian $g(Q,-i\lambda\partial_{Q})$. We consider the parameter range where the dissipation-induced widths of the intrawell levels $g_{n}$, which are determined by the dissipation rate, are much larger than the tunnel splitting $\sim|J_{0\pm}(g_{n})|$. In this case the statistical distribution over the intrawell states of the oscillator is formed much faster then the distribution over different wells. By symmetry arguments, the intrawell transition rates are the same for the different wells of $g(Q,P)$. For concreteness, we will consider transitions between the states in the well $\nu=0$ and use a simplified notation $|n\rangle\equiv|n\rangle_{\nu=0}$.

As mentioned above, we enumerate the intrawell states in such a way that $g_{n}>g_{m}$ if $n>m$. Even at zero temperature, where dissipation corresponds to emission of excitations into the thermal reservoir and thus the oscillator makes transitions only from its Fock states with higher (true) energy to the Fock states with lower energy, in terms of the intrawell states $|n\rangle$ the transitions occur both toward larger and lower $n$. Formally, this is a consequence of the states $|n\rangle$ being linear combinations of the Fock states. The transitions between the states $|n\rangle$ are sketched in Fig. 2.

The rate of a transition $|n\rangle\to|m\rangle$ is determined by the overlapping of the wave functions $|n\rangle,|m\rangle$. Therefore it falls off exponentially with the “distance” $|n-m|$. However, as we show in the next section, because the stationary state populations $\rho_{n}$ also fall off exponentially with increasing $n$, transitions from the remote small-$n$ states close to the bottom of the well can make a dominating contribution to the populations $\rho_{n}$ of large-$n$ states. This means strong non-locality of the distribution $\rho_{n}$ in the state-number space. Such effect is fairly general for systems far from equilibrium, whereas we are not aware of systems in thermal equilibrium where it would be known, unless there are selection rules that make the transition rates very different for different groups of states.

III.1 Balance equation for the state populations

To describe dissipation, we assume that the periodically driven nonlinear oscillator is weakly coupled to a thermal reservoir and that the coupling is linear in the oscillator dynamical variables. Then dissipation comes from the oscillator transitions between neighboring Fock states with emission or absorption of excitations in the reservoir. We further assume that the spectral density of these excitations weighted with the interaction is nearly constant around the oscillator eigenfrequency within a band with width given by the detuning $|\delta\omega|$ and the nonlinear shift of the oscillator eigenfrequency $|\gamma|\langle q^{2}\rangle/\omega_{0}$. Then the oscillator dynamics in slow time $\tau=|\delta\omega|t$ is Markovian and is described by a quantum kinetic equation of the Landau-Lindblad form (cf. Landau (1927); Walls and Milburn (2008)),

[TABLE]

Here, $\kappa$ is the dimensionless decay rate of the oscillator amplitude; it is related to the decay rate of the amplitude in the unscaled time (the coefficient of viscous friction, in the classical description) $\Gamma$ as $\kappa=\Gamma/\|\delta\omega|$; $\bar{n}$ is the Planck number, $\bar{n}=[\exp(\hbar\omega_{0}/k_{B}T)-1]^{-1}$.

Since the overlap integrals of the wave functions in different wells of $g(Q,P)$ are exponentially small, one can disregard mixing of the states in different wells by the relaxation superoperator ${\cal D}$. For $\kappa\gg|J_{0\pm}|$, dissipation-induced intrawell transitions occur much faster than interwell transitions. For the dimensionless time $\tau$ small compared to the long time of the interwell relaxation, the interwell transitions can be disregarded; they are discussed below in Sec. VI.

If the distance between the intrawell levels $g_{n+1}-g_{n}$ is large compared to the level broadening $\propto\lambda\kappa$, the off-diagonal matrix elements of the density matrix $\langle n|\rho|n^{\prime}\rangle$ ($n\neq n^{\prime}$) are small compared to the state populations $\rho_{n}\equiv\langle n|\rho|n\rangle$ (we recall that $|n\rangle$ refers to the intrawell states in the well $\nu=0$; for concreteness, we consider the matrix elements of $\rho$ for the states in this well). If we disregard the off-diagonal matrix elements, we obtain a balance equation for the state populations $\rho_{n}$,

[TABLE]

where $W_{nn^{\prime}}$ are the rates of the dissipation-induced transition $|n\rangle\to|n^{\prime}\rangle$. From Eq. (III.1),

[TABLE]

We note that in a nonequilibrium system the ratio of the transition rates is not given by the Einstein relation for equilibrium systems, $W_{nn^{\prime}}/W_{n^{\prime}n}\neq(\bar{n}+1)/\bar{n}$ Kubo (1957).

III.2 Transition rates in the semiclassical limit

The rates of the dissipation-induced intrawell transitions $W_{n\,n\pm m}$ are determined by the matrix elements

[TABLE]

For small $\lambda$ and for not too large $|m|\ll 1/\lambda$ such matrix elements are given by the Fourier components of the vibrations of the corresponding classical dynamical variable with the scaled RWA energy $g_{n}$ Landau and Lifshitz (1977),

[TABLE]

Here, $\omega(g)$ is the cyclic frequency of the periodic motion along the intrawell classical trajectory $g(Q,P)=g$, which is described by the Hamiltonian equations

[TABLE]

The function $a(\tau;g)$ is calculated as a function of time on this trajectory and is expressed in terms of the coordinate $Q$ and momentum $P$ on the trajectory as $a(\tau;g)=(2\lambda)^{-1/2}[Q(\tau;g)+iP(\tau;g)]$. Equations (14) - (16) allow one to calculate the hopping rates numerically in the WKB approximation.

The WKB ansatz (15) requires $n,n+m\gg 1$. For $n,n+m$ close to 1, the wave functions $|n\rangle,|n+m\rangle$ are close to the wave functions of a harmonic oscillator (see below) and the calculation of the matrix elements of $Q$ and $P$ is straightforward.

In terms of the matrix elements (15), the transition rates become

[TABLE]

Here, we have neglected corrections $\sim\lambda$ and, for large $n$, set $a_{-m}(g_{n+m})\approx a_{-m}(g_{n})$, in which case $W_{n\,n+m}\approx W_{n-m\,n}$.

The matrix elements $a_{m}(g_{n})$ slowly vary with $n$ for large $n$, but rapidly change with $m$. The structure of the solution of the balance equation depends on the asymptotic behavior of $a_{m}$ for large $|m|$, such that $1\ll|m|\ll 1/\lambda$. To find the large-$|m|$-behavior of $a_{m}(g)$, one can shift the integration contour in Eq. (15) into complex time. For positive and negative $m$ the contour is shifted into the lower and upper halfplane of the $\tau$-plane, respectively. The locations of the singularities of the function $a(\tau;g)$ in the complex time plane determine the exponential decay rate of the matrix element $a_{m}$. The details of this calculation are in Appendix C. For large negative $m$, we have

[TABLE]

For large positive $m$, we have

[TABLE]

The parameter $\tau_{\infty}$ is the smallest distance along the imaginary time $\tau$ to the singularity of $Q(\tau;g)$ and $P(\tau;g)$. It is defined by Eq. (56) and is the same in the upper and the lower half-plane of the complex $\tau$ plane. It is shown in Fig. 3.

Since the distance to the singularities is the same in the both half-planes, the matrix elements $|a_{m}(g_{n})|^{2}$, and thus the transition rates, fall off with $|m|$ exponentially with the same exponent proportional to $\tau_{\infty}>0$ for $m>0$ and $m<0$. This is qualitatively different from the rates of transitions between the Floquet states of an oscillator driven by a force with frequency close to the eigenfrequency or a parametrically modulated oscillator at frequency close to twice the oscillator eigenfrequency Marthaler and Dykman (2006); Guo et al. (2013b); Peano and Dykman (2014). In the both latter cases, for $T=0$, the rates of transitions away from the stable state fall off exponentially faster with the interstate distance than the rates of transitions toward the stable state.

We show in Fig. 3 the decrement $\omega(g)\tau_{\infty}(g)$ of the fall-off rate of the matrix elements Eqs. (III.2) and (III.2) with $|m|$. The decrement monotonically decreases with the increasing $g$, and interestingly, for a given $\Delta g{}$, it depends only weakly on the scaled drive amplitude $f$. For $g$ close to $g_{\rm min}$, it diverges logarithmically

[TABLE]

due to the divergence of $\tau_{\infty}(g)$ seen in the inset of Fig. 3, and weakly depends on $f$, for a given $\Delta g$. This is because classical vibrations near the minima of $g(Q,P)$ become almost harmonic. Therefore, higher-order Fourier components of the classical motion decay fast. For $g$ close to the saddle point value $g_{s}$, the period of the classical orbit diverges logarithmically and the motion becomes extremely anharmonic: very slow near the saddle points, but fast away from them. As a result, $\omega(g)\tau_{\infty}(g)$ approaches zero logarithmically for $g$ close to $g_{s}$.

Another distinctive and important feature of period tripling is that the singularities of $a(\tau;g)$ in the complex plane $\tau$ are branching points rather than poles, as for resonant driving or parametric modulation. The branching points at $i\tau_{\infty}$ and $-i\tau_{\infty}$ are different. Respectively, the prefactors in Eqs. (III.2) and (III.2) show a power-law dependence on $m$, but with different exponents.

The parameter used to obtain the large-$|m|$ asymptotic behavior of the coefficients $a_{m}(g)$ is $|m|\omega(g)\tau_{\infty}\gg 1$. For $f\lesssim 1$ we have $\tau_{\infty}\sim 1$, see Fig. 3 ($\tau_{\infty}\gg 1$ for $g\to g_{\min}$). Then Eqs. (III.2) and (III.2) apply provided $|m|\omega(g)\gg 1$. In this case for $\bar{n}=0$ the ratio of the rates of transitions away and toward $g_{\min}$ is $W_{n\,n+m}/W_{n\,n-m}\propto[f^{2}/m\omega(g)]^{2/3}\ll 1$ ($m>0$). This shows that, in the stationary regime, the system is mostly localized near the minimum of the well of $g(Q,P)$: it is much more probable to make a transition toward this minimum than away from it. From Eq. (C), for $f\gg 1$ the frequency $\omega(g)$ scales with $f$ as $f^{2}$, and therefore still $W_{n\,n+m}/W_{n\,n-m}\ll 1$ for $m>0$ as long as $m\omega(g)\tau_{\infty}(g)\gg 1$.

III.3 Breakdown of the detailed balance

A periodically driven oscillator is far from thermal equilibrium. Generally, it should not have detailed balance. The absence of detailed balance was shown for a resonantly driven oscillator in the classical limit Dykman and Krivoglaz (1979). However, as mentioned in the Introduction, in the deeply quantum regime of zero temperature, a resonantly driven Duffing oscillator with the relaxation described by the master equation (III.1) has detailed balance in the RWA Drummond and Walls (1980). Moreover, for the same type of relaxation, detailed balance for zero temperature was found also for a Duffing oscillator parametrically modulated close to twice its eigenfrequency Kryuchkyan and Kheruntsyan (1996). The physical reason for the detailed balance, which exists only for $T=0$ and only for this relaxation mechanism is not known (it breaks down in the presence of dephasing Marthaler and Dykman (2006)). It is interesting and, in a way, illuminating, to find out whether detailed balance persists for $T=0$ and relaxation described by Eq. (III.1) in the case of an oscillator driven close to triple the eigenfrequency.

The condition of detailed balance is the relation between the rates of interstate transitions that shows that the ratio of the rates of direct transitions back and forth between the states is equal to the ratio of the rates of transitions via an intermediate state,

[TABLE]

For $T=0$ we have $W_{nn^{\prime}}=2\kappa|a_{n^{\prime}-n}(g_{n})|^{2}$. As seen from Eqs. (III.2) and (III.2), for period tripling there is no detailed balance even for $T=0$. This shows that the detailed balance for modulation at frequencies $\approx\omega_{0}$ and $\approx 2\omega_{0}$ for $T=0$ is non-generic.

IV Quasistationary distribution over intrawell states

Over the dimensionless relaxation time $\sim\kappa^{-1}$, transitions between the intrawell states form a quasistationary distribution over these states. This distribution can be found from the stationary solution of the balance equation (13). An important property of the transition rates is

[TABLE]

that is, the system is more likely to go from a higher-lying to a lower-lying state than in the opposite direction. This property is seen from Eqs. (III.2) and (III.2) for large $m$ and from the numerical evaluation of the semiclassical rates $W_{nn^{\prime}}$ based on Eq. (15) for $|n-n^{\prime}|\sim 1$. It is clear then that the stationary distribution will be maximal near the bottom of the well of $g(Q,P)$. However, since the probabilities of transitions up along the $g_{n}$ axis are nonzero, the stationary distribution over the levels $g_{n}$ will have a nonzero width.

We are interested in the semiclassical parameter range $\lambda\ll 1$, where the number of intrawell bound states is large. The analysis of the distribution should be done differently in different ranges of the values of the scaled RWA energy $g$. We note that it does not apply in a narrow range of $g$ close to the saddle value $g_{s}$, where the interwell transitions have to be taken into account.

IV.1 The vicinity of the minima of the wells of $g(Q,P)$

The distribution over intrawell states can be found analytically near the minima of the RWA Hamiltonian function $g(Q,P)$. In the rotating frame, the vibrations of $Q,P$ about a minimum are almost harmonic and, respectively, the levels $g_{n}$ are almost equidistant. To find the distribution for the well $\nu=0$, we expand $g(Q,P)$ about the $\nu=0$-minimum,

[TABLE]

Here $(Q_{0},P_{0}=0)$ is the position of the minimum and $g_{\min}$ is the value of $g(Q,P)$ at the minimum, $g_{\min}=fQ_{0}^{2}(3f-4Q_{0})/12$.

It is convenient to change from $P$ and $Q-Q_{0}$ to the raising and lowering operator $b^{\dagger}$ and $b$ using a squeezing transformation,

[TABLE]

where $\omega_{\rm min}$ is the vibration frequency at the bottom of the well of $g(Q,P)$ and the squeezing parameter $\phi_{*}$ is given by the equation $\tanh\phi_{*}=(|g_{QQ}|^{1/2}-|g_{PP}|^{1/2})/(|g_{QQ}|^{1/2}+|g_{PP}|^{1/2})$. The intrawell states $|n\rangle$ near the bottom of the well are well approximated by the eigenstates of the operator $b^{\dagger}b$.

From Eq. (IV.1), the lowering operator of the oscillator $a=(2\lambda)^{-1/2}(Q+iP)$ is a linear combination of the operators $b$ and $b^{\dagger}$. Therefore the transitions rates $W_{n\,n+m}$, Eq. (14) are nonzero for $m=\pm 1$, and this is the case even for $\bar{n}=0$, where the oscillator makes dissipative transitions only to a lower Fock state. The rates $W_{n\,n+m}$ are easy to find from Eq. (IV.1), taking into account that $\langle n|b|n+m\rangle=(n+1)^{1/2}\delta_{m,1}$. Substituting the rates into the balance equation (13), we find that the stationary solution of this equation near the minimum of the well is

[TABLE]

where $n_{\rm eff}$ is the effective Planck number,

[TABLE]

Equation (IV.1) has the form of the Boltzmann distribution. As seen from Eq. (24), the effective temperature of the distribution over the intrawell states is nonzero even where the effective temperature of the thermal reservoir is equal to zero. This is the quantum heating effect Dykman and Smelyanskii (1988); Marthaler and Dykman (2006). For a resonantly driven oscillator it was observed experimentally in Ref. Ong et al., 2013.

The dependence of $\bar{n}_{\rm eff}$ on $f$ is different for different signs of the detuning of the drive frequency $\delta\omega=(\omega_{F}/3)-\omega_{0}$. For $\delta\omega<0$, as mentioned previously, $g(Q,P)$ has a minimum at $Q=P=0$, whereas the wells at nonzero $Q^{2}+P^{2}$ appear only for $f>2$. In this case $\bar{n}_{\rm eff}$ monotonically decreases as $f$ increases. Near the threshold $f=2$, $\bar{n}_{\rm eff}$ diverges as $(f-2)^{-1/4}$. For large $f\gg 1$, $\bar{n}_{\rm eff}$ asymptotically approaches $(2\bar{n}+1)/\sqrt{3}-1/2$.

For $\delta\omega>0$ function $g(Q,P)$ has the form of a Mexican hat, with the top at $Q=P=0$ and the rim at $Q^{2}+P^{2}=1$ for $f=0$. The wells $\nu=0,1,2$ emerge for $f>0$ with no threshold, as a modulation of the depth of the rim. The effective Planck number $\bar{n}_{\rm eff}$ is non-monotonic as a function of the scaled drive amplitude $f$. It diverges as $f^{-1/2}$ for $f\to 0$. As $f$ increases, $\bar{n}_{\rm eff}$ reaches a minimum at $f=1/\sqrt{2}$, where $\bar{n}_{\rm eff}=\bar{n}$. This corresponds to the disappearance of the squeezing, $\sinh\phi_{*}=0$ for $f=1/\sqrt{2}$. For $f>1/\sqrt{2}$, $\bar{n}_{\rm eff}$ increases with the increasing $f$ and approaches the same asymptotic value as in the case $\delta\omega<0$.

We show in Fig. 4 the scaled effective inverse temperature $\lambda\omega_{\rm min}/T_{\rm eff}$ as a function of the scaled drive amplitude $f$ for $\delta\omega>0$. For $\bar{n}=0$, the inverse temperature diverges at $f=1/\sqrt{2}$ as $\left|\log(f-1/\sqrt{2})^{2}\right|$.

For $f=1/\sqrt{2}$ and $\bar{n}=0$ the analysis of the distribution over the states $|n\rangle$ has to go beyond the harmonic approximation (IV.1) and take into account the terms of higher order in $P$ and $Q-Q_{0}$ in Eq. (IV.1). By perturbation theory in the parameter $\lambda$, the cubic terms in $P$ and $Q-Q_{0}$ lead to a finite transition rate $W_{02}\propto\lambda$, and more generally, $W_{n\,n+2}\propto\lambda$ for $n\sim 1$. The finite rates of the upward in $g$ transitions result in a nonzero population of the excited states; the distribution over the states falls off as a power series in $\lambda$.

IV.2 The eikonal approximation

To find the stationary solution of the balance equation (13) away from the minimum of $g(Q,P)$ it is convenient to seek this solution in the eikonal form,

[TABLE]

As it stands, Eq. (25) just uses the fact that $\rho_{n}\geq 0$. It applies near the minimum of the well of $g(Q,P)$ as well, as seen from Eq. (IV.1), but is not limited to this range of $g$.

We will assume, and then check a posteriori, that $R(g)$ is a smooth function of $g$. This is an analog of the eikonal approximation in optics or the WKB approximation in quantum mechanics. In this approximation, for $|n^{\prime}-n|\ll 1/\lambda$ we have

[TABLE]

where $R^{\prime}(g)\equiv dR/dg$; we have used here the semiclassical expression for the quasienergy level spacing in terms of the classical vibration frequency $\omega(g)$.

The approximation (IV.2) was used to find the stationary distribution of resonantly and parametrically driven quantum oscillators Dykman and Smelyanskii (1988); Marthaler and Dykman (2006) from the balance equation for the state populations (13). As we show, it allows one to solve the balance equation for the considered here period tripling problem in a broad parameter range, but it can break down for $\bar{n}=0$. This is qualitatively different from the cases considered earlier.

We emphasize that Eq. (IV.2) assumes that the logarithm of the distribution $\rho_{n}$ smoothly depends on $n$, whereas the distribution itself does not have to be a smooth function of $n$ and can vary by a factor ${\cal O}(1)$ when $n$ is incremented by 1. This is similar to the WKB approximation, where the wave function is oscillating fast, but the wavelength changes only a little on the distance given by the wavelength itself.

From Eqs. (25) and (IV.2), we see that $\rho_{n+m}/\rho_{n}=\exp[-m\omega(g_{n})R^{\prime}(g_{n})]$ for $|m|\ll 1/\lambda$. This means that $\omega(g_{n})R^{\prime}(g_{n})$ can be interpreted as the $g_{n}$–dependent inverse temperature of the Boltzmann-like distribution over the RWA energy levels $g_{n}$.

Equations (13) and (IV.2) reduce the problem of finding the distribution over the excited intrawell states ($n\gg 1$) to an algebraic equation

[TABLE]

Here the hopping rates $W_{nn^{\prime}}$ are given by the semiclassical expression (17). As indicated above, for $|n-n^{\prime}|\sim 1$ they can be calculated numerically from Eqs. (16), whereas for large $|n-n^{\prime}|$ we have the explicit expressions (17) - (III.2). In writing Eq. (27) we used $W_{n\,n+m}=W_{n-m\,n}$ for large $n\gg 1$ and $|m|\ll n$.

A crucial feature of Eq. (27) is that it is local in the level number $n$. It assumes that the stationary population $\rho_{n}$ is formed by the transitions from a few states $n+m$ surrounding a given state $n$, with $|m|\ll n,1/\lambda$. Formally, it requires that the sum over $m$ in Eq. (27) converges. If this is the case, Eq. (27) defines a smooth function $R^{\prime}(g)$ which depends on the single parameter, the scaled amplitude of the driving field $f$. However, the applicability of this approximation is not known a priori. We note that Eq. (27) does not contain the scaled Planck constant $\lambda$; in particular, $R^{\prime}$ is independent of $\lambda$.

IV.2.1 The vicinity of the minima of the Hamiltonian function

Typtically, the WKB approximation applies for highly excited states that correspond to large $n$ in Eq. (27). In the considered case the situation is different. Because near the minima of the quasienergy surface the dynamics in the rotating frame can be mapped onto the dynamics of an auxiliary weakly nonlinear oscillator, cf. Eq. (IV.1), the semiclassical approximation works for $g$ approaching $g_{\rm min}$. One can show that for $g-g_{\min}\ll 1$, the semiclassical transition rates $W_{n+m\,n}$ in Eq. (17) scale as $W_{n+m\,n}\propto(g-g_{\rm min})^{|m|}$. If one sets to zero the rates $W_{n+m\,n}$ with $|m|\geq 2$, Eq. (27) gives the same result for the distribution as in Eq. (IV.1), with $R^{\prime}(g_{\rm min})\omega(g_{\rm min})\to\lambda\omega_{\rm min}/T_{\rm eff}$ for $g\to g_{\min}$. Clearly, the relation $W_{n\,n+m}=W_{n-m\,n}$ does not apply for small $n\sim 1$, but the range of small $n$ corresponds to a very narrow range of $g$ for small $\lambda$.

In Fig. 4 we show the result of solving Eq. (27) for $g$ close to $g_{\rm min}$. It matches the solution (IV.1) obtained in the harmonic approximation well in a broad parameter range. However, it also shows that the transitions with $|m|>1$ lead to an occupation of the excited intrawell states for $\bar{n}=0$ even for the value of the driving amplitude $f$ where, in the harmonic approximation, such states remain unoccupied. This shows that $T_{\rm eff}$ becomes nonzero if the transitions with $|m|>1$ are taken into account.

IV.2.2 Classical limit

For large Planck numbers $\bar{n}$, Eq. (27) greatly simplifies, since $\omega(g)R^{\prime}(g)\ll 1$. One can then expand $\exp[-m\omega(g)R^{\prime}]$ in a series to second order in $R^{\prime}$. Then, from Eqs. (17) and (27),

[TABLE]

Using that the matrix elements $a_{m}(g)$ are the Fourier components of the dynamical variable $a(\tau;g)\propto Q(\tau;g)+iP(\tau;g)$ calculated for a classical orbit (16) with a given $g$, the sums in Eq. (28) can be written as integrals over the area enclosed by this orbit in the phase space Goldstein et al. (2001). Then Eq. (28) can be rewritten as,

[TABLE]

As indicated above, the integrals here are taken over the interior of the region $g(Q,P)=g$ within a well of $g(Q,P)$ and ${\bm{\nabla}}^{2}g(Q,P)\equiv\partial^{2}_{Q}g+\partial^{2}_{P}g=4(P^{2}+Q^{2})-2\,{\rm sgn}(\delta\omega)$.

In the high temperature range $k_{B}T\gg\hbar\omega_{0}$, $R^{\prime}$ is inversely proportional to $T$. The distribution over intrawell states goes over into the classical (but non-Boltzmann) distribution formed as a result of thermal-noise-induced diffusion over the RWA-energy.

Equation (28) applies for an arbitrary temperature near the threshold of the period tripling, but not too close to the threshold, so that the condition $\omega(g_{\min})\gg\kappa$ still holds. Being near the threshold means $f\ll 1$ for $\delta\omega>0$ and $f-2\ll 1$ for $\delta\omega<0$, to zeroth order in $\kappa$. Near the threshold, the vibration frequency $\omega(g)\ll 1$ for the whole range of the intrawell values of $g$. Indeed, near the minimum of $g(Q,P)$ we have $\omega_{\rm min}\approx(6f)^{1/2}\ll 1$ for $\delta\omega>0$ and $\omega_{\rm min}\approx 2\sqrt{3}(f-2)^{1/4}\ll 1$ for $\delta\omega<0$. For larger $g$, $\omega(g)$ monotonically decreases and approaches zero as $g$ approaches $g_{s}$.

We compare in Appendix D the result of the classical limit given by Eq. (28) with the semiclassical calculation based on Eq. (27). The classical results match well the semiclassical ones already for moderately small $\bar{n}\sim 1$.

V The breakdown of the local approximation

The eikonal approximation relies on the stationary distribution being formed locally by transitions between a few nearest states. The number of the states involved should be much smaller than $1/\lambda$, otherwise the expansion of the logarithm of the probability distribution (IV.2) does not apply. As a consequence, the polynomial equation for $R^{\prime}$ (27), which is a direct analog of the equation for the derivative of the logarithm of the wave function in the WKB approximation, does not apply either. It turns out that, for $\bar{n}=0$, the locality, and thus the eikonal approximation, can break down for the oscillator that displays period tripling. In this section, we study the condition of such a breakdown and show how to go beyond the local approximation. We focus on the case of a positive detuning $\delta\omega>0$.

V.1 The locality condition and its breakdown for zero temperature

The stationarity of the intrawell state populations $\rho_{n}$ means that the probability flux out of a state $n$ (the rate of leaving the state) and the probability flux into the state are equal. We can write the balance equation (13) in terms of the relative incoming and outgoing fluxes as

[TABLE]

The exponential decay of the rate $W_{n\,n+m}\propto\exp[-2|m|\omega(g_{n})\tau_{\infty}(g_{n})]$ for large $|m|$ guarantees that the total relative outgoing flux $\sum_{m}{\cal J}^{\rm out}_{m}$ is always finite. Since the stationary state populations $\rho_{n}$ decrease with the increasing $n$, that is, $\rho_{n+m}/\rho_{n}<1$ for $m>0$, the relative incoming flux from the states with larger RWA energy $\sum_{m>0}{\cal J}^{\rm in}_{m}(g_{n})$ is also finite. Moreover, it is local, only a few nearest states contribute to this flux.

In contrast, but for the same reason, the relative flux ${\cal J}^{\rm in}_{m<0}$ into a given state from the states with lower RWA energy can be non-local. If the increase of the population $\rho_{n-|m|}$ with the increasing $|m|$ is faster than the decrease of the transition rate $W_{n-|m|\,n}$, the incoming relative flux ${\cal J}^{\rm in}_{m<0}$ increases with the increasing distance between the states $|m|$.

If we use the eikonal approximation for the state populations, Eq. (IV.2), we find that, for large $|m|$ (but $|m|\ll 1/\lambda$) and for a given $n$, ${\cal J}^{\rm in}_{m<0}\propto\exp[-|m|\omega(2\tau_{\infty}-R^{\prime})]$. Thus, the locality condition for a given $n$ is

[TABLE]

The right-hand side of this inequality gives the fall-off rate of the relative influx with the distance to the lower-lying state, in the eikonal approximation. The left-hand side is proportional to the reciprocal number of states in the well of $g(Q,P)$, and thus puts a limit on the number of the lower-lying states. When the condition (31) is not met, the eikonal approximation (IV.2) does not allow one to find the stationary probability distribution over the intrawell states. One needs to use the full balance equation (13) (we remind that, in fact, we are discussing the quasi-stationary distribution, as we have so far disregarded interwell transitions). Also, one should not use Eq. (17) to calculate the transition rates, but rather the full WKB wavefunctions or the Wannier-type functions to calculate the relevant matrix elements.

We start with analyzing the locality breakdown in the semiclassical limit $\lambda\rightarrow 0$. The function $\omega(g)[2\tau_{\infty}(g)-R^{\prime}(g)]$ calculated in the “local” approximation (27) is shown in Fig. 5. For relatively large (purple line) or small (red line) values of the scaled driving field $f$, function $2\tau_{\infty}(g)-R^{\prime}(g)$ remains positive for all intrawell values of $g$, from $g=g_{\rm min}$ to $g=g_{s}$. Thus, the locality holds. However, for $f\sim 1$, function $2\tau_{\infty}(g)-R^{\prime}(g)$ becomes equal to zero at some $g=g_{\rm NL}$ as seen from the data shown by the blue and green lines in Fig. 5. For $g>g_{\rm NL}$, Eq. (27) does not have a nontrivial solution for $R^{\prime}$, signaling the breakdown of the locality.

We denote the minimum value of $g$ where the locality breaks down by $g_{\rm NL}$. In Fig. 6 we show the dependence of $g_{\rm NL}$ on the scaled driving amplitude $f$. As seen from the figure, $g_{\rm NL}$ first emerges with the increasing $f$ at the saddle point value $g_{s}$ of $g(Q,P)$ for $f\approx 0.2$. It then moves deeper into the well of $g(Q,P)$. The distance $\Delta g_{\rm NL}$ from $g_{\rm NL}$ to the bottom of the well scaled by the depth of the well $g_{s}-g_{\min}$ is minimal for $f\approx 0.5$. As $f$ further increases, $g_{\rm NL}$ goes back to $g_{s}$. For $f\gtrsim 1.4$ the nonlocality disappears.

To find accurately the location of $g_{\rm NL}$ from Eq. (27), it is important to take into account the large-$|m|$ contribution to the influx ${\cal J}^{\rm in}_{m<0}$. This is because for $g$ close to $g_{\rm NL}$, ${\cal J}^{\rm in}_{m<0}(g)$ falls off with $|m|$ rather slowly, with the fall-off exponent $\omega(g)[2\tau_{\infty}(g)-R^{\prime}(g)]$ close to zero as can be seen in Fig. 5. The large-$|m|$ contribution can be calculated using the asymptotic expression for the transition rates, Eqs. (14) and (III.2).

In the region of nonlocality we have $R^{\prime}(g)>2\tau_{\infty}(g)$. The stationary probability distribution in this region is found numerically and is discussed in Sec. V.3.

V.2 Effect of nonzero temperature

The very existence of the non-locality in the transitions among intrawell states depends sensitively on the interrelation between the temperature and the scaled Planck constant $\lambda$. In the limit $\lambda\rightarrow 0$ the locality is formally preserved at any nonzero temperature. This is seen from the expression (17) for the transition rates $W_{nn^{\prime}}$. At nonzero temperatures, the second term in $W_{nn^{\prime}}$, which comes from the absorption of energy from the thermal reservoir, becomes non-zero. Although this term is small for low temperatures, it can significantly change the distribution. We now provide analytical arguments to show that this is the case. The numerical results are given in the next section.

The relative influx into the state with a given RWA energy $g_{n}$ can be written as a sum of the contribution from nearby states [small $|m|$ in Eq. (V.1)] and that from the states that are further away (large $|m|$). We will use that, where the locality holds, $\rho_{n+m}/\rho_{n}=\exp[-mR^{\prime}(g_{n})\omega(g_{n})]$ and consider the “most dangerous” situation where $R^{\prime}\approx 2\tau_{\infty}$ and therefore the contribution of the remote states is significant. Using Eqs. (III.2) and (III.2) we find

[TABLE]

Here $C_{1},C_{2}$ are constants independent of $R^{\prime}$, whereas the parameter $C_{0}$ comes from the small-$|m|$ terms in Eq. (V.1) and depends on $2\tau_{\infty}-R^{\prime}$ in a nonsingular way.

As seen from Eq. (V.2), for $R^{\prime}(g)\to 2\tau_{\infty}(g)$ the part of the relative influx, which is proportional to $\bar{n}$, diverges. The other terms in Eq. (V.1) remain finite. Therefore, for any nonzero $\bar{n}$, the“local” balance equation (27) has a nontrivial stationary solution for $R^{\prime}(g)$ with $R^{\prime}(g)<2\tau_{\infty}(g)$, as the term $\propto\bar{n}$ can compensate other terms. In other words, the locality is preserved in the limit $\lambda\rightarrow 0$ but $\bar{n}\neq 0$. However, as we have seen, the locality can break down for $\bar{n}=0$.

Where $R^{\prime}(g_{n})$ approaches $2\tau_{\infty}(g_{n})$, the number of states that contribute to the influx into the state $n$ increases $\propto[2\tau_{\infty}(g_{n})-R^{\prime}(g_{n})]^{-1}$. For $\bar{n}\to 0$ this number diverges. Physically, the total number of states in a well of $g(Q,P)$ is always finite. It is $\propto 1/\lambda$. Therefore for a nonzero $\lambda$, the transition from the non-locality at $\bar{n}=0$ to the locality at $\bar{n}\neq 0$ is smeared out. The “non-locality” means that essentially all states with $g_{m}<g_{n}$ contribute to the stationary population $\rho_{n}$, whereas the “locality” means that the number of such states is much smaller. For a small nonzero $\bar{n}$, the number of the contributing states is given by the condition that the term $\propto\bar{n}[2\tau_{\infty}(g)-R^{\prime}(g)]^{-1/3}$ in Eq. (V.2) is of the same order as other terms in the balance equation. The number of the contributing states $\sim[2\tau_{\infty}(g)-R^{\prime}(g)]^{-1}$ is thus $\propto\bar{n}^{-3}$. Therefore the locality holds where $\bar{n}^{3}\gtrsim\lambda$.

We note that there is a profound difference between the non-locality for period tripling discussed here and that for an oscillator driven near its eigenfrequency Guo et al. (2013b); Peano and Dykman (2014) or parametrically modulated near its second overtone Marthaler and Dykman (2006). There, the locality holds for $\bar{n}=0$. The non-locality emerges only in a narrow range of nonzero $\bar{n}$ due to the asymmetry of the decay of the semiclassical matrix elements $a_{m>0}$ and $a_{m<0}$ with $|m|$. In contrast, in the period-tripling case, the semiclassical matrix elements $a_{m>0}$ and $a_{m<0}$ decay with $|m|$ with the same exponents, albeit with different prefactors. The non-locality then occurs at zero temperature. Another qualitative distinction is that, for period tripling, the logarithm of the probability distribution $\log\rho_{n}=-R(g_{n})/\lambda$ remains a smooth function of $g_{n}$. The function $R(g)$ can still be expanded in a series, it is just that keeping the first term in this expansion is not sufficient for finding the probability distribution where the locality breaks down.

V.3 Numerical analysis using the Wannier functions

To confirm the breakdown of the local approximation, we numerically find the Wannier-type intrawell wave functions defined by Eqs. (II.2) and (11), then calculate the transition rates using the matrix elements of the operators $a,a^{\dagger}$ on these wave functions, and then solve the balance equation (13). This method does not require the local approximation, but it significantly relies on the scaled Planck constant $\lambda$ being nonzero. An alternative approach is to calculate the matrix elements using the Landau method Landau and Lifshitz (1977) for states with significantly different quantum numbers. The Landau method was used Guo et al. (2013b) in the problem of an oscillator driven close to its eigenfrequency, where the singularities of the classical trajectories (16) were simple poles periodically located on the complex-time plane. In the present case the trajectories have multiple branching points, which significantly complicates the application of the method.

We compare in Fig. 7 the inverse temperature $R^{\prime}(g)$ obtained using the Wannier-type functions with that obtained from Eq. (27). For a broad range of parameters, the inverse temperatures $R^{\prime}(g_{n})$ obtained using the two methods coincide. When $2\tau_{\infty}$ becomes close to $R^{\prime}(g_{n})$, the result obtained from Eq. (27) deviates from the Wannier-functions method, indicating the breakdown of the local approximation that underlies Eq. (27).

To visualize the nonlocality, we show in Fig. 8 the normalized probability flux

[TABLE]

from the state $|n^{\prime}\rangle$ to the state $|n\rangle$ inside a well of $g(Q,P)$ in the stationary regime. The normalization factor is the maximal over $n^{\prime}$ value of the flux $\rho_{n^{\prime}}W_{n^{\prime}n}$. The flux is obtained from Eq. (13) using the transition rates calculated with the Wannier-type wave functions (11). At zero temperature (left panels), the fall-off of the flux (33) with the increasing interstate distance $n-n^{\prime}$ becomes strongly non-exponential as the RWA energy of the state into which the transition occurs $g_{n}$ approaches $g_{\rm NL}$. For relatively large values of $g_{n}$, the flux $\rho_{n^{\prime}}W_{n^{\prime}n}$ becomes nonmonotonic in $n-n^{\prime}$. This is a manifestation of the nonlocality: a significant portion of the probability flux into a high-lying state $|n\rangle$ comes from the states much deeper in the well, because the increase of their population with $n-n^{\prime}$ is faster than the fall-off of the transition rate $W_{n^{\prime}n}$. Moreover, the flux from the states that are further away can be larger than from at least some states closer to a given state $n$. At finite temperatures (right panels), such behavior goes away and $\rho_{n^{\prime}}W_{n^{\prime}n}$ falls off exponentially as $|n-n^{\prime}|$ increases.

VI Escape from the period-three states of a dissipative oscillator

In the discussions in Secs. IV and V, we have disregarded tunneling between the states localized in different wells of $g(Q,P)$. As explained in Sec. II.2, the tunneling causes splitting of the intrawell levels into triplets with the level spacing $\sim|J_{0\pm}|$, see Eq. (II.2). If the scaled relaxation rate $\kappa$ is smaller than $|J_{0\pm}|/\lambda$, dissipation leads to transitions within the triplets and between different triplets, and ultimately there is formed a stationary distribution over the broadened tunnel-split Floquet states of the oscillator. However, for small $\lambda$ the level splitting $\sim|J_{0\pm}(g_{n})|$ is exponentially small. Therefore such picture is relevant provided $\kappa$ is also exponentially small.

In this section we focus on the parameter regime where the tunnel splitting is small compared to the dissipative level broadening, $|J_{0\pm}|\ll\lambda\kappa$. In this case coherent tunneling is replaced by an interwell hopping (switching) with a rate $\propto|J_{0\pm}|^{2}/\lambda^{2}\kappa$, an analog of the quantum diffusion in solids Kagan and Leggett (1992). Since the intrawell states in the rotating frame correspond to the vibrations at frequency $\omega_{F}/3$ with different phases, such interwell switching corresponds to phase-flip transitions in the laboratory frame.

If one disregarded the population of the excited intrawell states, the switching would occur in the lowest state. The dissipation-induced transitions between the intrawell states significantly complicate the picture. The distribution over the intrawell states is formed fast, over the dimensionless time $\kappa^{-1}\ll\lambda^{2}\kappa/|J_{0\pm}|^{2}$. The populations of the excited intrawell states $\rho_{n}$ fall off exponentially with the increasing $g_{n}$, whereas the tunneling matrix elements $|J_{0\pm}(g_{n})|$ exponentially increase with $g_{n}$. The rate of the interwell hopping is then determined by some optimal $g_{n}$ where $\rho_{n}|J_{0\pm}(g_{n})|^{2}$ has a sharp maximum. Such competition of the exponential factors is similar to that in systems in thermal equilibrium, where the intrawell distribution is of the Boltzmann form Larkin and Ovchinnikov (1985).

As we show below, in our case the optimal value of the RWA energy for interwell switching corresponds to the height of the barrier between the wells of $g(Q,P)$, i.e., the saddle-point value $g_{s}$. In other words, escape from a well of $g(Q,P)$ occurs via an over-barrier transition. This is an analog of thermal activation in systems in thermal equilibrium and is called quantum activation, since it is due to quantum fluctuations and occurs even for $T=0$. For quantum oscillators driven close to their eigenfrequency or parametrically modulated close to twice the eigenfrequency, quantum activation was found earlier Dykman and Smelyanskii (1988); Marthaler and Dykman (2006). These systems have detailed balance for $T=0$ Drummond and Walls (1980); Kryuchkyan and Kheruntsyan (1996), whereas our system does not, and therefore the nature of the escape can be expected to be different.

We will focus on the case of positive detuning, $\delta\omega>0.$ For negative detuning, $\delta\omega<0$, along with interwell tunneling there also occurs tunneling between the intrawell states and the state localized near $Q=P=0$.

VI.1 The condition for the onset of quantum activation

For small $\lambda$ the tunneling matrix element $J_{0\pm}(g_{n})$ can be calculated in the WKB approximation. The problem is somewhat different from the problem of tunneling of a particle in a potential well, where the momentum as a function of coordinate and energy has only two values and is either real or imaginary. In the present case, the classical momentum $P(Q,g)$ defined by the equation $g(Q,P)=g$ has 4 branches as a function of the coordinate and is complex rather than purely imaginary in the classically inaccessible region. As a result the WKB wave function is not just decaying, but also oscillating in this region. Nevertheless, extending the analysis for the lowest intrawell state Zhang et al. (2017), one can show that

[TABLE]

Here, $S_{\rm tun}(g)$ is the imaginary part of the complex classical action calculated along the Hamiltonian trajectory (16) that goes in complex time from one well of $g(Q,P)$ to another well with the same RWA energy $g$. One has to choose the trajectory with the minimal $S_{\rm tun}$.

To find the optimal value of the RWA energy $g$ for interwell hopping, we will use a quasicontinuous approximation for the intrawell distribution and replace $\rho_{n}=\exp[-R(g_{n})/\lambda]$ with $\rho(g)=\exp[-R(g)/\lambda]$. In this approximation the rate of interwell hopping with the RWA energy between $g$ and $g+dg$ is $\propto dg\,\exp\{-[R(g)+2S_{\rm tun}(g)]/\lambda\}$. The optimal $g$ is determined by the extremum of the exponent.

If escape from a well of $g(Q,P)$ occurs via quantum activation, it means that the amplitude of the exponent decreases with the increasing $g$ for all values of $g$ inside the well of $g(Q,P)$,

[TABLE]

We found in Sec. V that, where the intrawell distribution is formed by transitions between comparatively close states (the locality condition), $R^{\prime}(g_{n})$ is limited by the decrement $2\tau_{\infty}(g_{n})$ of the decay of the intrawell transition rates $W_{nn^{\prime}}$ with $|n-n^{\prime}|$, i.e., $0<R^{\prime}(g)\leq 2\tau_{\infty}(g)$. This limits the first term in Eq. (35).

The second term in Eq. (35) is negative and is determined by the imaginary part of the time of moving along a Hamiltonian trajectory (16) through the classically inaccessible region between the wells,

[TABLE]

The derivative $\partial_{P}g$ is calculated here for $P=P(Q,g)$. The details of the calculation are given in Appendix E. We note that the calculation is done differently in different regions of $g$ where the trajectories in the classically allowed region are ellipse-like or horse-shoe-like.

We show in Fig. 9 the imaginary tunneling time $\tau_{\rm tun}$ for various values of the driving amplitude. For $g$ close to $g_{\rm min}$, $\tau_{\rm tun}$ diverges logarithmically,

[TABLE]

For $g$ close to $g_{s}$, one can show that $\tau_{\rm tun}$ is linear in $g_{s}-g$ and is finite at $g=g_{s}$; see Appendix E.

Also shown in Fig. 9 is the function $\tau_{\infty}(g)$, which limits $R^{\prime}(g)/2$ in the region of locality. As seen from Fig. 9, $|\tau_{\rm tun}|>\tau_{\infty}$, thus guaranteeing that the quantum-activation condition (35) holds true. Note that for $g$ close to $g_{\rm min}$, both $\tau_{\infty}$ and $\tau_{\rm tun}$ diverge logarithmically, but the prefactor of $\tau_{\rm tun}$ is twice as large as that of $\tau_{\infty}$ as can be seen from Eqs. (20) and (36).

In the range of the driving amplitude and the RWA energy where, for $\bar{n}=0$, the intrawell probability distribution is formed by transitions from remote states, we have $R^{\prime}(g)>2\tau_{\infty}(g)$. We have found that the numerically calculated value of $R^{\prime}$ in this range satisfies the condition (35). Therefore escape from a well of $g(Q,P)$ does occur via quantum activation.

VI.2 Quantum activation energy

The rate of escaping from a well of $g(Q,P)$ in the rotating frame is the rate of escaping from a period-3 vibrational state in the laboratory frame. To logarithmic accuracy this rate is

[TABLE]

The parameter $R_{A}$ is the quantum activation energy of escape. The quantum nature of the fluctuations leading to escape is clear from the fact that the escape rate displays activation dependence on $\hbar$, $\log W_{\rm esc}\propto\lambda^{-1}\propto\hbar^{-1}$. Yet, as we mentioned, the escape occurs via an overbarrier transition, as if it were thermally activated, with $\hbar$ playing the role of temperature.

We show in Fig. 10 the dependence of the activation energy $R_{A}$ on the scaled driving amplitude $f$. In contrast to the inverse effective temperature of the intrawell distribution near the minima of $g(Q,P)$, which is shown in Fig. 4, $R_{A}$ monotonically increases as $f$ increases. We present the results for $\bar{n}>\lambda^{1/3}$. In this case the intrawell distribution is formed by transitions $|n\rangle\to|n^{\prime}\rangle$ with $|n-n^{\prime}|\ll 1/\lambda$, and then $R_{A}$ is independent of $\lambda$. Already for the Planck number $\bar{n}\sim 1$, $R_{A}$ approaches the classical limit described by Eq. (28) and scales inversely proportional to $2\bar{n}+1$.

An important comment is due here. In classical systems, in a fluctuation-induced escape from a metastable state, with overwhelming probability the system goes over the saddle point on the boundary of the attraction basin to the metastable state, cf. Freidlin and Wentzell (1998). In the case of the oscillator we are considering, as seen from the phase portrait in Appendix B, Fig. 12, the classical switching occurs from a period-3 state to the zero-amplitude state, not to another period-3 state.

A quantum oscillator, however, can directly switch between the period-3 states. This is because the condition of small tunnel splitting compared to the level broadening can be violated near the saddle point, as the tunneling exponent $S_{\rm tun}$ falls off linearly with $g$ near the saddle-point value $g_{s}$, cf. Eq. (36) and Fig. 9. If the tunnel splitting becomes comparable to the level broadening, the oscillator will tunnel between the wells once it has reached the values of $g$ close to $g_{s}$. The overall rate of interwell switching will still be close to that given by Eq. (37), as this equation gives the rate of reaching $g_{s}$, to the logarithmic accuracy.

VI.3 Escape near the bifurcation point where the period-3 state emerges

In the previous sections we considered the weak damping limit, where the intrawell level spacing is much larger than the level widths. If the level broadening is disregarded, the onset of the period-3 vibrational states has no threshold in the amplitude of the driving field for $\delta\omega>0$. In the presence of dissipation, the driving field must have a finite amplitude for the period-3 states to emerge. In terms of the theory of classical dynamical systems, the onset of period-3 vibrations is a bifurcation where the system acquires new stable states, in this case, period-3 states.

Generically, one of the dynamical variables of the system becomes slow near a bifurcation point (the bifurcation parameter values). This variable is an analog of a soft mode. It controls the overall dynamics, as other dynamical variables follow it adiabatically Guckenheimer and Holmes (1997). In the case of a driven quantum oscillator this significantly simplifies the problem of fluctuations Dykman (2007); *Dykman2012. Indeed, since the oscillator dynamics in the rotating frame is Markovian, and since it is controlled by a single dynamical variable, the operator nature of this variable becomes irrelevant. The variable commutes with itself, moreover, it commutes with itself at different instants of time. The fluctuations are still quantum, their intensity is nonzero even for $T\to 0$, but otherwise the dynamics is fully classical.

Classical dynamics of an oscillator resonantly driven close to triple its eigenfrequency is well understood Nayfeh and Mook (2004); Jordan and Smith (2007). This dynamics is briefly summarized and the way to incorporate fluctuations consistent with the kinetic equation (III.1) is described in Appendix B. The important feature for the discussion here is that the period-3 states emerge as a result of three simultaneous saddle-node bifurcations occurring at three points $(Q_{B}^{(i)},P_{B}^{(i)})$ $(i=1,2,3)$ located at the vortices of an equilateral triangle on the $(Q,P)$-plane. At each point $(Q_{B}^{(i)},P_{B}^{(i)})$ there merge a stable and an unstable stationary state in the rotating frame, which correspond to the appropriate stable and unstable period-3 states in the laboratory frame.

If the parameter $\kappa$ is close to its bifurcational value $\kappa_{B}\equiv\kappa_{B}(f)$, one can expand the classical equations of motion near a point $(Q_{B}^{(i)},P_{B}^{(i)})$ in $\delta Q^{(i)}=Q-Q_{B}^{(i)},\;\delta P^{(i)}=P-P_{B}^{(i)}$ and rotate the variables so as to single out the slow variable $z=\delta Q^{(i)}\cos\phi_{B}^{(i)}+\delta P^{(i)}\sin\phi_{B}^{(i)}$. In the quantum Langevin equation for this variable one should keep the leading-order terms in $z,\kappa-\kappa_{B}$. This equation generically has the form Dykman (2012)

[TABLE]

where $\langle\xi_{z}(\tau)\xi_{z}(\tau^{\prime})\rangle=\lambda\kappa(2\bar{n}+1)\delta(\tau-\tau^{\prime})$ is the quantum noise, which is $\delta$-correlated in the slow time $\tau$ (such structure of the noise corresponds to the Markovian dynamics in slow time). The explicit expressions for $Q_{B}^{(i)},P_{B}^{(i)},\phi_{B}^{(i)}$, and the parameters $a_{B},b_{B}$ are given in Appendix B. Here we note that $a_{B}b_{B}<0$.

For $\kappa<\kappa_{B}(f)$ the system (38) has a stable and an unstable stationary state, if we disregard fluctuations. They are located at $\mp z_{\rm st}\,{\rm sgn}\,a_{B}$, where

[TABLE]

At $\kappa=\kappa_{B}(f)$ these states merge together, and for the scaled decay rate $\kappa>\kappa_{B}(f)$ the oscillator does not have period-3 states in the laboratory frame. Equivalently, period-3 states exist in the range of the scaled field amplitude where

[TABLE]

(see Apendix B). In Fig. 11 we plot $\tilde{f}_{B}^{2}=f_{B}^{2}/\kappa$; the parameter $\tilde{f}^{2}=f^{2}/\kappa=F_{0}^{2}/24\omega_{0}\gamma\Gamma$ is advantageous as it does not depend on the frequency detuning of the drive from $3\omega_{0}$.

Fluctuations caused by the noise $\xi_{z}(t)$ in Eq. (38) can lead to escape from the dynamically stable state. Even though our system is far away from thermal equilibrium, the escape rate $W_{\rm esc}$ near the bifurcation point can be found in the same way as for a Brownian particle in a potential well. The dynamics of such a particle in the case where inertial effects can be disregarded (the overdamped limit) is described by the Langevin equation that has the same form as Eq. (38), except that the noise is thermal. Using the Kramers solution of the thermal-equilibrium problem Kramers (1940) and taking into account the explicit form of the intensity of the quantum noise $\xi_{z}(t)$, we obtain

[TABLE]

The exponent of the escape rate scales with the distance to the bifurcation point as $|\kappa-\kappa_{B}|^{3/2}$ and scales with temperature as $(2\bar{n}+1)^{-1}$. In terms of the parameters of the driving field, near the bifurcation point we have $\kappa-\kappa_{B}\approx\Gamma[(\delta\omega)_{B}-\delta\omega]/(\delta\omega_{B})^{2}\approx(f-f_{B})/(df_{B}/d\kappa)$, where $(\delta\omega)_{B}$ is the bifurcational value of the frequency detuning and $f_{B}(\kappa)$ is given by Eq. (39). The scaling (40) with the parameters of the driving field and with temperature holds also for a resonantly driven oscillator near a bifurcation point Dykman (2007); *Dykman2012 and has been seen for such an oscillator in the experiment with Josephson junction oscillators Vijay et al. (2009, 2012).

We note that, even though the zero-amplitude state of the oscillator that displays period tripling is stable, for sufficiently strong driving quantum fluctuations lead to escape from this state into period-3 states. For a strong field, $f^{2}\gtrsim(1+\kappa^{2})/\lambda(2\bar{n}+1)$, one can show that quantum fluctuations essentially wash away the zero-amplitude state. The full analysis of escape from the zero-amplitude state is beyond the scope of this paper.

VII Summary and outlook

This paper describes the dissipative quantum dynamics of a nonlinear oscillator resonantly driven into period tripling. The analysis is based on a minimalistic model. The scaled RWA Hamiltonian in the rotating frame has one parameter, the scaled amplitude of the drive $f$, the dissipation is described by one parameter, the scaled linear friction coefficient $\kappa$, and the only other parameter is the scaled Planck constant $\lambda\propto\hbar$. The model has allowed revealing generic features of the dissipative Floquet dynamics of a multi-state system. Of particular importance is that it revealed the features that remained hidden in the previously explored models.

For period tripling, the RWA Hamiltonian as a function of the phase-space variables $g(Q,P)$ has three symmetrically located wells. In the most interesting regime, the wells contain many bound states. The number of the states is $\propto 1/\lambda\gg 1$. Because of the three-fold symmetry, the intrawell states would be degenerate if there were no interwell tunneling. The tunneling splits the eigenvalues of $g$ into triplets.

Coupling to a thermal reservoir leads to transitions between the states. Where the transition rates ($\propto\kappa$) exceed the exponentially small tunnel splitting (divided by $\hbar$), the distribution over the states is formed in two stages. First, over the dimensionless time $\sim 1/\kappa$ there is formed a distribution over the intrawell states. Then over an exponentially longer time, there is formed a distribution over the wells.

Arguably, the most unexpected feature of the dynamics is that the intrawell distribution is formed by transitions between remote states within the well, i.e., by a nonlocal walk over the intrawell states. This is despite the corresponding transition rates being exponentially small. The effect is a consequence of a very steep decay of the distribution with the increasing RWA energy.

A broader implication of this effect is related to the fact that, even though the distribution is very steep, its logarithm is a smooth function of the number $n$ of an intrawell state [or, equivalently, of the eigenvalue $g_{n}$ of the RWA Hamiltonian]. A standard approach to the analysis of the distribution in this case is based on the eikonal approximation. Here, the logarithm of the distribution is expanded in the state number and the balance equation for the state populations is reduced to the equation for the derivative of the logarithm of the distribution over $g_{n}$, cf. Dykman and Smelyanskii (1988). Such approach is a complete analog of the WKB approximation where the Schrödinger equation is reduced to the equation for the derivative of the logarithm of the wave function. It has been used to analyze numerous problems of chemical kinetics and population dynamics Kamenev (2011); Assaf and Meerson (2017), in which case $n$ stands for the number of species and often has several components that correspond to different species.

Usually, the WKB approximation in quantum mechanics breaks down locally near singularities of the wave front, cf. Landau and Lifshitz (1977); Berry (1989). A similar behavior was found for the distribution of the resonantly driven oscillator Guo et al. (2013b); Peano and Dykman (2014). In contrast, in the present case, the logarithm of the probability distribution remains smooth. However, the local approximation, that underlies the eikonal method, breaks down. Such behavior may be expected also in other types of physical systems away from thermal equilibrium as well as in the problems of chemical kinetics and population dynamics.

The other generic feature of the dynamics of period tripling is the lack of detailed balance. This contrasts the previous results Drummond and Walls (1980); Kryuchkyan and Kheruntsyan (1996); Marthaler and Dykman (2006) for the oscillator driven close to its eigenfrequency or twice the eigenfrequency, where the detailed balance holds at $T=0$. At the technical level, the semiclassical dultiple singularities of the classical Hamiltonian trajectories in the complex time and the branching at the singularities. This significantly complicates the analysis, but also broadens the applicability of the results.

In common with the previously considered driven oscillators, the switching between the period-3 states occurs via transitions over the barrier of $g(Q,P)$ rather than via interwell tunneling, if the decay rate exceeds the tunnel splitting divided by $\hbar$. We have also found the scaling behavior of the switching rate with the distance to the bifurcation point where the period-3 states disappear.

We have not provided the results for the switching from the quiet state. The rate of switching from this state is a peculiar and important problem on its own, starting with the very possibility to switch to any of the period-3 states. The switching rate can be made comparatively large in this case, while still remaining much smaller than the decay rate. A larger switching rate should grossly simplify the observation of the switching in the experiment, thus providing the means to study quantum-fluctuations induced switching where the detailed balance condition is strongly violated. Such an observation should be contrasted with the important experiment on the switching of a resonantly driven quantum oscillator in the regime where the detailed balance was effectively present Vijay et al. (2009). A detailed analysis of the switching from the quiet state is left for the future work.

In conclusion, our results show that the dissipative dynamics of resonant period tripling provides new insights into the effects of dissipation and quantum fluctuations in Floquet systems. Experimentally, these effects can be studied with various nonlinear microwave and optical cavities and also with nanomechanical systems, where quantum regime has been already recently reached Satzinger et al. (2018); Chu et al. (2018).

Acknowledgements.

We are grateful to J. Ankerhold, S. M. Girvin and J. Gosner, who participated in the work at the early stage. We also acknowledge valuable discussions with C. Bruder, N. Lörch and V. Peano. Y.Z. was supported by the National Science Foundation (DMR-1609326). M.I.D. acknowledges partial support from the National Science Foundation (Grants No. DMR-1806473 and and CMMI-1661618).

Appendix A Fixed points of the Hamiltonian $g$

The Hamiltonian function $g(Q,P)$, Eq. (II), which describes the motion in the rotating frame in the rotating-wave approximation, has three minima located at

[TABLE]

Here $Q_{0}$ and the minimal value of $g$ are given by the expressions

[TABLE]

For the detuning of the drive frequency $\delta\omega\equiv\frac{1}{3}\omega_{F}-\omega_{0}>0$, function $g(Q,P)$ has a local maximum at $Q=P=0$, whereas for $\delta\omega<0$ it has a local minimum at $Q=P=0$. For $\delta\omega<0$, as seen from Eq. (A), there is a threshold in the drive amplitude for the local minima to emerge.

Function $g(Q,P)$ also has three saddle points located at

[TABLE]

Here $Q_{s}$ and the saddle-point value of $g$ are given by the expressions

[TABLE]

For $\delta\omega>0$, the saddle points lie on the boundaries between the wells of $g(Q,P)$. For $\delta\omega<0$ the saddle points lie at the boundaries between the wells of $g(Q,P)$ located at nonzero $Q^{2}+P^{2}$, see Eq. (41), and the well at $Q=P=0$.

In the presence of weak dissipation, the three minima of $g(Q,P)$ given by Eq. (41) become stable period-three vibrational states, whereas the local maximum or minimum at the origin becomes a stable zero-amplitude state, see Appendix B.

The scaling used in Sec. II does not allow one to consider the driving on exact resonance. An alternative scaling of the dynamical variables in the rotating frame is

[TABLE]

In these variables, the Hamiltonian in the rotating frame in the RWA becomes $U^{\dagger}HU-i\hbar U^{\dagger}\dot{U}=(2F^{4}/27\gamma^{3})g^{\prime}{}(Q^{\prime}{},P^{\prime}{})$ with

[TABLE]

Clearly, the dynamics is still described by the single parameter, in this case $\zeta^{\prime}{}={\rm sgn}(\delta\omega)/f^{2}$, and by the dimensionless Planck constant $\lambda^{\prime}{}$, which is independent of the frequency detuning $\delta\omega$. The tristability requires $\zeta^{\prime}{}>-1/4$.

Appendix B Classical dissipative dynamics

The Markov quantum kinetic equation in “slow” time (III.1) corresponds to the quantum Langevin equations of motion of variables $Q,P$ of the form

[TABLE]

where $\xi_{Q}(\tau)$ and $\xi_{P}(\tau)$ are Gaussian quantum noises, which are $\delta$-correlated in slow time, $\langle\xi_{Q}(\tau)\xi_{Q}(\tau^{\prime})\rangle=\langle\xi_{P}(\tau)\xi_{P}(\tau^{\prime})\rangle=\lambda\kappa(2\bar{n}+1)\delta(\tau-\tau^{\prime})$, and satisfy the commutation relation $\langle[\xi_{Q}(\tau),\xi_{P}(\tau^{\prime})]\rangle=2i\lambda\kappa\delta(\tau-\tau^{\prime})$, cf. Dykman (2007); *Dykman2012.

In the classical limit the non-commutativity of the noise components can be disregarded and the noise correlators become independent of $\hbar$, $\langle\xi_{Q}(\tau)\xi_{Q}(\tau^{\prime})\rangle=\langle\xi_{P}(\tau)\xi_{P}(\tau^{\prime})\rangle=(3\gamma\kappa k_{B}T/4\omega_{0}^{3}|\delta\omega|)\delta(\tau-\tau^{\prime})$. The classical equations of motion read

[TABLE]

The dynamics described by Eq. (B) in the absence of noise is well understood Nayfeh and Mook (2004); Jordan and Smith (2007). For a sufficiently strong driving (or a small decay rate) these equations can have 7 stationary states. One state is trivial, $Q_{\rm st}=P_{\rm st}=0$. This is a stable state, a focus on the phase plane $(Q,P)$, with the eigenvalues of the linearized motion equal to $-\kappa\pm i$. It corresponds to the local maximum (for $\delta\omega>0$) or minimum (for $\delta\omega<0$) of $g(Q,P)$.

Six other roots have the form $Q^{(j)}_{\rm st}=r^{(j)}_{\rm st}\cos\phi_{\rm st}^{(j)}$, $P^{j}_{\rm st}=r^{(j)}_{\rm st}\sin\phi_{\rm st}^{(j)}$ ($j=0,...,5$). The radii $r^{(j)}_{\rm st}$ take on one of the two values,

[TABLE]

whereas the phases are given by the equation

[TABLE]

We set $j=0,1,2$ for the stationary states $r^{(j)}_{\rm st}=r_{+}$ and $j=3,4,5$ for the states $r^{(j)}_{\rm st}=r_{-}$. A straightforward calculation shows that the eigenvalues of the linearized motion for $j=0,1,2$ have negative real parts, which shows that these states are stable. For the states $j=3,4,5$ the eigenvalues are real and have opposite signs, which means that these states are saddle points on the $(Q,P)$-plane.

From Eq. (49), the triples of the stable states and of the saddle points are located at the vertices of equilateral triangles on the $(Q,P)$ plane . The stable states are further away from the origin $Q=P=0$ , as $r_{+}>r_{-}$. The phase portrait is illustrated in Fig. 12. For comparatively weak damping, the separatrices form spirals that loop around stable states. In the absence of noise, the oscillator will move inside the spiral of the separatrices of the same color toward the stable state inside the corresponding loop. If initially the oscillator is between the loops of different colors, it will move toward the zero-amplitude state.

As seen from Eqs. (B) and (49), the stable and unstable states approach each other pairwise with the increasing $\kappa$ or decreasing $f$. In Fig. 12, the loops of the separatrices “squeeze” so that the squares inside the loops merge with circles on the loops. The merging of a stable state and a saddle is a saddle-node bifurcation. The bifurcational value of the scaled damping rate is given by the condition $r_{+}=r_{-}$,

[TABLE]

Equivalently, one can express the scaled bifurcational field amplitude $f_{B}$ in terms of $\kappa$, see Eq. (39). It is important that the amplitude of the vibrations $r_{B}=(1+\frac{1}{2}f^{2})^{1/2}$ is nonzero at the bifurcation point, i.e., even though there is a threshold in terms of the driving field strength, when the period-3 vibrations emerge in the laboratory frame, they have a nonzero amplitude.

The analysis of the dynamics near the bifurcation point is the same for the three period-3 states. We will choose one of them and introduce a scaled complex vibration amplitude at the bifurcation point $x_{B}=(Q_{B}+iP_{B})$ with $Q_{B},P_{B}$ given by one of the solutions of Eq. (49) for $\kappa=\kappa_{B}$; by construction, $Q_{B}^{2}+P_{B}^{2}={\rm sgn}(\delta\omega)+f^{2}/2$. If we now set in the equations of motion (B) $Q=Q_{B}+\delta Q,\,P=P_{B}+\delta P$, expand the right hand sides in $\delta Q,\delta P$, and write the equation of motion for the combination

[TABLE]

we find that, if we choose the rotation angle as

[TABLE]

the linear in $\delta Q,\,\delta P$ terms drop out from this equation.

The combination of the variables (51) is the slow variable of the system. The linearly independent combination $y=-\delta Q\sin\phi_{B}+\delta P\cos\phi_{B}$ decays with the rate $2\kappa_{B}$ for $\kappa=\kappa_{B}$. Over the dimensionless time $\sim(2\kappa_{B})^{-1}$, this variable approaches its quasi-stationary (adiabatic) value $y_{\rm ad}\equiv y_{\rm ad}(z)$ for a given $z$. One can show that, to the leading order,

[TABLE]

The resulting equation for the slow variable $z$, to the leading order, has the form (38) with the coefficients

[TABLE]

where $X_{B}=x_{B}\exp(-i\phi_{B})$. It is important that the parameters $a_{B}$ and $b_{B}$ are the same for all three period-three states. Indeed, changing from one of the three stable states to another corresponds to the rotation $x_{B}\to x_{B}\exp(\pm 2\pi i/3)$. As seen from Eq. (52), such rotation leads to the angle $\phi_{B}$ being incremented by $\pm 2\pi/3$, as expected from the symmetry arguments. Therefore $X_{B}$ is the same for all three states as is also $\exp(3i\phi_{B})$.

Appendix C Asymptotic behavior of the semiclassical matrix elements

Finding the asymptotic behavior of the matrix elements $a_{m}(g_{n})$ in Eq. (15) for large $|m|$ requires understanding the analytical properties of the function $Q(\tau;g)+iP(\tau;g)\equiv(2\lambda)^{1/2}a(\tau;g)$. This function is determined by the classical Hamiltonian dynamics with the Hamiltonian $g(Q,P)$. In the explicit form, the Hamiltonian equations (16) for $Q(\tau;g)$ and $P(\tau;g)$ read

[TABLE]

Here and in what follows we consider the case $\delta\omega>0$.

For real time, Im $\tau=0$, the functions $Q(\tau;g)$ and $P(\tau;g)$ describe smooth Hamiltonian trajectories, which in the range of $g$ we are interested in, $g_{\min}<g<g_{s}$, are shown in Fig. 1 (c). However, $Q(\tau;g)$ and $P(\tau;g)$ generally become complex for complex time and have singularities where $|Q|$ and $|P|$ go to infinity. Because $g(Q,P)$ is quartic in $Q,P$, these singularities can be reached in a complex time $\tau_{\rm sing}$ with finite real and imaginary parts. From Eq. (15), for large $|m|$ the Fourier components $a_{m>0}$ and $a_{m<0}$ are given by the values of $\tau_{\rm sing}$ with the smallest in the absolute value negative and positive imaginary part, respectively.

By symmetry, the values of Im $\tau_{\rm sing}$ are the same for all three wells of $g(Q,P)$. It is most convenient to find them by considering the classical trajectory for the $\nu=0$-well of $g(Q,P)$, which is centered on the axis $P=0$ in Figs. 1 (b) and (c). This trajectory in real time is shown by the red line in Fig. 13. For a given $g$, the momentum $P$ as a function of $Q$ on this trajectory has the form

[TABLE]

as seen from the explicit expression (II) for $g(Q,P)$.

The time on the classical trajectory can be continued to the complex plane in such a way that $Q(\tau;g)$ will remain real and will go to infinity. To this end, we start from the point on the real trajectory where $Q$ is maximal, $Q=Q_{\max}$ (the utmost right point on the red trajectory in Fig. 13). We choose $\tau=0$ for $Q=Q_{\max}$ and then move along the imaginary axis on the $\tau$-plane. As seen from Eq. (C), in this case $Q(\tau;g)$ remains real, whereas $P(\tau;g)$ becomes purely imaginary.

The real trajectory (C) spins clockwise, and therefore $\dot{P}(\tau;g)<0$ for $\tau=0$ (i.e., for $Q=Q_{\max}$). Therefore the sign of Im $P$ is opposite to the sign of Im $\tau$, when we start incrementing Im $\tau$. Since $P(\tau=0;g)=0$, function $P(\tau;g)$ is purely imaginary for small $|{\rm Im}~{}\tau|$,

[TABLE]

From Eqs. (C), $Q(\tau;g)$ remains real and $P(\tau;g)$ remains purely imaginary for purely imaginary $\tau$. Moreover, both $Q$ and $P$ increase in the absolute value as $|\tau|$ increases independent of the sign of Im $\tau$.

In the upper half-plane of the $\tau$-plane, the value of $\tau$ at which $Q\rightarrow\infty$ is given by $i\tau_{\infty}$,

[TABLE]

where we have used Eq. (C) and chose Im $P<0$, in agreement with Im $\tau>0$; note that $B(Q{})>0$ and Re $\dot{Q}=0$ for $Q>Q_{\max}$.

When $Q$ goes to infinity, $|P|$ also goes to infinity and so does $|Q+iP|$ as well. Indeed, as seen from Eq. (C),

[TABLE]

From this expression, $Q+iP\propto Q$ or $\sqrt{Q}$ for $Q\to\infty$ depending on the sign of Im $P$ or, equivalently, on the sign of -Im $\tau$.

Once we know the locations of the singularities of the function $Q(\tau;g)+iP(\tau;g)$, we are ready to perform the integral in Eq. (15) in the limit of large $m$. For the case of $m<0(m>0)$, we move the integration contour to the upper (lower) half-plane of the complex-$\tau$ plane; the case where the integral is shifted to the upper half-plane is shown in Fig. 13. With our choice of the time on the orbit, the singularities of $Q(\tau;g)+iP(\tau;g)$ are on the imaginary-$\tau$ axis at $\tau=\pm i\tau_{\infty}$.

To evaluate the contour integral, we need to know the function $Q+iP$ near the singularities. From Eqs. (56) and (57), for large $Q$ we have $\dot{Q}=P(Q{})\sqrt{B(Q{})}\approx\mp iQ\sqrt{16fQ^{3}/3}$, where the upper and lower signs correspond to the upper and lower half-planes of the $\tau$-plane. Then for $|\tau-i\tau_{\infty}|\gg\tau_{\infty}$

[TABLE]

whereas for $|\tau+i\tau_{\infty}|\ll\tau_{\infty}$

[TABLE]

The above equations show that the singularities at $\pm i\tau_{\infty}$ are branching points. We choose the branching cuts so that the contour integral does not cross them. For the chosen contour, the integral of $\exp[-im\omega(g)\tau]\,[Q(\tau;g)+iP(\tau;g)]d\tau$ in Eq. (15) is reduced to the integral around the branch cut. Other contributions to the contour integral can be disregarded. Indeed, since $Q(\tau;g)$ and $P(\tau;g)$ are periodic with period $2\pi/\omega(g)$, the integrals along the vertical lines that go from $\pm\pi/\omega(g)$ in Fig. 13 cancel each other. The integral along the horizontal sections of the contour with large $|{\rm Im}~{}\tau|$ is exponentially small for large $|m|$.

The calculation of the integral around the branch cut is straightforward with the account taken of Eqs. (58) and (59). It gives Eqs. (III.2) and (III.2). The different dependence of the prefactors in these equations on $m$ for $m>0$ and $m<0$ comes from the different dependence of $a(\tau;g)\propto Q(\tau;g)+iP(\tau;g)$ on $\tau$ near $\pm i\tau_{\infty}$ in Eqs. (58) and (59).

Implied in the above calculation was an assumption that the singularities at $\pm i\tau_{\infty}$ are the closest singularities to the real-$\tau$ axis. Generally, there are other singularities of $Q(\tau;g)$ and $P(\tau;g)$ for complex $\tau$. If we start from the leftmost point $Q_{\rm min}$ on a real trajectory in the $\nu=0$-well of $g(Q,P)$ (the red trajectory in Fig. 13) and move along the imaginary time axis, $Q$ diverges for Im $\tau=\pm(\tau_{\rm tun}-\tau_{\infty})$. Here, $\tau_{\rm tun}\equiv\tau_{\rm tun}(g)$ is the imaginary part of the tunneling time to go from the classical real-time trajectory of the $\nu=0$ well of $g(Q,P)$ with a given $g$ to the classical real-time trajectory with the same $g$ in the $\nu=1$ or $\nu=2$ well; see Sec. VI for details. However, an explicit calculation shows that $\tau_{\infty}<\tau_{\rm tun}-\tau_{\infty}$, and therefore the singularities at Im $\tau=\pm(\tau_{\rm tun}-\tau_{\infty})$ do not affect the large-$|m|$ asymptotic of $a_{m}(g)$. Starting from various points on the real-time $\nu=0$-trajectory and moving along the imaginary time axis, we checked numerically that the singularities at $\pm i\tau_{\infty}$ are indeed the closest to the real time axis.

We show in Fig. 14 a comparison between the asymptotic expressions for the semiclassical matrix elements in Eq. (III.2) and (III.2) and the results obtained using numerical integration of the Hamiltonian equations of motion inside a well of $g(Q,P)$. Both the decay exponents and the prefactors of the matrix elements agree well with the numerics. We note that the matrix elements always satisfy $|a_{-m}|>|a_{m}|$ for $m>0$. This means that the rate $W_{n\,n-|m|}$ of transitions down in quasienergy is always larger than the rate $W_{n\,n+|m|}$ of transitions up, so that the stationary distribution over the intrawell states monotonically falls off with the increasing $g$.

The discussed here behavior of $Q+iP$ as a function of complex time for period tripling is qualitatively different from the behavior of the analogous function for a parametric modulated oscillator Marthaler and Dykman (2006) and for a nonlinear oscillator driven close to its eigenfrequency Dykman (2012); Guo et al. (2013b). This function has different forms for the parametric and resonantly driven oscillators, but in the both cases it is described by the Jacobi elliptic functions. It is double periodic in complex time, with two poles within the parallelogram of periods and with no branch points. This difference leads to a qualitative difference of the intrawell distribution and the probability of interstate switching, as discussed in the main text.

Appendix D The intrawell distribution in the classical limit

In this section, we compare the RWA-energy dependent inverse effective temperatures $R^{\prime}(g)$ and the scaled logarithms of the probability distribution $R(g)$ obtained by solving the full semiclassical balance equation (27) and by solving Eq. (28) that formally refers to the classical limit. The results are shown in Fig. 15.

For a weak drive (the upper left panel of Fig. 15), the behavior of $R^{\prime}(g)$ is similar to that in the classical limit even at zero temperature. Here, the function $R^{\prime}(g)$ monotonically increases from $g=g_{\rm min}$ to $g=g_{s}$ and the function $(2\bar{n}+1)R(g)$ weakly changes with temperature. In contrast, for a relatively large drive (the lower left panel of Fig. 15), $R^{\prime}(g)$ at low temperature decreases as $g$ increases, which is significantly different from the behavior in the classical limit. As the temperature increases, both $R^{\prime}(g)$ and $R(g)$ converge to the classical limit.

Appendix E The imaginary tunneling time

The imaginary tunneling time $\tau_{\rm tun}(g)$ is the imaginary part of the time of moving, with a given $g$, along a Hamiltonian trajectory in the classically inaccessible region of the $(Q,P)$ plane between the wells of $g(Q,P)$. It is given by Eq. (36), which we repeat here to simplify the reading,

[TABLE]

where $P(Q{})\equiv P(Q,g)$ and $B(Q{})\equiv B(Q,g)$ are given in Eq. (C). For concreteness, we will consider tunneling from the $\nu=0$ well of $g(Q,P)$, which has the minimum on the $P=0$-axis and is symmetric with respect to the reflection $P\to-P$. We assume that $\delta\omega>0$.

To calculate $\tau_{\rm tun}$ one has to establish the limits of the integral in Eq. (E). For a given RWA energy $g$, these limits are determined by the boundaries of the classical real-time orbits $g(Q,P)=g$ shown in Fig. 1 (c). We have to consider the orbits in the wells $\nu=0$ and $\nu=1$ or 2. We start with the orbits in the $\nu=0$-well of $g(Q,P)$. Such orbits have two turning points $\dot{Q}\equiv P\sqrt{B(Q)}=0$ on the $P=0$ axis. These points, $Q_{\min}(g)$ and $Q_{\max}(g)$ are given by the equation $g(Q,0)=g$.

The points $Q_{\min}(g)$ and $Q_{\max}(g)$ are the only turning points in the region of the RWA energy $g<g_{\rm cr}$, because $B(Q,g)>0$ in this region. The shape of an orbit is close to elliptical, as seen in the left panel of Fig. 16. The point $Q_{\min}(g)$ is the leftmost point on the orbit.

For $g_{\rm cr}<g<g_{s}$, function $B(Q,g)$ becomes zero on the classical orbit at a point $Q=Q_{B}(g)<Q_{\rm min}(g)$. At $Q=Q_{B}$ we have $\dot{Q}=0$, and thus $|dP/dQ|\to\infty$, but $P\neq 0$. Respectively, the orbit becomes horse-shoe like as illustrated in the right panel of Fig. 16. The point $Q_{B}(g)$ is the leftmost point on the orbit for $g>g_{\rm cr}$.

At $g=g_{\rm cr}$ the points $Q_{B}$ and $Q_{\min}$ merge, that is, the function $B(Q,g_{\rm cr})$ is zero at $Q=Q_{\min}(g_{\rm cr})$. From Eq. (C) it follows that at this point $A(Q)=0$. These conditions give for $g_{\rm cr}$ and $Q_{\rm cr}\equiv Q_{\min}(g_{\rm cr})$ the expressions

[TABLE]

The integration over $Q$ in Eq. (E) for tunneling from the $\nu=0$-well goes in the negative direction from the leftmost point on the classical orbit of this well to the rightmost point $Q_{\rm fin}$ of the classical orbit of the well $\nu=1$ or, equivalently, $\nu=2$. By construction, $\dot{Q}=0$ at the rightmost point $Q_{\rm fin}$ of the classical orbit of the well $\nu=1$ (or $\nu=2$). Since the orbits of the wells $\nu=1$ and $\nu=2$ do not cross, $P\neq 0$ at this point. Therefore the point $Q_{\rm fin}\equiv Q_{\rm fin}(g)$ is given by the intermediate root of the equation $B(Q,g)=0$.

The equation $B(Q,g)=0$ has three real roots. The smallest one gives the leftmost point on the orbits of the $\nu=1,2$ wells where $\dot{Q}=0$. The largest root is $Q_{B}$; for $g>g_{\rm cr}$ it lies on the classical orbit of the $\nu=0$ well, whereas for $g<g_{\rm cr}$ it lies in the classically forbidden region between the wells. In the region between $Q_{\rm fin}(g)$ and $Q_{B}(g)$ we have $B(Q,g)<0$, and thus the momentum $P(Q,g)$ is complex, as seen from Eq. (C).

Taking into account the difference between the regions $g<g_{\rm cr}$ and $g>g_{\rm cr}$, we write the expression for the imaginary tunneling time as

[TABLE]

Here we have set $B(Q,g)^{1/2}=i|B(Q,g)|^{1/2}$ for $B(Q,g)<0$. Then, taking into account that Im $P(Q,g)<0$ on the tail of the wave function of the $\nu=0$-well, we find from Eq. (C) that, for $Q_{\rm fin}(g)<Q<Q_{B}(g)$ the momentum in the classically forbidden region takes the form

[TABLE]

For $g<g_{\rm cr}$ and $Q_{B}(g)<Q<Q_{\min}(g)$, we have $A(Q)<0,B(Q,g)>0$, and $P(Q,g)$ is purely imaginary:

[TABLE]

Clearly, $\tau_{\rm tun}(g)<0$.

The integrand in Eq. (E) is singular at the saddle-point value of $g$ but the integral is finite. For $g=g_{s}$, the classical trajectory of the well $\nu=0$ has common points $(Q=-Q_{s}/2,P=\pm\sqrt{3}Q_{s}/2)$ with the classical trajectories of the wells $\nu=1,2$, where $Q_{s}$ is the position of the saddle point on the axis $P=0$; as seen from Eq. (C), $Q_{s}=[f-(f^{2}+4)^{1/2}]/2$. Clearly, $S_{\rm tun}(g)\to 0$ for $g\to g_{s}$. From Eq. (E), near $g=g_{s}$ we have, to the leading order in $g_{s}-g$

[TABLE]

One can show that $\tau_{\rm tun}(g)$ is linear in $g-g_{s}$ near $g_{s}$.

The integrands in Eq. (E) are also singular at $g_{\rm cr}$. In a narrow region around $Q=Q_{B}(g)$ the WKB approximation breaks down. However, the contribution of the corresponding region is small for small $\lambda$; a somewhat cumbersome analysis shows that $\tau_{\rm tun}(g)$ as given by Eq. (E) is smooth and its first derivative is at least continuous at $g_{\rm cr}$.