Enhanced excitation of a driven bistable system induced by spectrum degeneracy

TL;DR

This paper investigates how spectrum degeneracy enhances excitation in driven bistable systems, revealing new critical parameters, state superpositions, and linking multi-photon transitions with tunneling effects.

Contribution

It introduces the concept of spectrum degeneracy's impact on excitation and generalizes Keldysh theory for bistable systems, connecting multi-photon transitions with tunneling.

Findings

Spectrum degeneracy drastically alters system statistics and kinetics.

A new critical quasienergy parameter emerges in nearly degenerate states.

Multi-photon transition amplitude equals tunneling amplitude.

Abstract

The non-equilibrium statistics and kinetics of a simple bistable system (resonantly driven nonlinear oscillator coupled to reservoir) have been investigated by means of master equation for the density matrix and quasiclassical Fokker-Planck equation in quasienergy space. We found out that the system's statistical and kinetic properties drastically change when the quasienergy states become nearly degenerate and the occupation of the most excited state is strongly enhanced. It has been revealed that in nearly degenerate case a new critical quasienergy parameter emerges. Below the critical quasienergy value the eigenstates are superpositions of the quasiclassical states from different phase space regions, while above this value the eigenstates correspond to only one particular region of the phase space. We have also generalized Keldysh theory for ionization of atoms in the electromagnetic field for bistable systems. It has been demonstrated that Keldysh parameter in bistability region is large when pumping intensity is smaller than the critical value. It has been shown by direct calculations that multi-photon transition amplitude coincides with the tunneling amplitude. So, multi-photon transitions and tunneling between the regions of the phase space are just the same effects. We also demonstrated that for bistable systems the Keldysh parameter logarithmically depends on the external field amplitude.

Equations

$$\hat{H}=-\Delta {\hat{a}}^{†}\hat{a}+\frac{\alpha }{2}({\hat{a}}^{†}\hat{a}{)}^2+f(\hat{a}+{\hat{a}}^{†}).$$

$$H_{\mathrm{int}}={\hat{\xi }}^{†}\hat{a}+\hat{\xi }{\hat{a}}^{†},$$

$$\begin{array}{c}⟨\hat{\xi }(t){\hat{\xi }}^{†}(t^{\prime })⟩=\gamma (N+1)\delta (t-t^{\prime }),\\ ⟨{\hat{\xi }}^{†}(t)\hat{\xi }(t^{\prime })⟩=\gamma N\delta (t-t^{\prime }),\end{array}$$

$$\begin{array}{c}\frac{d{P}_n}{dt}=\sum _{n^{\prime }}{w}_{n{n}^{\prime }}{P}_{n^{\prime }}-w_{n^{\prime }n}{P}_n,\\ w_{n{n}^{\prime }}=\gamma \left[(N+1)|⟨n|\hat{a}|n^{\prime }⟩{|}^2+N|⟨n^{\prime }|\hat{a}|n⟩{|}^2\right].\end{array}$$

$$\begin{array}{c}\frac{\partial P_i(E)}{\partial t}=\frac{1}{T_i(E)}\frac{\partial J_i(E)}{\partial E},\\ J_i(E)=\vartheta K_i(E)P_i(E)+Q{D}_i(E)\frac{\partial P_i}{\partial E}.\end{array}$$

$$\frac{1}{2\pi }\oint pdq=2\pi n.$$

$$\begin{array}{c}{ϵ}_n^{(0)}=-\Delta n+\frac{\alpha }{2}{n}^2=\frac{2{\Delta }^2}{\alpha }\frac{n}{m}\left(\frac{n}{m}-1\right)\\ m\equiv \frac{2\Delta }{\alpha }.\end{array}$$

$${ϵ}_1(-n,f)={ϵ}_3(m-n,f),$$

$$n_3(ϵ,f)-n_1(ϵ,f)=m.$$

$${\omega }_R(E)\sim \frac{1}{T(E)}{e}^{-S_{\mathrm{tunn}}(E)},$$

$$\begin{gathered} S_{\mathrm{tunn}}(E)=\frac{1}{2}\oint p_{\mathrm{im}}dq= \\ =\frac{m}{2}{\int }_{s_1}^{s_2}\mathrm{acosh}\left\{\frac{E+\frac{s^2}{2}-\frac{s^4}{8}}{\sqrt{2\beta }s}\right\}sds, \end{gathered}$$

$$S_{\mathrm{tunn}}(E)=\frac{m}{2}\ln \frac{1}{\beta }\sqrt{1+2E}+O(1).$$

$$S_{\mathrm{tunn}}(E_n)=\frac{m}{2}\ln \frac{1}{\beta }\left(1-\frac{2n}{m}\right),$$

$$H_n=\left(\begin{array}{cc}{ϵ}_{n_1}& -{\omega }_R^n\\ -{\omega }_R^n& {ϵ}_{n_3}\end{array}\right),$$

$$\delta {ϵ}_n={ϵ}_{n_1}-{ϵ}_{n_3}=-\frac{2\pi \delta m}{T(E_n)}.$$

$$\begin{array}{c}|n_{+}⟩=c_n|n_1⟩+s_n|n_3⟩,\\ |n_{-}⟩=-s_n|n_1⟩+c_n|n_3⟩.\end{array}$$

$$\begin{array}{c}{c}_n\equiv \cos {\theta }_n,s_n\equiv \sin {\theta }_n,\\ \tan 2{\theta }_n=\frac{2{\omega }_R}{\delta {ϵ}_n}=\frac{1}{\pi \delta m}{e}^{-S_{\mathrm{tunn}}(E_n)}.\end{array}$$

$$\frac{1}{\pi \delta m}{e}^{-S_{\mathrm{tunn}}(E_c)}=1.$$

$$\delta {ϵ}_n^{(2)}=\frac{f^2}{\alpha }\cdot \frac{(m+1)}{(m-2n{)}^2-1}.$$

$${\omega }_R^{n,m-n}=\frac{V_{n,n+1}\dots V_{m-n-1,m-n}}{({ϵ}_n^{(0)}-{ϵ}_{n-1}^{(0)})\dots ({ϵ}_n^{(0)}-{ϵ}_{m-n-1}^{(0)})}\propto \Delta {\beta }^{\frac{m}{2}-n},$$

$$\begin{array}{c}\frac{\partial P_n^s}{\partial t}=\sum _{n^{\prime }{s}^{\prime }}{w}_{n{n}^{\prime }}^{s{s}^{\prime }}{P}_{n^{\prime }}^{s^{\prime }}-w_{n^{\prime }n}^{s^{\prime }s}{P}_n^s.\end{array}$$

