Effects of tunneling and multiphoton transitions on squeezed states generation in bistable driven systems

TL;DR

This paper investigates how tunneling and multiphoton transitions influence the generation of squeezed states in bistable driven nonlinear oscillators, revealing enhanced squeezing and relaxation effects due to quantum phenomena.

Contribution

It demonstrates the impact of tunneling and multiphoton transitions on squeezing and relaxation in bistable systems, linking classical and quantum descriptions through Fokker-Planck and master equations.

Findings

Tunneling enhances slow relaxation rates between stable states.

Quasienergy states in bistability regions are squeezed with uncertainty below quantum limit.

Tunneling effects can significantly improve the generation of squeezed states.

Abstract

Bistability of nonlinear resonantly driven oscillator in the presence of external noise is analyzed by means of classical Fokker-Planck equation in quasienergy space with account for tunneling effects and by quantum master equation in quasienergy states representation. Two time scales responsible for different stages of bistable system relaxation have been obtained. We found out that the slow relaxation rate caused by fluctuation--induced transitions between different stable states can be enhanced by several orders due to tunneling effects. It was also revealed that tunneling between nearly degenerate quasienergy states and resonant multiphoton transitions between the genuine eigenstates of the nonlinear oscillator are just the similar effects. It was demonstrated that the quasienergy states in the bistability region corresponding to higher amplitude are squeezed. The degree of squeezing is determined by the ratio between nonlinearity and detuning, so the uncertainty of one quadrature can be considerably smaller than the quantum limit. It was found out that tunneling effects can enhance the generation of output oscillator squeezed states. It is shown that 1D Fokker--Planck equation is a quasiclassical limit of a quantum master equation.

Figures

Equations

$$H_{\mathrm{eff}}=-\Delta |a{|}^2+\frac{\alpha }{2}|a{|}^4-f(a+a^{\ast }),$$

$$\mathcal{H}=-|A{|}^2+\frac{1}{2}|A{|}^4+\sqrt{\beta }(A+A^{\ast }),$$

$$i\frac{\partial A}{\partial \tau }=-A+A|A{|}^2+\sqrt{\beta }.$$

$$\begin{array}{c}i\frac{\partial A}{\partial \tau }=-A(1+i\vartheta )+A|A{|}^2+\sqrt{\beta }+\xi ,\\ ⟨\xi (\tau ){\xi }^{\ast }({\tau }^{\prime })⟩=Q\delta (\tau -{\tau }^{\prime }),\vartheta \equiv \frac{\gamma }{\Delta }.\end{array}$$

$$\frac{\partial \mathcal{P}}{\partial t}=\frac{\partial }{\partial A}\left(i\mathcal{P}\frac{\partial H}{\partial A^{\ast }}+\vartheta A\mathcal{P}+\frac{Q}{2}\frac{\partial \mathcal{P}}{\partial A^{\ast }}\right)+c.c.$$

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

$$\begin{array}{c}{P}_1(E_{\mathrm{sep}})=P_2(E_{\mathrm{sep}})=P_3(E_{\mathrm{sep}}),\\ J_2(E_{\mathrm{sep}})=J_1(E_{\mathrm{sep}})+J_3(E_{\mathrm{sep}}).\end{array}$$

$$P_i^{\mathrm{st}}(E)=C{e}^{S_i-S_i(E)},$$

$$\begin{array}{c}{S}_i(E)\equiv \frac{\vartheta }{Q}{\int }_{E_i}^E\frac{K_i(E)}{D_i(E)}dE,S_i\equiv S_i(E_{\mathrm{sep}}).\end{array}$$

$$P_i(t,E)=\sum _{\lambda }{P}_i^{\lambda }(E)e^{-\lambda t},$$

$$-\lambda P_i^{\lambda }=\frac{1}{T_i(E)}\frac{\partial }{\partial E}\left[\vartheta K_i(E)+Q{D}_i(E)\frac{\partial }{\partial E}\right]P_i^{\lambda }.$$

$$\begin{array}{c}{P}_i^{\lambda }(E)=P_i^{\mathrm{st}}(E)\left[1+\frac{1}{Q}{\int }_{E_{\mathrm{sep}}}^E\frac{{\Phi }_i(E^{\prime })d{E}^{\prime }}{D_i(E^{\prime })P_i^{\mathrm{st}}(E^{\prime })}\right],\\ {\Phi }_i(E)=-\lambda {\int }_{E_i}^{E}d{E}^{\prime }{T}_i(E^{\prime })P^{\mathrm{st}}(E^{\prime }),\end{array}$$

$$\lambda =\frac{{\vartheta }^2}{Q}\cdot \frac{K_1(E_{\mathrm{sep}})K_2(E_{\mathrm{sep}})}{K_2(E_{\mathrm{sep}})-K_1(E_{\mathrm{sep}})}\left[\frac{e^{-S_2}}{D_2^{\prime }(E_2)}-\frac{e^{-S_1}}{D_1^{\prime }(E_1)}\right].$$

$$\begin{array}{c}\frac{\partial P_{(1,3)}}{\partial t}=\frac{1}{T(E)}\frac{\partial J_{(1,3)}}{\partial E}+{\lambda }_T(E)(P_{(3,1)}-P_{(1,3)}).\end{array}$$

$${\lambda }_T/\vartheta \sim e^{-\frac{2\Delta }{\hslash \alpha }{S}_{\mathrm{tunn}}(E)}.$$