$$\begin{array}{c}\frac{\partial P_n^{+}}{\partial t}=\sum _{n^{\prime }}{W}_{n{n}^{\prime }}^{+}{P}_{n^{\prime }}^{+}-W_{n^{\prime }n}^{+}{P}_n^{+}+w_{n{n}^{\prime }}^{+-}{P}_{n^{\prime }}^{-}-w_{n{n}^{\prime }}^{-+}{P}_{n^{\prime }}^{-},\\ \frac{\partial P_n^{-}}{\partial t}=\sum _{n^{\prime }}{W}_{n{n}^{\prime }}^{-}{P}_{n^{\prime }}^{-}-W_{n^{\prime }n}^{-}{P}_n^{-}+w_{n{n}^{\prime }}^{-+}{P}_{n^{\prime }}^{+}-w_{n{n}^{\prime }}^{+-}{P}_{n^{\prime }}^{+},\end{array}$$

$$\begin{array}{c}{W}_{n{n}^{\prime }}^{+}=(N+1)\left(c_{n^{\prime }}^2|a_{n_1{n}_1^{\prime }}{|}^2+s_{n^{\prime }}^2|a_{n_3{n}_3^{\prime }}{|}^2\right)+N\left(c_{n^{\prime }}^2|a_{n_1^{\prime }{n}_1}{|}^2+s_{n^{\prime }}^2|a_{n_3^{\prime }{n}_3}{|}^2\right),\\ W_{n{n}^{\prime }}^{-}=(N+1)\left(s_{n^{\prime }}^2|a_{n_1{n}_1^{\prime }}{|}^2+c_{n^{\prime }}^2|a_{n_3{n}_3^{\prime }}{|}^2\right)+N\left(s_{n^{\prime }}^2|a_{n_1^{\prime }{n}_1}{|}^2+c_{n^{\prime }}^2|a_{n_3^{\prime }{n}_3}{|}^2\right),\\ w_{n{n}^{\prime }}^{+-}=(N+1)(c_n^2{s}_{n^{\prime }}^2|a_{n_1{n}_1^{\prime }}{|}^2+s_n^2{c}_{n^{\prime }}^2|a_{n_3{n}_3^{\prime }}{|}^2-2c_n{c}_{n^{\prime }}{s}_n{s}_{n^{\prime }}\mathrm{Re}{a}_{n_1{n}_1^{\prime }}{a}_{n_3^{\prime }{n}_3}^{\ast })\\ +N(c_n^2{s}_{n^{\prime }}^2|a_{n_1^{\prime }{n}_1}{|}^2+s_n^2{c}_{n^{\prime }}^2|a_{n_3^{\prime }{n}_3}{|}^2-2c_n{c}_{n^{\prime }}{s}_n{s}_{n^{\prime }}\mathrm{Re}{a}_{n_1^{\prime }{n}_1}^{\ast }{a}_{n_3{n}_3^{\prime }}),\\ w_{n{n}^{\prime }}^{-+}=(N+1)(s_n^2{c}_{n^{\prime }}^2|a_{n_1{n}_1^{\prime }}{|}^2+c_n^2{s}_{n^{\prime }}^2|a_{n_3{n}_3^{\prime }}{|}^2-2s_n{s}_{n^{\prime }}{c}_n{c}_{n^{\prime }}\mathrm{Re}{a}_{n_1{n}_1^{\prime }}{a}_{n_3^{\prime }{n}_3}^{\ast })\\ +N(s_n^2{c}_{n^{\prime }}^2|a_{n_1^{\prime }{n}_1}{|}^2+c_n^2{s}_{n^{\prime }}^2|a_{n_3^{\prime }{n}_3}{|}^2-2s_n{s}_{n^{\prime }}{c}_n{c}_{n^{\prime }}\mathrm{Re}{a}_{n_1^{\prime }{n}_1}^{\ast }{a}_{n_3{n}_3^{\prime }}).\end{array}$$

$$\begin{array}{c}\frac{\partial P^{\pm }}{\partial t}=\frac{1}{T}\frac{\partial J_{\pm }}{\partial E}+{\Lambda }_{\mathrm{tunn}}^{\pm }.\end{array}$$

$$J_{\pm }=\vartheta K^{\pm }+Q{D}^{\pm }\frac{\partial P^{\pm }}{\partial E},$$

$$\begin{gathered} {\Lambda }^{+(-)}(E)\equiv {\int }_{E_{\mathrm{sep}}}^{E_1}d{E}^{\prime }T(E^{\prime })[w_{E{E}^{\prime }}^{+-(-+)}{P}^{-(+)}(E^{\prime }) \\ -w_{E^{\prime }E}^{-+(+-)}{P}^{+(-)}(E)]. \end{gathered}$$

$$\bar{P}(E)=\exp \left\{-\frac{\vartheta }{Q}{\int }_{E_{\mathrm{sep}}}^E\frac{K_1+K_3}{D_1+D_3}d{E}^{\prime }\right\},E_{\mathrm{sep}}<E<E_c.$$

$$\begin{array}{c}{P}^{\pm }(E)=P_0^{\pm }(E)+\delta P^{\pm }(E),\end{array}$$

$$\begin{array}{c}{P}_0^{+}(E)=\bar{P}(E_c)\exp \left\{-\frac{\vartheta }{Q}{\int }_{E_c}^E\frac{K_1}{D_1}d{E}^{\prime }\right\},E_c<E<E_1\\ P_0^{-}(E)=\bar{P}(E_c)\exp \left\{-\frac{\vartheta }{Q}{\int }_{E_c}^E\frac{K_3}{D_3}d{E}^{\prime }\right\}{E}_c<E<\infty .\end{array}$$

$$\begin{array}{c}\delta P^{+}=-P_0^{+}{\int }_{E_{\mathrm{sep}}}^E\frac{d{E}^{\prime }}{Q{D}_1(E^{\prime })P_0^{+}}{\int }_{E_1}^{E^{\prime }}{\Lambda }_{\mathrm{tunn}}^{+}(E^{\prime \prime })T(E^{\prime \prime })d{E}^{\prime \prime },\\ \delta P^{-}=P_0^{-}{\int }_{E_{\mathrm{sep}}}^E\frac{d{E}^{\prime }}{Q{D}_3(E^{\prime })P_0^{-}}{\int }_{E^{\prime }}^{\infty }{\Lambda }_{\mathrm{tunn}}^{-}(E^{\prime \prime })T(E^{\prime \prime })d{E}^{\prime \prime },\end{array}$$

Full text

I Introduction