$$\begin{array}{c}{P}_{(1,3)}^T(E)=P_{(1,3)}^{\mathrm{st}}(E)\left[1+\frac{1}{Q}{\int }_{E_{\mathrm{sep}}}^E\frac{{\tilde{\Phi }}_{(1,3)}(E^{\prime })d{E}^{\prime }}{D(E^{\prime })P_{(1,3)}^{\mathrm{st}}}\right],\\ {\tilde{\Phi }}_{(1,3)}(E^{\prime })={\int }_{E_{(1,\mathrm{sep})}}^{E^{\prime }}d{E}^{\prime \prime }T(E^{\prime \prime }){\lambda }_T(E^{\prime \prime })(P_{(1,3)}^{\mathrm{st}}-P_{(3,1)}^{\mathrm{st}}).\end{array}$$

$$S_{\mathrm{total}}=\frac{2\Delta }{\hslash \alpha }{S}_{\mathrm{tunn}}(E)+\frac{\vartheta }{Q}{S}_i(E).$$

$$\mathcal{H}=-\left(\frac{q^2+p^2}{2}\right)+\frac{1}{2}{\left(\frac{q^2+p^2}{2}\right)}^2-\sqrt{2\beta }q.$$

$$S_{\mathrm{tunn}}={\int }_{q_1}^{q_c}|p_1|dq+{\int }_{q_c}^{q_2}|p_3|dq.$$

$$p_{1,3}^2=2-q^2\pm 2\sqrt{1+2E-2\sqrt{2\beta }q}.$$

$$S_{\mathrm{tunn}}={\int }_{q_1}^{q_2}\mathrm{acosh}\left\{\frac{E+\frac{x^2}{2}-\frac{x^4}{8}}{\sqrt{2\beta }x}\right\}xdx.$$

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

$$S_{\mathrm{tunn}}(E_1)=\ln \frac{1}{\beta }+O(1).$$

$${\lambda }_T\propto \frac{{\Delta }^2}{\gamma }{\beta }^{\frac{2\Delta }{\hslash \alpha }}.$$

$${\lambda }_T\propto \gamma {\beta }^{\frac{2\Delta }{\hslash \alpha }}.$$

$$A_{k,m-k}=\frac{\Delta }{m}{\left(\frac{\beta m^3}{2}\right)}^{\frac{m}{2}-k}\frac{\sqrt{(m-k)!}}{(m-2k-1){!}^2\sqrt{k!}}.$$

$$|A_{0,m}{|}^2\propto {\Delta }^2{\beta }^m.$$

$$\frac{P_2^{\mathrm{st}}(E_2)}{P_1^{\mathrm{st}}(E_1)}=e^{S_2}\left(e^{-S_1}+\frac{{\lambda }_T(E_1)Q}{{\vartheta }^2}\frac{D_1^{\prime }(E_1)}{K_1(E_{\mathrm{sep}})}\right).$$

$${\lambda }_{\mathrm{total}}={\lambda }_{\mathrm{noise}}+{\lambda }_T(E_1).$$

$$\begin{array}{c}{\hat{H}}_0=-\Delta {\hat{a}}^{†}\hat{a}+\frac{\alpha }{2}{\hat{a}}^{†}{\hat{a}}^{†}\hat{a}\hat{a}+f(\hat{a}+{\hat{a}}^{†}),\\ [\hat{a},{\hat{a}}^{†}]=1.\end{array}$$

Full text

I Introduction

Complex systems with two or more stable states appear in many fields of science from biology and chemistry to quantum optics and electronics Innes et al. (2013) Semenov et al. (2016) Hoang et al. (2018) Ma et al. (2018).

Ability to control and manipulate these complex systems relies on one’s knowledge of their stable states, the extent of their robustness with respect to environmental fluctuations, and on ability to make specific perturbations inducing transitions between these states.

First, one needs to understand the behavior of bistable systems interacting with the environment. Bistable systems in optics and electronics are widely used as switching elements in communications systems, basic elements of memory devices, logic gates, optical turnstiles, etc. So, the investigation of fluctuation–induced transitions between different stable states is crucial to improve the stability of optical and electronic devices and to control their switching rates.

Due to unprecedented miniaturization of optical and electronics devices quantum effects became very important for their operation Kasprzak et al. (2010) Albert et al. (2013). It is impossible to study bistability without accounting for quantum effects. Thus, it is necessary to trace the correspondence between the classical and quantum regimes of the system Kasprzak et al. (2013). In the quantum regime it is important to understand whether quantum fluctuations impose a fundamental limit on stability of optical and electronic devices.

It is well established that in many nonlinear optical and electronic interface systems there exist a set of quantum states — squeezed states — which have less uncertainty in one quadrature than a coherent state Walls (1983). Generation of squeezed states is a key for implementation of quantum information protocols and for stability enhancement of quantum optics devices Savage and Walls (1986). Bistable quantum optics systems are promising candidates to realize the squeezed states. Recently the squeezed exciton–polariton field has been observed in pillar–shaped semiconductor microcavities in the bistable regime near the critical point of the bistable curve Boulier et al. (2014).

A driven nonlinear oscillator interacting with a thermal bath is the minimal model describing fluctuation–induced transitions in bistable systems out of equilibrium. The dynamics of various microcavities coupled with nonlinear media and coherently driven by an external field including exciton–polaritons in semiconductor microresonators with external pumping can exhibit a bistable behavior and can be described by the model of a driven nonlinear oscillator. Recent experiments demonstrated that as external coherent pumping is increased the occupied exciton–polariton mode shows strong sudden jumps from one state to another. Such behavior is caused by the fluctuation–induced transitions between the stationary states. These transitions could also lead to decrease of the hysteresis area of an internal microcavity field under the S–shaped response curve with respect to the external pumping Rodriguez et al. (2017).