Nowadays bistability is one of the most pronounced phenomena in modern optics and electronics. It has broad applications in all–optical logic and memories performance. So, controllable changing of different stable states occupation and the control of transition rates between them are among the most important problems. The solution of these problems relies on one’s knowledge of optimal perturbation which transfers the system from one stable state to another one. Another problem is how to control the structure of stable states by changing the system parameters such as external field frequency and intensity. Bistability has been widely studied in different experimental setups: cold clouds of atoms inside the optical cavity Rosenberger et al. (1991), Wang et al. (2001), exciton–polariton modes in microcavities with external pumping Gippius et al. (2007), fiber ring cavities Li et al. (2017) and mesoscopic Josephson junction array resonators Muppalla et al. (2018) in the external field, etc. The results of many of these experiments can be understood by investigation of bistable single–mode system with Kerr–like nonlinearity. Resonance response of a bistable single–mode system with a Kerr–like nonlinearity to the external field is described by a model of a driven nonlinear oscillator interacting with dissipative environment Drummond and Walls (1980). This model also describes an atomic system with several energy levels coupled to the cavity mode after adiabatically excluding atomic variables Gothe et al. (2019), Shirai et al. (2018), Dong et al. (2019).

Moreover, the model of a driven nonlinear oscillator is a minimal model of a bistable driven system out of equilibrium. In the quasiclassical limit, its statistics and kinetics can be described by a Fokker–Planck equation (FPE) Maslova et al. (2007). By means of FPE, it is possible to find the stationary distributions, relaxation rates and the occupations of two classical stable states of the oscillator. For the quantum oscillator, the non–equilibrium statistics and relaxation kinetics at different temperatures were studied numerically using the rate equation Risken et al. (1987), Risken and Vogel (1988). Also for the quantum oscillator the exact Glauber–Sudarshan function of the steady state was obtained Drummond and Walls (1980) in the case of zero bath temperature.

In the described model of the bistable driven system, it has been shown Maslova et al. (2007) that there exist different regions of the classical phase space with degenerate energies. Among the quantum effects, the effects of tunneling transitions between these regions are of particular interest. As mentioned in Maslova et al. (2019), tunneling increases the population of the higher amplitude state, which is also a squeezed state. Also it increases the relaxation rate to the stationary distribution. From the quasiclassical point of view, tunneling can lead to hybridization of quasiclassical states from different regions of the classical phase plane. It can be shown to be very strong in the case of integer or half–integer detuning–nonlinearity ratio, when the quasiclassical states from the different regions of the phase space become degenerate. This also corresponds to the multi–photon resonance between the real energy levels of the driven system. The hybridization of the quasienergy states from different regions of the phase space in the case of multi–photon resonance can strongly change the kinetics of the considered system. However, the kinetics of the considered bistable system with the eigenstates which are superpositions of the states from different regions of the classical phase space drastically differs from the non–degenerate case investigated previously Maslova et al. (2019) and are not studied yet.

In the case of strong hybridization between the states from different regions of the phase plane, the transition rate between them can be explained by generalization of Keldysh theory for ionization of atoms in electromagnetic field Keldysh (1965) for bistable systems. The Keldysh theory explains the interplay between tunneling and multi–photon ionization. However, for bistable systems the correspondence between multi–photon transitions and tunneling effects is not clear. To understand it, one should define the Keldysh parameter $\gamma_{K}$ as the ratio of <<tunneling time>> to the period of motion along the classical trajectory. It will be demonstrated that when the field intensity is much smaller than the critical value defining the range of bistability, the Keldysh parameter is large, $\gamma_{K}\gg 1$ and tunneling probability is just the same as the probability of multi–photon excitation. It will be also shown that $\gamma_{K}$ logarithmically depends on the strength of the external field $f$ while in the case of ionization of atoms it is proportional to $f^{-1}$.

II The simple model of bistable driven system

The Hamiltonian of a driven resonant mode with Kerr–like nonlinearity in the rotating–frame approximation reads

[TABLE]

Here $\Delta$ is the detuning between the driving force and the resonant oscillator frequency, $\alpha$ is the blue shift due to nonlinearity, and $f$ is proportional to the amplitude of the driving force.

In the classical limit, one should replace the operators $\hat{a}$, $\hat{a}^{\dagger}$ in (1) with classical field amplitudes $a,a^{*}$ to obtain the classical Hamiltonian. The classical phase portrait of the system is shown on Fig. 1: the classical trajectories in the $a$ plane are given by the contour lines of the classical Hamiltonian. Each classical trajectory corresponds to a certain quasienergy value $\epsilon$. The only dimensionless parameter governing the system classical dynamics is $\alpha f^{2}/\Delta^{3}\equiv\beta$,

which can be treated as the rephasing rate of the nonlinear driven oscillator Borenstein and Lamb (1972). This parameter can also be identified with the Dicke cooperation parameter determining the typical rate of the intensity growth of a superradiance pulse. Bistability appears when $0<\beta<4/27$.

For the quantum Hamiltonian, there exists another dimensionless parameter $m\equiv 2\Delta/\alpha$. The quasiclassical limit is acquired at large $m$.

A key feature of the driven nonlinear oscillator is the presence of two stable stationary states which means bistability. Another important feature is the presence of a self–intersecting trajectory called separatrix, which divides the phase plane into three regions $1$, $2$ and $3$ and passes through the unstable stationary point $S$. The regions $1$ and $2$ contain the stationary states with smaller and larger amplitude respectively. The quasienergies of the states 1,2,S are denoted by $\epsilon_{1}$, $\epsilon_{2}$ and $\epsilon_{\mathrm{sep}}$, and they always obey the inequality $\epsilon_{2}<\epsilon_{\mathrm{sep}}<\epsilon_{1}$. In further discussion, we will also use the dimensionless quasienergy defined as $E=\alpha\epsilon/\Delta^{2}$, $E_{1}=\alpha\epsilon_{1}/\Delta^{2}$, $E_{2}=\alpha\epsilon_{2}/\Delta^{2}$, $E_{\mathrm{sep}}=\alpha\epsilon_{\mathrm{sep}}/\Delta^{2}$,

The states of the quantum Hamiltonian in the limit of large numbers of excitation quanta correspond to a discrete set of classical trajectories on the phase portrait. The corresponding energies can be obtained from the Bohr–Sommerfeld quantization rule.

The interaction with the environment can be described by the following interaction Hamiltonian:

[TABLE]

where $\hat{\xi},\hat{\xi}^{\dagger}$ are the operators of random force with correlation functions defined as

[TABLE]

where $N$ is the number of noise quanta and $\gamma$ is the damping caused by the interaction with the environment. The kinetics of the quantum bistable driven oscillator in the limit of weak coupling with the environment can be treated in the diagonal approximation for the density matrix. In this approximation, one obtains the rate equation dealing with probabilities $P_{n}$ of occupation of the $n$–th quasienergy state:

[TABLE]

If each quasienergy state can be uniquely attributed to one of the regions of the phase space, in the limit of large number of excitation quanta the rate equation transforms to the classical Fokker–Planck equation in the quasienergy space:

[TABLE]