Another experimental realization which can be analyzed using the nonlinear oscillator model is a mesoscopic Josephson junction array resonator Muppalla et al. (2018). In such a device, the anharmonicity can be of the same order as the linewidth, and the dynamics of bistability has been experimentally measured by observing the jumps between different stable states. It was shown experimentally that the switching rate strongly depends on the pumping intensity.

In addition, the model of a driven nonlinear oscillator is applicable to highly excited molecular vibration modes in the presence of anharmonicity.

The model of a driven nonlinear oscillator has been extensively theoretically studied since 1980’s. However, fluctuation–induced transitions between two stable states were traditionally analyzed using the classical 1D Fokker–Planck equation (FPE) in quasienergy space without accounting for quantum tunneling Vogel and Risken (1988), Maslova (1986), Dmitriev and D’yakonov (1986).

The ultra–quantum limit of dispersive bistability was analyzed by Drummond and Walls Drummond and Walls (1980), where the kinetic equation for generalized Glauber function was solved analytically for the case of zero bath temperature. The same model was analyzed numerically using the technique of quantum master equation Risken et al. (1987).

Nevertheless, there has been no detailed analysis of kinetics of the nonlinear driven oscillator allowing one to trace the transition between the classical and quantum descriptions of this system. Moreover, the structure of quasienergy states of a driven quantum nonlinear oscillator and the influence of their degeneracy Dykman and Fistul (2005) on kinetics still remains a relatively unexplored area of research. In order to understand rich physical properties of bistable systems one could start by considering the minimal model of a driven nonlinear oscillator.

In present work, we derive the quasiclassical kinetic equations taking into account the tunneling effects. These equations are a quasiclassical limit of the quantum master equation for the density matrix of a quantum driven nonlinear oscillator. We show that in the quasiclassical limit, tunneling transitions reduce the threshold value of intensity of the external field responsible for switching between the most probable states of the system. We also show that tunneling between trajectories in different regions of the phase space and multi–photon resonant transitions between the states of the nonlinear oscillator are the same effects. In the quantum case, we explore the structure of eigenstates and show that the quasienergy states corresponding to the higher amplitude stable state are squeezed, and the uncertainty in one of the quadratures can be much lower than the usual quantum limit.

II Classical bistability

II.1 The basic model

We consider a model system consisting of a single oscillator mode with Kerr–like nonlinearity excited by a resonant field. Its key feature is the bistability in a certain range of external pumping intensity: the presence of two different classical stable states.

The effective Hamiltonian for such a model is Maslova (1986), Vogel and Risken (1990), Maslova et al. (2007)

[TABLE]

where $a,a^{*}$ are slowly varying amplitudes of the internal oscillator field; $\Delta\equiv\Omega-\omega_{0}$ is the detuning between the external field frequency $\Omega$ and the frequency of the resonance $\omega_{0}$; $\alpha$ is the anharmonicity parameter; $f$ is the interaction strength with external field (proportional to its amplitude). Such a model can arise for various systems in the rotating–wave approximation such as microcavity with a nonlinear medium coherently driven by an external field. For example, this effective Hamiltonian can be derived for the Janes–Cunnings model after adiabatically excluding the atomic variables. It also describes the microcavity exciton–polaritons driven by an external field as well as strongly excited vibration modes of molecules in the presence of an external resonant field. Here we use the normalized field amplitude, $A\equiv a\sqrt{\alpha/\Delta}$, and a dimensionless time, $\tau=\Delta t$. The only dimensionless parameter which governs the system dynamics is $\beta\equiv\alpha f^{2}/\Delta^{3}$.

In dimensionless variables the parameter $\beta$ 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. Note that the original Dicke model deals with collective superradiance of the system of quantum two–level emitters interacting with the cavity field. However, as it was shown in Refs. Il’inskii and Maslova (1988) and Kocharovsky et al. (2017), a superradiance pulse can also arise in a classical system of nonlinear oscillators coupled to the cavity field due to rephasing processes.

In terms of new variables the dimensionless Hamiltonian is given by $\mathcal{H}=(\alpha/\Delta^{2})H_{\mathrm{eff}}$:

[TABLE]

while the equation of motion reads

[TABLE]

The classical phase trajectories of the nonlinear oscillator in the plane $(A,A^{*})$ are the contour lines of the classical Hamiltonian function (1) [Fig. 1(a)]. Let us focus on the structure of (1) as the function of two variables, $\operatorname{Re}A$ and $\operatorname{Im}A$.

At $\beta=0$, the function has a shape of the Mexican hat potential. It is radially symmetric, and its contour lines are concentric circles. At nonzero $\beta$, $0<\beta<\beta_{\mathrm{crit}}=4/27$, the hat is deformed, as shown in Fig. 1(b). Instead of infinitely many local minima, two extrema arise: a true local minimum and a saddle point.

The stationary values of $a$ are given by the stationary solutions of Eq. (3), which defines the S–shaped response curve (Fig. 2) of the internal field amplitude to the external field.

In the bistability region $0<\beta<\beta_{\mathrm{crit}}$, there are two stable stationary states 1, 2 and one unstable state S, which lies on a self–intersecting trajectory called separatrix. It divides the phase plane into three regions: the two inner regions 1 and 2 with the corresponding stable states inside them and the outer region 3. The stable state in region $1$ has a lower field amplitude, while the stable state in region $2$ has a higher field amplitude.