where $T_{i}(E)$ is the period of motion along the classical trajectory, the coefficients $K_{i}(E)$ and $D_{i}(E)$ are the drift and diffusion coefficients in quasienergy space. The coefficients $T_{i}(E),K_{i}(E)$ and $D_{i}(E)$ are defined in the Appendix A. The dimensionless parameters $\vartheta,Q$ are defined as $\vartheta=\gamma/\Delta$ and $Q=\vartheta\alpha(N+1/2)/\Delta$. The probability densities $P_{i}(E)$ are the continuous limits of $P_{n}$, where the probabilities $P_{n}$ are considered as functions of dimensionless quasienergy $E\equiv\alpha\epsilon_{n}/\Delta^{2}$, and $i$ denotes the region of the phase space according to Fig. 1. The function $P_{2}(E)$ is defined for $E_{2}<E<E_{\mathrm{sep}}$, $P_{1}(E)$ for $E_{\mathrm{sep}}<E<E_{1}$, and $P_{3}(E)$ for $E>E_{\mathrm{sep}}$.

However, as the quasienergy states in regions 1 and 3 can have same energies, it is possible that the true quasienergy states are superpositions of the quasiclassical states from regions 1 and 3 due to quantum tunneling. So, the classical Fokker–Planck equation should be generalized to take this effect into account. It will be shown that hybridization of quasiclassical states from regions 1 and 3 strongly modifies the non–equilibrium statistics and kinetics of the system.

III The structure of quasienergy states

In this section, we will consider the model quasi–classically, although the results obtained here are valid beyond the quasiclassical approximation. Within the quasiclassical approach, the eigenstates of the quantum Hamiltonian correspond to a discrete set of trajectories which are obtained using the Bohr-Sommerfeld rule:

[TABLE]

The variables $q,p$ are the canonical coordinate and momentum defined by $a,a^{*}=\frac{q\pm ip}{\sqrt{2}}$. For the states from region 2, the value of quasienergy uniquely defines the classical trajectory. However, if a state has quasienergy $\epsilon_{\mathrm{sep}}<\epsilon<\epsilon_{1}$, it can lie either in the region 1 or region 3. Therefore if the quasienergies of some states obtained from the Bohr–Sommerfeld rule are close enough, the true eigenstates of the quantum Hamiltonian can be superpositions of the quasiclassical states due to quantum tunneling.

Let the set of quasiclassical states from regions 1 and 3 be $|n_{1}\rangle$, $|n_{3}\rangle$ with corresponding quasienergies $\epsilon_{n_{1}}$, $\epsilon_{n_{3}}$. If for some $n$, $n^{\prime}$ the quasienergies are almost equal, $\epsilon_{n_{1}}=\epsilon_{n^{\prime}_{3}}$, the quasiclassical states $|n_{1}\rangle$ and $|n^{\prime}_{3}\rangle$ form superpositions. Formally, for each $n,n^{\prime}$ it is possible to find the parameters $\beta$, $m$ such as $\epsilon_{n_{1}}=\epsilon_{n^{\prime}_{3}}$. In principle, these values of $\beta$ and $m$ could be different for different pairs $n,n^{\prime}$. However, a special feature of the Hamiltonian (1) is that the degeneracy of the quasiclassical states occurs exactly at integer values of $m$ for all $\beta$ in the region of bistability. Therefore at each integer value of $m$ all quasiclassical states from regions 1 and 3 can be grouped in pairs so that within each pair the values of quasienergy are equal. When $m$ is close to an integer, the states still can be grouped in pairs with close values of quasienergies. This is clearly seen from the exact diagonalization of the Hamiltonian (1). On the Fig. 2, the eigenvalues of the Hamiltonian (1) are shown as functions of $m$ at constant $\beta$. It can be seen that for $E_{\mathrm{sep}}<E<E_{1}$, there exist two families of lines which correspond to the states from regions 1 and 3. The anti–crossings of these lines indicate the degeneracy of the quasiclassical states. It is evident that all these anti–crossings occur exactly at integer values of $m$.

To put the statements of the previous paragraph on the theoretical ground, let us first consider the case $f=0$. At $f=0$, the Hamiltonian commutes with the number of excitation quanta operator $a^{\dagger}a$, and eigenstates have the form $|n\rangle$ with corresponding quasienergies

[TABLE]

At integer $m$, it is clear that all quasienergy levels split into pairs with same energy, as $\epsilon_{n}^{(0)}=\epsilon_{m-n}^{(0)}$.

Let us proceed to the case $f>0$. Let us consider the pair of states $|n\rangle$ and $|m-n\rangle$, $n<m/2$ which are degenerate at $f=0$. From the quasiclassical point of view, these states corresponds to circular trajectories on the phase plane with different radii. From the Bohr–Sommerfeld rule, it follows that the values of adiabatic invariant $\frac{1}{2\pi}\oint pdq$ for these states are $-n$ and $m-n$. Then let us switch on the external field adiabatically. After that, the quasienergies of the states change, but the values of the adiabatic invariant remain the same. The corresponding trajectories lie in regions 1 and 3 respectively unless $f$ is so large that they merge into a single trajectory from region 2. Therefore the quasienergies of the resulting states are $\epsilon_{1}(-n,\beta)$ and $\epsilon_{3}(m-n,\beta)$, where the energies are understood as functions of the adiabatic invariant.

For the considered model, it is possible to prove the identity

[TABLE]

or equivalently

[TABLE]

This happens because $n_{1}(E,f)$ and $n_{3}(E,f)$ have an analytic expression through the same elliptic integral with different contours (see Appendix A). With the same reasoning, it is also easy to verify that the classical periods of the trajectories from regions 1 and 3 with equal quasienergies are the same Dykman and Fistul (2005). Thus, in the quasiclassical limit the energies of the states with numbers $n$ and $m-n$ obtained from Bohr–Sommerfeld rule remain the same even at finite $f$. This supports the statement that the degeneracy of the quasiclassical states from regions 1 and 3 happens simultaneously for all states at integer values of $m$.

To find the eigenstates in the case when the quasiclassical states from regions 1 and 3 are degenerate, it is necessary to find the tunneling amplitude, which can also be obtained quasiclassically. For that, one should consider the motion of the system in the classically forbidden area. The quasiclassical tunneling exponent equals half of the action of a classical closed trajectory in imaginary time. The resulting amplitude is

[TABLE]

[TABLE]

where $p_{\mathrm{im}}$ is defined from the equation $H(q,ip_{\mathrm{im}})=\epsilon$. When $\beta\ll\beta_{\mathrm{crit}}$,

[TABLE]

For nearly degenerate states, when $E_{n}\approx E_{m-n}$, the tunneling exponent takes the form

[TABLE]

which is especially convenient for comparison with multi–photon transition amplitude obtained by perturbation theory.

The eigenstates with account for tunneling between the nearly–degenerate eigenstates can be found from an effective two–level Hamiltonian for two near–degenerate quasiclassical states. It has form

[TABLE]

where $\epsilon_{n,1}$ and $\epsilon_{n,3}$ are the quasienergies of the quasiclassical states $|n_{1}\rangle$, $|n_{3}\rangle$, and $\omega^{n}_{R}$ is the tunneling amplitude. When $m$ slightly deviates from an integer, $m=m_{0}+\delta m$, the difference of quasiclassical energies is

[TABLE]

The eigenstates of the effective Hamiltonian are

[TABLE]

where the coefficients $c_{n},s_{n}$ are defined by the ratio between $\omega_{R}$ and $\delta\epsilon_{n}$:

[TABLE]

Because of exponential dependence of $\omega_{R}$ on $E$, the coefficients $c_{n}$ and $s_{n}$ have a step–like dependence on $n$. When $\omega_{R}(E_{n})\ll\delta\epsilon_{n}$, $s_{n}\approx 0,\quad c_{n}\approx 1$, and the eigenstates correspond to distinct trajectories in the regions 1 and 3 of the phase portrait. On the contrary, when $\omega_{R}(E_{n})\gg\delta\epsilon_{n}$, $s_{n}\approx c_{n}\approx 1/\sqrt{2}$, and the states are very close to symmetric and antisymmetric superpositions of trajectories. These alternatives are separated by the critical value of quasienergy $E_{c}$ for which $\omega_{R}(E_{c})\approx\delta\epsilon(E_{c})$. Combining this with Eq. (15), one get the equation for $E_{c}$

[TABLE]

The resulting structure of quasienergy states is schematically depicted on Fig. 3. Below $E_{\mathrm{sep}}$, there are only states from classical region 2. Between $E_{\mathrm{sep}}$ and $E_{c}$, the states from regions 1 and 3 form superpositions, and above $E_{c}$, the states from regions 1 and 3 don’t hybridize.

IV The symmetry of the Hamiltonian

The quasiclassical arguments of the previous section are valid only in the leading order in the quasiclassicity parameter $m$. In particular, the quasienergy values of the system can be expanded in asymptotic series in $1/m$, and the Bohr–Sommerfeld rule gives only the leading term of these series. Therefore additional arguments are needed to explain the simultaneous anti–crossings of the quasienergy levels at integer $m$.

The rigorous proof can be given using the perturbation theory in $f$ for the quantum Hamiltonian (1). At $f=0$, the Hamiltonian commutes with $\hat{n}=a^{\dagger}a$, so the quasienergies of the Hamiltonian (1) are given by (7). At small $f$, one can use the perturbation theory to find the corrections to the energy of the state $|n\rangle$. For example, the second–order correction is

[TABLE]

The perturbative correction $\epsilon^{(2)}_{n}$ is symmetric with respect to replacement $n\to m-n$. This is in agreement with the quasiclassical arguments of Sec. III: it was shown that the changes of quasiclassical quasienergy due to adiabatic change of $f$ are the same for states $|n\rangle$ and $|m-n\rangle$. This statement is not only valid in the quasiclassical framework but also holds exactly in the second order of perturbation theory. Moreover, we obtained a rigorous proof that the same holds for higher orders $k$ of perturbation theory, $\epsilon_{n}^{(k)}=\epsilon_{m-n}^{(k)}$ up to order $m-2n$. The complete proof is given in the Appendix B. For integer $m$ one should utilize the degenerate perturbation theory which also takes into account multi–photon transitions between the states $|n\rangle$ and $|m-n\rangle$. These transitions occur only in order $m-2n$, and up to this order, the non–degenerate perturbation theory remains valid. Therefore the quasienergy splitting between the states $|n\rangle$ and $|m-n\rangle$ occurs only in the order $m-2n$.

V The generalization of Keldysh theory

for ionization of atoms in electromagnetic field

Splitting between the quasienergy states and the transitions between the regions of the phase space can be treated not only as tunneling between the regions of the phase space but also in $(m-2n)$–th order of perturbation theory in $f$, as mentioned in Sec. IV. Both approaches lead to just the same effects. Such behavior can be understood in the frame of Keldysh theory for ionization of atoms in electromagnetic field generalized for bistable systems with discrete spectrum. Transition amplitude depends on Keldysh parameter $\gamma_{K}\equiv T_{\mathrm{im}}(E)/T(E)$ where $T_{\mathrm{im}}$ is the <<tunneling time>> which is defined as $T_{\mathrm{im}}\equiv(\alpha/\Delta^{2})\partial_{E}S_{\mathrm{tunn}}$, where $S_{\mathrm{tunn}}(E)$ is the tunneling action defined by (12). The <<tunneling time>> has the meaning of the half–period of motion along the closed trajectory in the imaginary time. From the expression (13), we obtain that $T_{\mathrm{im}}(E)\sim\Delta^{-1}\ln{\frac{1}{\beta}}$. The period in the real time is always $\sim\Delta^{-1}$ unless the quasienergy is close to $E_{\mathrm{sep}}$. So, $\gamma_{K}\sim\ln{\frac{1}{\beta}}\gg 1$ for $\beta\ll\beta_{\mathrm{crit}}$. In this limit, the tunneling amplitude coincides with $(m-2n)$–order perturbation theory multi–photon transition amplitude. Indeed, as it was shown in Maslova et al. (2019), the multi–photon transition amplitude reads

[TABLE]

where $\hat{V}=f(\hat{a}+\hat{a}^{\dagger})$. This coincides with the expression $\omega_{R}(E)\propto\Delta e^{-S_{\mathrm{tunn}}(E)}$ where $S_{\mathrm{tunn}}(E)$ is defined by (13). According to both formulas, $\omega_{R}\sim\Delta\beta^{\frac{m}{2}-n}$. Moreover, not only the power–law dependence on $\beta$ coincides but also the numerical coefficient which is accurately derived in the Appendix C. Therefore tunneling and multi–photon transitions are just the same effects. Also let us note that for bistable driven systems the Keldysh parameter $\gamma_{K}$ logarithmically depends on the external field amplitude $f$ whereas in the case of ionization of atoms by strong electromagnetic field $\gamma_{K}\sim f^{-1}$.

VI The effect of degeneracy on kinetics

The fact that the eigenstates of the quantum driven nonlinear oscillator can be superpositions of quasiclassical states from regions $1$ and $3$ has a strong effect on kinetics described by the rate equation. When each state can be uniquely attributed to a single region of the phase space, in the limit of large number of excitation quanta the rate equation transforms to the classical Fokker–Planck equation in quasienergy representation. However, this is not the case when $m$ is close to an integer. The actual eigenstates are the superpositions of the quasiclassical states from regions $1$ and $3$, which breaks the crucial assumption under which the Fokker–Planck equation is derived.