II.2 Fokker–Planck equation in the presence of white noise

In any realistic system, noise and damping due to interactions with the environment are always present. They result in appearance of the damping term with dimensionless damping constant $\vartheta\equiv\gamma/\Delta$ and additional random field $\xi$ in the right–hand side of the equations of motion.

[TABLE]

The effect of damping is that the field amplitude relaxes to one of the stable stationary states. Noise has the opposite effect. First, it results in small random deviations from the stationary states. Second, it can induce transitions between the stationary states. At weak noise intensity, these transitions are exponentially rare.

In the case of the white noise (4), it is possible to derive the FPE for the probability density Johne et al. (2009).

[TABLE]

Since the transitions between the stationary states are very rare, the relaxation consists of two stages. At first, the relaxation to the quasi–stationary distribution occurs in each region of the phase space at time scales determined by the inverse damping constant. Then, at a much slower rate, the probability distribution evolves to the true stationary distribution due to noise–induced transitions between the stable states.

At small damping ($\vartheta\ll 1$) and weak noise ($\vartheta/Q\gg 1$), a significant simplification of the 2D FPE is possible. Weak damping and noise give only a small correction to the motion along the phase trajectories. So, it is natural to average the distribution function in each region of the phase space along the trajectory and define the approximate $P_{i}(t,H(A,A^{*}))$, $i=1,2,3$.

Different trajectories with the same quasienergies can exist in regions 1 and 3 [Fig. 1]. By averaging the full FPE, one gets the 1D FPE in quasienergy space Vogel and Risken (1990), Dmitriev and D’yakonov (1986):

[TABLE]

The expressions for $K_{i}(E),D_{i}(E),T_{i}(E)$ were derived in Refs.Vogel and Risken (1990), Maslova (1986), Dmitriev and D’yakonov (1986), Maslova et al. (2007) and are reproduced in the Appendix. $T_{i}(E)$ is the period of motion along the trajectory with quasienergy $E$ in the region $i$, and $K_{i}(E)$ and $D_{i}(E)$ are the drift and the diffusion coefficients in quasienergy space in the region $i$.

This Fokker–Planck equation should be solved in every region of the phase space. The full solution should be obtained by applying the boundary conditions near the separatrix, which include the continuity of the probability distribution and the conservation of the flow:

[TABLE]

The stationary distribution can be obtained by setting the flow $J_{i}(E)$ to zero, if the tunneling effects are neglected (the discussion of the tunneling effects is given below).

II.3 Relative occupation of two stable states

The general formula for the stationary distribution function follows immediately from (6):

[TABLE]

[TABLE]

The distribution has maxima in the vicinity of states $1$, $2$, i.e., at the corresponding quasienergies $E_{1}$, $E_{2}$ Maslova et al. (2007), Vogel and Risken (1990). Outside the neighborhood of $E_{1}$ and $E_{2}$, $P^{\mathrm{st}}_{i}(E)$ is exponentially small. Depending on whether $S_{1}>S_{2}$ or $S_{1}<S_{2}$, the probability density is mostly concentrated around either state $1$ or state $2$.

Numerical evaluation of $S_{1}$ and $S_{2}$ shows that $S_{1}=S_{2}$ at $\sqrt{\beta_{0}/\beta_{\mathrm{crit}}}=0.29$. Therefore $\beta_{0}$ corresponds to the threshold pumping intensity: at $\beta<\beta_{0}$ the oscillator mostly remains in state 1 with a small amplitude, and at $\beta>\beta_{0}$ it mostly remains in state 2 with a large amplitude. Thus the choice of the most probable state is defined by a single parameter $\beta$, and the switching from one most probable state to another occurs at the universal threshold value $\beta=\beta_{0}$. The width of the threshold region is determined by the characteristics of the noise. When $|\beta-\beta_{0}|\sim Q/\vartheta$, both states have comparable probabilities.

II.4 Transition rates between different stable states

The relaxation of a nonlinear driven oscillator happens in two stages. The first stage is the fast relaxation to the quasi–stationary distribution which occurs independently in regions $1$ and $2$. After that, the slow relaxation to the real stationary state occurs, which is governed by rare fluctuation–induced transitions between the stable states.

Every solution of the FPE can be expressed as a sum over eigenfunctions:

[TABLE]

where $\lambda$ and $P^{\lambda}$ are the solutions of the eigenvalue problem

[TABLE]

The eigenvalues of the FPE provide an important information about the kinetics of the system. As shown in Fig. 4, in the bistability region the lowest nonzero eigenvalue is several orders of magnitude smaller than the rest of eigenvalues. It determines the last stage of the relaxation process which was described above.

At small $Q/\vartheta$, the lowest eigenvalue $\lambda$ is exponentially small. Thus it is possible to use the perturbation theory for $P^{\lambda}_{i}$ Vogel and Risken (1990). In each region of the phase space, the distribution function up to the first order in $\lambda$ is given by

[TABLE]

where $P^{\mathrm{st}}(E)$ is the stationary distribution (8). Using the continuity of the probability distribution and the conservation of the flow, one gets the following expression for the lowest eigenvalue $\lambda$:

[TABLE]

It is clear that the analytical expression (13) fits well the numerical results everywhere in the bistability region except in the vicinity of its edges [Fig. 4].