However, even when $m$ is close to integer, there exists a quasiclassical limit of the rate equation corresponding to large $m$ which has also the form of the classical Fokker–Planck equation in the quasienergy space. It is obtained under assumption that the occupations $P^{+}_{n}$ and $P^{-}_{n}$ of states $|n_{+}\rangle$ and $|n_{-}\rangle$ (see Eq. (16)) still depend smoothly on $n$. The rate equation in terms of $P^{s}_{n}$, $s=\pm$, takes form

[TABLE]

where the transition probabilities $w_{nn^{\prime}}^{ss^{\prime}}$ between the states $|n_{s}\rangle$, $|n^{\prime}_{s^{\prime}}\rangle$ are defined by general formula (4), and the states $|n_{s}\rangle$ are defined by (16). In the equation (21), it is possible to perform a gradient expansion of $P_{n}^{\pm}$. To perform such a procedure, it is convenient to rewrite the equation (21) in a slightly different form:

[TABLE]

with the newly defined coefficients $W^{+}_{nn^{\prime}}=w^{++}_{nn^{\prime}}+w^{-+}_{nn^{\prime}}$ and $W^{-}_{nn^{\prime}}=w^{--}_{nn^{\prime}}+w^{+-}_{nn^{\prime}}$. They are expressed via matrix elements between the states from regions 1 and 3 as

[TABLE]

Using the rate equations in the form (22), it is easy to understand the structure of the quasiclassical limit of equations (22). Following the derivation of the Fokker–Planck equation which was in detail described in Maslova et al. (2019), it is clear that after gradient expansion the terms $\sum_{n^{\prime}}W^{\pm}_{nn^{\prime}}P_{n^{\prime}}^{\pm}-W^{\pm}_{n^{\prime}n}P_{n}^{\pm}$ transform to the expressions $\frac{1}{T}\frac{\partial}{\partial E}\left[\vartheta K^{\pm}P^{\pm}+QD^{\pm}\frac{\partial P^{\pm}}{\partial E}\right]$ where $K^{\pm}(E)$ and $D^{\pm}(E)$ describe the drift and the diffusion correspondingly in quasienergy space. The terms $\pm(w^{+-}_{nn^{\prime}}P_{n^{\prime}}^{-}-w^{-+}_{nn^{\prime}}P_{n^{\prime}}^{+})$ describe tunneling transitions between the states $|n_{+}\rangle$ and $|n_{-}\rangle$. Thus, the whole system of equations takes the form

[TABLE]

where

[TABLE]

and $\Lambda_{\mathrm{tunn}}^{+(-)}$ are the terms responsible for tunneling:

[TABLE]

The system (24) is exactly the Fokker–Planck equation with tunneling term from Maslova et al. (2019). The features of the equations (24) can be understood from the structure of the eigenstates which was described in Sec. III, see Eq. (16) and Fig. 3. As was mentioned in Sec. III, the hybridization of the quasiclassical states from different regions of the phase portrait is strong for quasienergies below $E_{c}$ and it is very small for quasienergies above $E_{c}$. This defines the behavior of the coefficients $K^{\pm}(E)$, $D^{\pm}(E)$ and terms $\Lambda_{\mathrm{tunn}}^{\pm}(E)$ in (24) which is different for $E<E_{c}$ and $E>E_{c}$. For $E>E_{c}$ there is almost no hybridization of states from regions 1 and 3. Therefore, the drift and diffusion coefficients of the Fokker–Planck equation $K^{\pm}=K_{1,3}$, $D^{\pm}=D_{1,3}$ and the tunneling transition rate is small: $\Lambda^{\pm}_{\mathrm{tunn}}T\ll 1$. On the contrary, for $E>E_{c}$ the actual eigenstates are symmetric and antisymmetric superpositions of the quasiclassical states from different regions. Using the expressions (23) and keeping in mind that in the considered case $s_{n}\approx c_{n}\approx 1/\sqrt{2}$, it is obvious that the drift and diffusion coefficients for $|n_{\pm}\rangle$ states are $K^{+}\approx K^{-}\approx\frac{1}{2}(K_{1}+K_{3})$, $D^{+}\approx D^{-}\approx\frac{1}{2}(D_{1}+D_{3})$, and the tunneling rates are large and almost equal: $\Lambda_{\mathrm{tunn}}^{+}T\approx\Lambda_{\mathrm{tunn}}^{-}T\sim 1$.

From the considerations of the previous paragraph, it is easy to obtain the stationary distribution. As for $E<E_{c}$ the drift and diffusion coefficients and tunneling transition rates are almost equal for states $|n_{\pm}\rangle$ and the tunneling transition rate is large, the stationary probability densities $P^{+}$ and $P^{-}$ are almost equal too, $P^{+}\approx P^{-}\approx\bar{P}$. Thus, the stationary probability density $\bar{P}$ is obtained from the condition of zero flow $J^{\pm}(E)$ and is given by the expression

[TABLE]

For $E_{c}<E<E_{1}$, the tunneling term in (24) is small, and the stationary distributions $P^{+}$($P^{-}$) can be obtained by perturbation theory in $\Lambda_{\mathrm{tunn}}^{\pm}$ as it was done in Maslova et al. (2019).

Let us define

[TABLE]

where $P^{+(-)}_{0}(E)$ are stationary distributions without tunneling term and $\delta P^{+(-)}$ are the first–order corrections caused by tunneling terms. Then

[TABLE]

The coefficients in (29) are defined by the continuity condition at $E=E_{c}$. The tunneling corrections have the form

[TABLE]

where $\Lambda_{\mathrm{tunn}}^{\pm}$ are expressed by equations (26) with $P_{0}^{\pm}(E)$.

On Fig. 4, the distribution functions obtained from analytical formulas (29), (30) are compared with those obtained by numerical solution of the equation (21), (22) for different values of $\delta m$. It is evident that the critical quasienergy $E_{c}$ shifts towards $E_{\mathrm{sep}}$ with increasing $\delta m$, according to the Eq. (18). Also it can be seen that for exactly degenerate quasienergy levels $\delta m=0$, the system remains close to the stable state 2 which is squeezed Maslova et al. (2019). In this case, the states corresponding to the regions 1 and 3 are equally occupied and the occupation probabilities are exponentially small.

The ratio between the probability densities in two stable stationary states equals

[TABLE]

The tunneling correction for $E>E_{c}$ doesn’t lead to any qualitative effects because it is of order $\alpha/\Delta$ comparing to $P_{0}^{\pm}$.