The lowest eigenvalue non–monotonically depends on the value of $\beta$ and achieves its minimum at $\beta=\beta_{0}$. At $\beta<\beta_{0}$ ($\beta>\beta_{0}$), it corresponds to the escape rate $\lambda_{21}$ ($\lambda_{12}$) from the higher (lower) amplitude state to the lower (higher) one, which drops (rises) with the growing external field intensity. At the threshold intensity $\beta_{0}$, $\lambda_{12}$ and $\lambda_{21}$ have the same values.

II.5 A tunneling term in the Fokker–Planck equation

The trajectories in regions $1$ and $3$ can have the same quasienergy. Thus there is a possibility of quantum tunneling between them. In the quasiclassical language, it can be described as the tunneling term in the FPE:

[TABLE]

Here, $\lambda_{T}$ is the tunneling rate. It can be calculated in the quasiclassical limit with the tunneling amplitude proportional to $e^{-\frac{\Delta}{\hbar\alpha}S_{\mathrm{tunn}}(E)}$ using the Fermi’s golden rule:

[TABLE]

At small $\lambda_{T}$, the stationary solution $P^{T}_{i}(E)$ can be obtained by the perturbation approach similar to (12):

[TABLE]

From this equation, we can determine the most probable quasienergy states in regions 1 and 3 from which the tunneling occurs. It is defined by the minimum of $S_{\mathrm{total}}$,

[TABLE]

For calculation of the quasiclassical tunneling exponent $S_{\mathrm{tunn}}$ it is necessary to rewrite the classical Hamiltonian (1) in real variables $p$ and $q$: $A=\frac{q+ip}{\sqrt{2}}$:

[TABLE]

The tunneling action is defined as an integral over the classically inaccessible area $p^{2}<0$:

[TABLE]

The functions $p_{1,3}(q)$ for a specific quasienergy $E$ are determined as the solutions of the equation $\mathcal{H}(p,q)=E$:

[TABLE]

The turning points $q_{1}$, $q_{2},q_{c}$ are defined by the conditions $p_{1}^{2}(q_{1})=0$, $p_{3}^{2}(q_{2})=0$, and $p_{1}^{2}(q_{c})=p_{3}^{2}(q_{c})$.

The resulting quasiclassical tunneling exponent has an integral representation

[TABLE]

At small $\beta/\beta_{\mathrm{crit}}$, it can be approximated as

[TABLE]

Now, from the expression (17) we can estimate the quasienergy state which is optimal for tunneling.

At $\vartheta/Q\gg\Delta/\hbar\alpha$, $S_{\mathrm{total}}$ has a minimum near $E\sim E_{1}$, Therefore tunneling transitions occurs directly between the lower–amplitude stable state and the corresponding state from region 3. On the contrary, at $\vartheta/Q\ll\Delta/\hbar\alpha$ the ”total action” $S_{\mathrm{total}}$ has a minimum near $E\sim E_{\mathrm{sep}}$. So, tunneling occurs between the states with quasienergy close to $E_{\mathrm{sep}}$, and the noise–induced transitions dominate.

We concentrate on the case $\vartheta/Q\gg\Delta/\hbar\alpha$. In this limit, the leading term in the tunneling action at $E=E_{1}$ is

[TABLE]

The pre–exponential factor has the order of $\Delta$. It can be evaluated by matching the quasiclassical solutions near the turning points.

Tunneling between the quasiclassical trajectories effectively occurs when the quasienergies obtained from the Bohr–Sommerfeld quantization rule become almost equal. In this case, tunneling leads to an exponentially small splitting between them. As will be shown below, even at finite $\beta$, this occurs when $2\Delta/\hbar\alpha$ is exactly integer. In this case, the tunneling rate between the classical trajectories in regions 1 and 3 with closest quasienergies is estimated as

[TABLE]

When $2\Delta/\hbar\alpha$ is not integer, and $\gamma\ll\Delta$, one has

[TABLE]

In both cases, the tunneling rate is proportional to $\beta^{\frac{2\Delta}{\hbar\alpha}}$. At integer $2\Delta/\hbar\alpha=m$, the tunneling rate can be treated as the probability of $m$–photon resonant transition between the real energy states of the nonlinear oscillator. So, the tunneling processes in the presence of a resonant external field and the multi–photon transitions between the energy states of a nonlinear oscillator are the similar effects Keldysh (1965).

The same expression for the tunneling amplitude in the lowest non–vanishing order can be also obtained in the framework of the quantum–mechanical perturbation theory for multi–photon transitions Larsen and Bloembergen (1976):

[TABLE]

For $k=0$,

[TABLE]

The state with $k=0$ corresponds to the point 1 on the phase portrait. So, for a driven bistable system, the probability of $m$–photon transition calculated quantum–mechanically (26) is the same as the tunneling probability between the degenerate quasienergy states in the quasiclassical treatment.

The same nature of tunneling effects and multi–photon ionization of atoms in a strong electromagnetic field was first demonstrated by L. V. Keldysh Keldysh (1965).

The presence of tunneling modifies both the distribution function and the relaxation rate. If $\lambda_{T}$ is small, its effect can be taken into account within the perturbation theory. The ratio of probability densities of states $1$, $2$ modifies as follows:

[TABLE]

According to this formula, tunneling leads to a decreasing probability to be in state $1$. Tunneling also changes the total transition rate between the stable states:

[TABLE]

Here, $\lambda_{\mathrm{noise}}$ is defined by (13).

The behavior of the transition rate between the stable states in the presence of tunneling is depicted in Fig. 5. Tunneling transitions shift the threshold value of the external field intensity towards lower values and increase the threshold values of the transition rate.

III Bistability in quantum oscillator

III.1 Quantum quasienergy states and squeezing