On Fig. 5, the analytical formula (31) is compared with the numerical result for the dependence of $P_{1}(E_{1})/P_{2}(E_{2})$ on $\delta m$. The analytical formula fits the numerical result quite well. However, the discreteness of the quasienergy levels manifests itself in the smooth steps in the dependence $P_{1}(E_{1})/P_{2}(E_{2})$ on $\delta m$ which are not reproduced by (31). These steps can be explained by the fact that the crossover energy $\epsilon_{c}$ can take only discrete values. Thus, for

[TABLE]

the effective position of $E_{c}$ remains the same. When $\delta m\approx e^{-S_{\mathrm{tunn}}(E_{n+1})}$, the value of $E_{c}$ abruptly changes from $E_{n}$ to $E_{n+1}$. This explains the presence of steps on Fig. 5. The width of the steps in the logarithmic scale is $S_{\mathrm{tunn}}(E_{n})-S_{\mathrm{tunn}}(E_{n+1})\approx 2\pi T_{\mathrm{im}}(E_{n})/T(E_{n})$. In the latter expression, we recognize the previously defined Keldysh parameter $\gamma_{K}$. For $n$ much smaller than $m$ and $\beta\ll\beta_{\mathrm{crit}}$, $\gamma_{K}\sim\ln{\frac{1}{\beta}}$.

On the inset in Fig. 5, the dependence of probabilities ratio $P_{1}(E_{1})/P_{2}(E_{2})$ on $m$ is shown in linear scale for $\sqrt{\beta/\beta_{\mathrm{crit}}}=0.2$. For non–degenerate case, when $m$ is far from an integer, the state with smaller amplitude and quasienergy $E_{1}$ is the most probable. However, when $m$ becomes close to an integer, the occupation of the state with quasienergy $E_{1}$ abruptly drops, and the value at the minima is exponentially small for large $m$ (see the dips on the inset of Fig. 5). The width of the dips in linear scale is also exponentially small.

VII Conclusions

We considered the non–equilibrium statistics and kinetics of the model of resonantly driven quantum nonlinear oscillator interacting with dissipative environment. We found out that the non–equilibrium statistics and kinetics are strongly modified when quasienergy states are nearly degenerate which occurs at integer or half–integer detuning–nonlinearity ratio. In particular, the occupation of the classical stable state with smaller amplitude is strongly reduced. So, in the case of exactly degenerate quasienergy levels the system always occupies the state with higher amplitude which is squeezed.

The coefficients of the Fokker–Planck equation which describes the kinetics in the quasiclassical limit are very sensitive to the structure of eigenstates of the system’s Hamiltonian. We found out that in the case of integer or half–integer detuning–nonlinearity ratio, which corresponds to the exact multi–photon resonance between the genuine energy levels of the unperturbed nonlinear oscillator, the quasienergy states from different regions of the phase space simultaneously hybridize and form symmetric and antisymmetric superpositions. This fact can be proven by applying a special symmetry transformation to the Hamiltonian. Also we revealed that when the quasienergy levels of the system are nearly degenerate, a new important critical quasienergy parameter $\epsilon_{c}$ emerges. Below $\epsilon_{c}$, all quasienergy states are superpositions of the quasiclassical states from regions 1 and 3, and above $\epsilon_{c}$, the quasienergy states correspond to either region 1 or 3. We found out that the coefficients of the Fokker–Planck equation which describes the quasiclassical kinetics of the oscillator in almost–degenerate case have different behavior above and below $\epsilon_{c}$. In particular, tunneling term is large below $\epsilon_{c}$ and exponentially small above $\epsilon_{c}$. Also the drift and diffusion coefficients are not affected by tunneling above $\epsilon_{c}$ whereas below $\epsilon_{c}$ they are strongly modified. The distribution functions and the ratio between occupations of the classical stable states calculated analytically fit well the numerical results.

We generalized Keldysh theory for ionization of atoms in electromagnetic field for bistable systems. It was demonstrated that Keldysh parameter defined as the ratio of <<tunneling time>> to the quasiclassical period of motion along the phase trajectory is large in the bistability region for external field intensity smaller than the critical value. So the multi–photon transition and tunneling between different regions of the phase space can be treated as the same effects. This fact was proved by direct calculation of transition amplitude using both tunneling and perturbation theory approach. Also we revealed that the Keldysh parameter for the considered system depends logarithmically on the amplitude $f$ of the external field. On the contrast, in the case of multi–photon ionization of atoms the Keldysh parameter is inversely proportional to the amplitude of the external field.

Acknowledgements.

This work was supported by RFBR grants 19–02–000–87a and 18–29–20032mk and by a grant of the Foundation for the Advancement of Theoretical Physics and Mathematics ’’Basis’’.

Appendix A The coefficients of the classical Fokker–Planck equation

The coefficients of the classical FPE are defined as line integrals along the classical trajectories of the nonlinear oscillator:

[TABLE]

The classical trajectory is a contour line of the classical Hamiltonian

[TABLE]

where $a=\frac{q+ip}{\sqrt{2}}$, and $q,p$ is a pair of canonically conjugate variables. The index $i$ denotes the region of the phase space according to Fig. 1. The coefficient $K_{i}(E)$ is proportional to the adiabatic invariant of the trajectory defined as

[TABLE]

The expressions (33) and (6) can be rewritten as two–dimensional integrals in $p,q$ plane with a Dirac delta function as in expression for $T_{i}$. Then, it is convenient to use variables $q,t$ instead of $q,p$ where $t=q^{2}+p^{2}$. Then the coefficients are transformed to one–dimensional integrals by $t$. Now let us focus on the expression for adiabatic invariant:

[TABLE]

Using this expression, we will prove the identity of Eq. (9). The contour of integration in (36) depends on the region of the phase space in which the trajectory lies. In the range of energies corresponding to the region 2 of the phase space, the polynomial has only two real roots, and the contour of integration encloses them. In the range of energies corresponding to the region $1$, there are 4 real roots: $t_{1}<t_{2}<t_{3}<t_{4}$ (see Fig. 6). The range $t_{1}<t<t_{2}$ ($t_{3}<t<t_{4}$) corresponds to the trajectories from the region 1 (3). Thus, the contour of integration for $n_{1}$ ($n_{3}$) encloses $t_{1}$ and $t_{2}$ ($t_{3}$ and $t_{4}$).

By deformation of the contour, it is easy to show that $n_{3}(\epsilon,f)-n_{1}(\epsilon,f)$ is expressed through residue of integrand (36) on infinity. Expanding the integrand of (36) in $t^{-1}$, one gets the desired identity

[TABLE]

Appendix B The proof of the symmetry of the perturbative corrections

Here we prove that the perturbation theory corrections in $f$ to the quasienergies $\epsilon_{n}^{(0)}$ of the eigenstate with $n$ excitation quanta of t he model (1) at $f=0$ are symmetric with respect to replacement $n\to m-n$. This fact was mentioned in Risken and Vogel (1988) and Dykman and Fistul (2005), but the authors didn’t give any proof to this fact.