The Hamiltonian for a quantum bistable oscillator in the rotating–wave approximation is given by:

[TABLE]

The operators $\hat{a}$ and $\hat{a}^{\dagger}$ are the creation and annihilation operators of the internal oscillator field. In the quasiclassical limit, $\sqrt{\hbar}\hat{a}$ and $\sqrt{\hbar}\hat{a}^{\dagger}$ correspond to the classical field amplitudes. In the following, we set $\hbar=1$.

The exact eigenstates of (30) should be obtained numerically by diagonalization of the Hamiltonian matrix. However, qualitatively, the structure of eigenstates can be understood using the classical analogy.

From the Bohr–Sommerfeld quantization rule, one concludes that the eigenstates of the Hamiltonian correspond to the discrete set of trajectories in the classical phase portrait (Fig. 1). However, the real picture is a bit more complicated because the quantum tunneling should also be taken into account. This is because the classical phase portrait has different regions with the same quasienergy, i.e., regions I and III. So, the real eigenstates may correspond not only to single trajectories but also to superpositions of two trajectories with the same quasienergy.

The possibility of quantum tunneling is closely connected to the degeneracy of eigenstates in the Hamiltonian (30) at $f=0$. At $f=0$, the Hamiltonian commutes with $a^{\dagger}a$, and the states with $k$ excitation quanta are the eigenstates of the Hamiltonian. Their quasienergy is

[TABLE]

For integer $2\Delta/\alpha$, the states with $k$ and $m-k$ excitation quanta become degenerate. At small but nonzero $f$, these states can mix: the true eigenstates are the superpositions of $|k\rangle$ and $|m-k\rangle$. In the quasiclassical language, this corresponds to tunneling between degenerate classical trajectories. Numerical diagonalization shows that such mixing occurs only when $2\Delta/\alpha$ is very close to an integer.

To provide some illustration to this qualitative picture, we calculated the eigenstates of the Hamiltonian (30) in the coherent basis:

[TABLE]

where $|z\rangle$ is a normalized coherent state. The function $\psi_{n}(z)$ corresponding to the $n$–th eigenstate has a maximum near the contour line of the classical Hamiltonian $H(a,a^{*})=E_{n}$. This means that the quantum state $|n\rangle$ corresponds to the classical motion along the trajectory $H(a,a^{*})=E_{n}$.

An important property of the quasienergy states is that the states corresponding to the higher amplitude stable point are squeezed. This can be shown by the mean–field expansion:

[TABLE]

The mean value of $\hat{a}$ is defined from the equation $\frac{\partial H}{\partial a}(\left\langle a\right\rangle_{2},\left\langle a^{*}\right\rangle_{2})=0$, which corresponds to the classical stable states $i$. For small $\beta$, the mean field value in the higher amplitude stable state 2 is $\langle a\rangle_{2}\approx\sqrt{\Delta/\alpha}(1+\sqrt{\beta}/2)$.

The quadratic part of the Hamiltonian takes the form

[TABLE]

We diagonalize this Hamiltonian using the Bogolyubov transformation:

[TABLE]

Let us consider the uncertainties in two quadratures $\hat{q}$ and $\hat{p}$, $\hat{a}=\frac{\hat{q}+i\hat{p}}{\sqrt{2}}$. Squeezing is more pronounced in the higher amplitude stable states $\langle a\rangle_{2}$.

[TABLE]

The quadratic approximation is correct when $\sqrt{\beta}$ is larger than $\alpha/\Delta$. When $\alpha/\Delta\ll 1$, this is fulfilled almost in the entire region of bistability, and the relations (36) are valid.

The minimum possible uncertainty of $\hat{q}$ is at $\beta\sim(\alpha/\Delta)^{2}$, where it can be estimated as

[TABLE]

Thus, the uncertainty in $\hat{q}$ quadrature can be far beyond the quantum limit.

As we have shown in the previous section (28), the tunneling effects increase the occupation of the stable state 2 with a higher amplitude and therefore enhance the generation of squeezed states.

III.2 Quantum kinetic equation

Let us assume that the system is weakly interacting with the environment:

[TABLE]

We assume that the correlation functions of damping operators are delta–correlated:

[TABLE]

where $N$ is the number of noise quanta.

With such assumptions, the density matrix evolution can be described by the master equation Haken (1965), Risken (1965), Graham and Haken (1970). Drummond and Walls (1980), Risken et al. (1987):

[TABLE]

If $\gamma$ is small compared to $\Delta$, the density matrix is almost diagonal in the basis of eigenstates $|n\rangle$, and the master equation reduces to the rate equation for probabilities $P_{n}$ to be in the $n$–th eigenstate:

[TABLE]

This equation is a quantum analog of (6). The evolution of the density matrix has the same features as the evolution of the distribution function for a classical oscillator with bistability. At infinite time, the density matrix evolves to the stationary distribution. The relaxation to $P^{\mathrm{st}}_{n}$ consists of two stages. The first stage corresponds to the relaxation to the quasi–stationary distribution. Its typical time is $\gamma^{-1}$. The second stage is the relaxation to the true stationary state. This stage is very slow and happens due to transitions between the classical stationary states. These transitions can be induced by quantum fluctuations as well as by thermal noise.

Formally, the general solution of (41) reads

[TABLE]

The lowest nonzero eigenvalue is much smaller than all other eigenvalues. Therefore, at large $t$ only the term with the lowest nonzero $\lambda$ should be retained in Eq. (42).

The density matrix relaxes to the true stationary distribution with the rate $\lambda_{\mathrm{min}}$, which can be interpreted as the rate of fluctuation–induced transitions between the stable states.