For several low–order corrections, this can be verified by straightforward calculation, as for the second–order correction (19):

[TABLE]

However, it is necessary to clarify what the expression $\epsilon_{m-n}^{(k)}$ means for non–integer $m$ because the corrections $\epsilon_{n}^{(k)}$ are defined only for integer $n$ which has the meaning of the number of excitation quanta. Thus, for non–integer $m$ the identity $\epsilon_{n}^{(k)}=\epsilon_{m-n}^{(k)}$ holds only for formal expressions. Up to now, we didn’t give any meaning to $\epsilon_{\nu}^{(k)}$ for non–integer $\nu$ except as analytic continuation of perturbation theory formulas.

However, it is in fact possible to give direct meaning to $\epsilon_{\nu}$ and $\epsilon_{\nu}^{(k)}$ at non–integer $\nu$. For that, we should formally assume that the Hamiltonian (1) acts on the space of all possible real <<numbers of excitation quanta>> $\nu$ with operators $a$, $a^{\dagger}$ defined as follows:

[TABLE]

Then, for $f=0$ each state $|\nu\rangle$ is an eigenstate with energy $\epsilon_{\nu}^{(0)}=\nu(m-\nu)$. For $f\neq 0$, they become coupled with $|\nu\pm 1\rangle,|\nu\pm 2\rangle,\dots,|\nu\pm k\rangle\dots$. If none of the states $|\nu\pm k\rangle$ are degenerate with $\nu$ (equivalently, $2\nu-m$ is non–integer), the amplitudes of $|\nu\pm k\rangle$ remain small at small $f$, and it is possible to define continuous $f$–dependent energy $\epsilon_{\nu}(f)$. This is the energy of the eigenstate which evolves from $|\nu\rangle$ after adiabatic switching of the perturbation. For integer $2\nu-m$, the energy $\epsilon_{\nu}(f)$ can’t be defined that way because of degeneracy between $|\nu\rangle$ and $|m-\nu\rangle$.

The series of perturbation theory for $\epsilon_{\nu}$ in the cases of integer and non–integer $\nu$ are completely identical because of definition (39). Thus, the claim that $\epsilon^{(k)}_{n}=\epsilon^{(k)}_{m-n}$ follows from even more general statement $\epsilon_{\nu}=\epsilon_{m-\nu}$.

We prove the identity $\epsilon_{\nu}=\epsilon_{m-\nu}$ in several steps. First, it is obvious from the previous considerations that $\epsilon_{\nu}(f)$ is an eigenvalue of the operator

[TABLE]

which corresponds to the state $|\nu\rangle$. Analogously, $\epsilon_{m-\nu}(f)$ arises from the operator $\mathcal{H}_{m-\nu}$. For convenience in the later discussion, we change the numeration of basis vectors in $\mathcal{H}_{m-\nu}$ so that $|\sigma\rangle$ becomes $|m-\sigma\rangle$. After such relabeling,

[TABLE]

Both $\mathcal{H}_{\nu}$ and $\mathcal{H}_{m-\nu}$ act on a single space with a set of basis vectors $|\sigma\rangle$ with such $\sigma$ that $\sigma-\nu$ is integer. We should emphasize that $\mathcal{H}_{\nu}$ and $\mathcal{H}_{m-\nu}$ are substantially different and could not be transformed to each other by any permutation of eigenvectors.

However, there exists a nontrivial linear operator $T$ which transforms $\mathcal{H}_{\nu}$ to $\mathcal{H}_{m-\nu}$:

[TABLE]

It has the form

[TABLE]

where

[TABLE]

The identities (42), (43), (44) are checked by direct calculation.

The existence of the operator $\mathcal{T}$ is possible only because of special form of $\epsilon_{\nu}=\frac{\alpha}{2}\nu(m-\nu)$. For any other dependence of $\epsilon_{\nu}$ on $\nu$, no analogous operator can be found this way. So, the symmetry property expressed by $\mathcal{T}$ is a special feature of Kerr–like nonlinearity.

The equivalence of Hamiltonians $\mathcal{H}_{\nu}$ and $\mathcal{H}_{m-\nu}$ proves that the energies $\epsilon_{\nu}(f)$ and $\epsilon_{m-\nu}(f)$ are equal when $2\nu-m$ is not integer. However, we are interested in the case of integer $m$ and integer numbers of excitation quanta. For this case, one should utilize degenerate perturbation theory to find the energies. Nevertheless, the corrections to energies of $|n\rangle$ and $|m-n\rangle$ obtained by degenerate perturbation theory are just the same as in non–degenerate perturbation theory up to the order $f^{|m-2n|}$. This happens because the leading contribution to composite matrix element (multi–photon Rabi frequency) between $|n\rangle$ and $|m-n\rangle$ is a product of $|m-2n|$ matrix elements of the perturbation $\hat{V}$ (see (20), (45)). For series of non–degenerate perturbation theory the identity for $k$–th order corrections $\epsilon_{n}^{(k)}=\epsilon_{m-n}^{(k)}$ holds even for integer $m$ and $n$, if $k<2|m-2n|$. For $k\geq 2|m-2n|$, the corrections of non–degenerate perturbation theory don’t make sense because of a zero in denominator, which is a manifestation of degeneracy. This means that the degeneracy of $|n\rangle$ and $|m-n\rangle$ is lifted only in the order $|m-2n|$, and the energy splitting happens only due to multi–photon Rabi oscillations: $\Delta\epsilon_{n,m-n}=2\omega_{R}^{n,m-n}+o(f^{|m-2n|})$.

Appendix C The identity of tunneling splitting and multi–photon transition amplitude

The multi–photon Rabi splitting between the quasienergy states $|n\rangle$ and $|m-n\rangle$ is given by formula

[TABLE]

At large values of $n$, $m$, it is possible to approximate the factorials using the Stirling formula. Then one gets the following expression for $\omega_{R}^{n,m-n}$:

[TABLE]

where $r=2n/m$. In this form, it is easy to compare it with tunneling splitting given by Eq. (10) and Eq. (11). At small $\beta$, the tunneling action (11) can be approximated as

[TABLE]

For small external force the quasienergy $E$ is related with number of excitation quanta by formula $E=\frac{2n}{m}(\frac{n}{m}-1)=r^{2}/2-r$. Evaluating the integral in (47) and substituting the expression for $E$ via $r$, it is easy to obtain that in current approximations $S_{\mathrm{tunn}}=\ln{(\omega_{R}^{n,m-n}/\Delta)}$. Thus, at $\beta\ll\beta_{\mathrm{crit}}$ perturbation theory is consistent with tunneling approach.