III.3 The quasiclassical limit

One can show that the continuous limit of (41) is the (6). As it was mentioned above, every eigenstate corresponds to a trajectory on the classical phase portrait, and the hybridization of the trajectories from regions $1$ and $3$ can be neglected unless $2\Delta/\alpha$ is very close to an integer. Thus, in the quasiclassical limit, the distribution function $P_{n}$ weakly depends on $n$ in each of the regions of the phase space. Moreover, the transition rates $w_{n,n^{\prime}}\equiv\tilde{w}_{\bar{n},k}$, $\bar{n}=\frac{n+n^{\prime}}{2}$, $k=n^{\prime}-n$ decrease fast with an increasing value of $|k|$ and weakly depend on $\bar{n}$, which is close to $n$. In this case, it is possible to perform a gradient expansion of $P_{n}$, $w_{nn^{\prime}}$ in (41):

[TABLE]

[TABLE]

In (44), we took into account that $w_{n,n+k}=\tilde{w}_{n+\frac{k}{2},k}=\tilde{w}_{\bar{n},k}$, $\bar{n}=n+\frac{k}{2}$. Keeping the terms up to the second order in $k$, one obtains the differential equation for $P_{n}$:

[TABLE]

where the coefficients $A(n)$ and $B(n)$ are given by the expressions (41) for probabilities $w_{nn^{\prime}}$:

[TABLE]

[TABLE]

Here $\tilde{T}(\epsilon_{n})$ is the period of the classical motion with quasienergy $\epsilon_{n}$.

In the quasiclassical limit, the averages over the quantum quasienergy states transform to time–averages over the classical trajectories. Thus, in the quasiclassical limit, $A(n)$ and $B(n)$ are expressed as line integrals over the classical trajectories:

[TABLE]

After a change of variables $\frac{\Delta}{\alpha}T(E)dE=2\pi dn$, $t\Delta=\tau$, and $\gamma/\Delta=\vartheta$ the equation transforms to the classical FPE (6). The coefficient $A(n)$ transforms to $\vartheta K(E)$, and $B(n)$ transforms to $QT(E)D(E)$, where

[TABLE]

$E=\alpha\epsilon/\Delta^{2}$ is the dimensionless quasienergy, and $T(E)=\Delta\tilde{T}(E\Delta^{2}/\alpha)$ is the dimensionless period as in (6).

III.4 Results and discussion

Qualitatively, the behavior of $P^{\mathrm{st}}_{n}$ in the diagonal approximation resembles the behavior of $P^{\mathrm{st}}_{i}(E)$ of a classical oscillator, as $P^{\mathrm{st}}_{i}(E)$ is the classical limit of $P^{\mathrm{st}}_{n}$ (here $i$ indicates the classical region of the phase space). As $P^{\mathrm{st}}_{i}(E)$, it consists of two sharp peaks which can be attributed to the classical stable states $1$ and $2$. Below (above) the threshold value of the external field, the state $1$ ($2$) dominates.

We directly compared the distributions obtained from the classical FPE and from the quantum master equation. In the classical limit $P^{\mathrm{st}}_{n}$ equals $(2\pi\alpha/\Delta)P_{i}^{\mathrm{st}}(E_{n})$, $E_{n}=\alpha\epsilon_{n}/\Delta^{2}$, where $P^{\mathrm{st}}_{i}(E)$ is the classical distribution function for a dimensionless Hamiltonian (2) with number of noise quanta defined by (49). The index $i$ corresponding to the classical region of the phase space is uniquely defined for each eigenstate unless $2\Delta/\alpha$ is an integer. In the latter case, the classical FPE should be derived from the quantum master equation more carefully. It can be obtained only after choosing the proper basic quasienergy states. One should deal with the quantum states corresponding to the trajectories in the regions of phase space 1 and 3, but not with their superposition.

For the classical case, it was shown that the change of the most probable stable state takes place at $\sqrt{\beta}\equiv\sqrt{\alpha f^{2}/\Delta^{3}}\approx 0.29$. From Fig. 7, it is clear that for a rather high number of noise quanta $N\gg\alpha/\Delta$ and for a large non–integer $2\Delta/\alpha\gg 1$ the probability distribution for a quantum oscillator obtained from the master equation (41) coincides with the classical distribution over quasienergies. For a small number of noise quanta, the situation is more complicated even for large $\Delta/\alpha$. Even though at large $\Delta/\alpha$, the quasiclassical approximation for matrix elements of $\hat{a}$ is valid, the quantum distribution function doesn’t coincide with the classical one due to quantum fluctuations and tunneling effects. For a quantum oscillator described by the rate equation (41), $\beta_{0}$ is no more the universal threshold parameter. Not only the parameter $\beta$ matters, but also $\Delta/\alpha$ and $N$ matter as well. The threshold value $\beta_{q}$ in the quantum limit can exceed the classical value for a non–integer $2\Delta/\alpha$, and it is below the classical value at integer $2\Delta/\alpha$. At a small number of noise quanta, when $N\lesssim\alpha/\Delta$, the quantum effects including the tunneling processes govern the transitions between different stable states. In this case, the distribution function and the relaxation rate strongly depend on the value of $\Delta/\alpha$.

When $\Delta/\alpha$ becomes close to an integer, the quantum distribution even at $N>\alpha/\Delta$, $\Delta/\alpha\gg 1$ doesn’t coincide with the classical one due to tunneling effects. This is demonstrated in Figs. 7(b) and 8(b). Tunneling between the classical regions $1$ and $3$ leads to a decreasing threshold intensity $\beta_{q}$ and an increasing relaxation rate, as was mentioned in Sec. II.5.

Moreover, as shown in Fig. 8, the value of $N$ also influences the threshold external field intensity and the threshold relaxation rate. As $N$ is decreased, the threshold external field intensity rises and reaches $\sqrt{\beta/\beta_{\mathrm{crit}}}\approx 0.5$ at $N=0$. Such behavior does not appear in the quasiclassical FPE solutions.

To clarify the influence of the fluctuation–induced transitions on statistical properties of the internal oscillator field, one should calculate the second–order correlation function $g^{2}(0)$:

[TABLE]

When one of the stable states dominates, $g^{(2)}(0)\sim 1$. In a narrow region of fluctuation–induced transitions $g^{(2)}(0)$ is significantly larger. For several values of $N$, we have calculated $g^{2}(0)$ as a function of two parameters: $\sqrt{\beta}$ and $\Delta/\alpha$.

For each $\Delta/\alpha$, the correlation function has similar dependence on $\sqrt{\beta}$: there is a sharp peak at some value of $\sqrt{\beta/\beta_{\mathrm{crit}}}$ between [math] and $1$, which indicates the change of the most probable state. However, the peak is not placed at $\sqrt{\beta/\beta_{\mathrm{crit}}}\approx 0.29$ as the classical FPE predicts. Its position is an oscillating function of $\Delta/\alpha$ with sharp minima when $2\Delta/\alpha$ is an integer, which corresponds to enhanced tunneling between degenerate quasienergy states. The amplitude of oscillations decreases as $2\Delta/\alpha$ increases, when one approaches the quasiclassical limit. This fact points to the quantum nature of these oscillations.

IV Conclusions

We have derived the quasiclassical kinetic equations for the probability distribution over quasienergy states of a nonlinear driven oscillator, taking into account the tunneling effects. The stationary distribution for a wide range of system parameters and the typical relaxation rate have been determined. It was shown that the relaxation consists of two stages. At first, the relaxation to the quasi–stationary distribution occurs in each region of the phase space at time scales determined by the inverse damping constant. Then, at exponentially large times, the probability distribution evolves to the true stationary one. The relaxation to the true stationary state happens due to fluctuation–induced transitions between the quasiclassical stable states. In the classical limit, if tunneling is neglected, there exists a universal threshold value of the field intensity responsible for switching between the most probable stable states of the system. Taking into account the tunneling effects renders this value non–universal. The tunneling transitions lead to a decrease of the threshold value of external field intensity and to an up to one order of magnitude increase of the fluctuation–induced transition rate between the stable states in the threshold area.

For a driven quantum nonlinear oscillator, we demonstrated that the quasienergy state corresponding to the classical higher amplitude stable state is squeezed. The degree of squeezing is determined by the ratio of nonlinearity and detuning, and the uncertainty of one of the oscillator quadratures can be much lower than the usual quantum limit.

As tunneling transitions increase the occupation of the higher amplitude stable state, the generation of squeezed states can be enhanced in the presence of tunneling effects.

Also we demonstrated that the quasienergy states become superpositions of trajectories from different regions of the phase space. This happens whenever the detuning is an integer or half–integer multiple of the nonlinear shift per quantum. This happens due to a multi–photon resonance between the real eigenstates of the nonlinear oscillator. It was shown that such resonance can be described in terms of tunneling between the quasienergy states in different regions of the classical phase space.

The kinetics of the quantum oscillator was investigated using the quantum master equation. It was shown that in the limit of large detuning–nonlinearity ratio, large number of thermal photons and in the absence of multi–quantum resonance, the classical FPE in the quasienergy space is a continuous limit of the quantum master equation. Importantly, a large value of the detuning–nonlinearity ratio is not sufficient for the validity of the classical FPE, because at a weak noise the quantum effects become especially pronounced. The relaxation rate and the threshold intensity of the external field are both very sensitive to the detuning–nonlinearity ratio. At an integer or half–integer detuning–nonlinearity ratio, the relaxation rate can increase up to several orders of magnitude and the threshold value of the external field intensity shifts towards lower values. In this case, tunneling between degenerate quasienergy states and the multi–photon resonant transitions between the original states of the nonlinear oscillator can be treated as the same effect.

Finally, it was shown that the second–order correlation function of the internal field strongly rises near the threshold pumping intensity, which indicates super–Poissonian statistics of the internal oscillator field.

Acknowledgements.

This work was supported by RFBR grants 19–02–000–87a and 18–29–20032mk.

Appendix A The coefficients of the classical Fokker–Planck equation

The coefficients of the classical FPE are defined as follows:

[TABLE]

For them, we obtained the following integral representations:

[TABLE]

[TABLE]

[TABLE]

The limits of the integration are the roots or the equation

[TABLE]

For energies $E_{2}<E<E_{\mathrm{sep}}$ corresponding to the classical region $2$, this equation has only two real roots. For energies $E_{\mathrm{sep}}<E<E_{1}$ corresponding to regions 1 and 3, there are four real roots. $t_{1}<t_{2}<t_{3}<t_{4}$. To obtain $T_{1}$, $K_{1}$ and $D_{1}$, the limits of integration should be $t_{1}$, $t_{2}$, and for $T_{3}$, $K_{3}$, $D_{3}$, they should be $t_{3}$, $t_{4}$. Finally, for $E>E_{1}$, which corresponds only to region 3, there are two real roots once again.