Hybrid phase-space--Fock-space approach to evolution of a driven nonlinear resonator
HHybrid phase-space–Fock-space approach to evolution of a driven nonlinear resonator
Mostafa Khezri
1, 2, ∗ and Alexander N. Korotkov Department of Electrical and Computer Engineering,University of California, Riverside, California 92521, USA Department of Physics, University of California, Riverside, California 92521, USA (Dated: November 9, 2018)We analyze the quantum evolution of a weakly nonlinear resonator due to a classical near-resonantdrive and damping. The resonator nonlinearity leads to squeezing and heating of the resonator state.Using a hybrid phase-space–Fock-space representation for the resonator state within the Gaussianapproximation, we derive evolution equations for the four parameters characterizing the Gaussianstate. Numerical solution of these four ordinary differential equations is much simpler and fasterthan simulation of the full density matrix evolution, while providing good accuracy for the systemanalysis during transients and in the steady state. We show that steady-state squeezing of theresonator state is limited by 3 dB; however, this limit can be exceeded during transients.
I. INTRODUCTION
Nonlinear quantum oscillators have been a subject ofvarious studies for a long time [1–5]. The renewed inter-est in this system is caused by the wide use of microwaveresonators in superconducting quantum computing cir-cuits [6, 7], as well as reaching a quantum regime fornanomechanical resonators [8–11]. In particular, dur-ing dispersive measurement of superconducting qubits[6, 12–15], nonlinearity of the measurement resonator isinduced by its coupling with the qubit; this nonlinear-ity causes significant deviations from the standard dis-persive regime in the case of a moderately or stronglydriven resonator [16–18]. The nonlinearity of Josephson-junction-based resonators is used in experiments for near-quantum-limited microwave signal amplification [19–22].Driven nonlinear resonators can produce squeezedstates [3, 4, 23–26] (note that quantum squeezing isclosely related to classical fluctuations, e.g., [27–29]).Even though squeezed states are usually discussed forparametrically driven linear resonators [30, 31] (in opticsa nonlinear material can be used to produce a paramet-ric drive at a doubled frequency), there is a similaritybetween these two systems [19, 32, 33]. In particular, itcan be shown that a nonlinear resonator near the bifur-cation point at large photon numbers is equivalent to adegenerate parametric amplifier driven with a detunedpump [33]. Squeezed states can be used to improve mea-surement accuracy [34, 35] in a range of applications,such as gravitational wave detectors [36], superconduct-ing qubit readout [37–42], and nano/micromechanical po-sition measurement [11, 43, 44]. There is currently asignificant experimental interest in producing squeezedmicrowave states with Josephson parametric amplifiers[21, 45–49]; the self-developing squeezing due to the non-linearity of a microwave resonator (with revival and for-mation of “cat” states) has also been demonstrated ex-perimentally [50]. ∗ email: [email protected] It is well known that the steady state of a parametri-cally driven resonator cannot be squeezed beyond 3 dB[30, 51, 52]; in other words, any (instantaneous) quadra-ture variance is not less than 1 / P -representation, linearization of evolution equation, for-malism of quasienergies, etc. Usually the transientsare neglected and only the steady state is analyzed.Moreover, most of the research has been focused onthe regimes close to bifurcation or within the bistabilityrange, in particular, with the goals to analyze switch-ing between the quasistable states and to analyze am-plification properties near the bifurcation point. In thispaper we are mainly interested in the opposite regime:far from the bifurcation and/or bistability, so that theeffects of nonlinearity are not yet very strong. This a r X i v : . [ qu a n t - ph ] O c t regime is relevant to the measurement of superconductingqubits, in which the weak nonlinearity of the microwaveresonator is induced by its interaction with the qubit.Nevertheless, this weak nonlinearity may lead to a sig-nificant self-developing squeezing of the microwave field[37], which affects qubit measurement fidelity. Anotherdifference of our analysis from most of the previous stud-ies is that we are mainly interested in transients, notthe steady state. This is also motivated by the impor-tance of transients in fast measurement of superconduct-ing qubits. Even though our motivation mainly comesfrom the use of weakly nonlinear microwave resonatorsfor qubit measurement, our results are equally applica-ble to the quantum dynamics of driven nanomechanicalresonators, which always show some nonlinearity [65].In this paper, we analyze the evolution of a coherentlydriven weakly nonlinear resonator using a hybrid phase-space–Fock-space approach [18]. This approach is basedon the observation that quantum state evolution due tononlinearity can be easily described in Fock space, whilethe effect of the drive and dissipation for a linear res-onator is well described in phase space. We show thatfor large photon numbers, a Gaussian state [66] in phasespace has also an approximately Gaussian form in Fockspace, thus obtaining a rather simple conversion betweenthe Fock-space and phase-space representations withinthe Gaussian-state approximation. The conversion equa-tions are then used to derive reasonably simple first-orderordinary differential equations, describing state evolutiondue to drive, dissipation, and weak nonlinearity.These evolution equations are for one complex andthree real parameters, which characterize the Gaussianstate of the resonator. The complex parameter describesthe center of the Gaussian state in the phase plane; itsevolution is given by an essentially classical equation,which takes into account nonlinearity. The three realparameters are Fock-space parameters, which after con-version into the phase space correspond to the minimumand maximum quadrature variances (therefore to squeez-ing and “unsqueezing”) and to the phase of the minimum-variance quadrature. The product of the minimum andmaximum variances (ratio of unsqueezing and squeezing)corresponds to an effective temperature, which can besignificantly higher [26] than the bath temperature. Wenote that our approach is physically similar to lineariza-tion of fluctuations around the classical trajectory withinthe Gaussian approximation [29], even though it is basedon a different framework.After deriving the hybrid phase-Fock-space evolutionequations, we numerically compare their results with themaster (Lindblad) equation simulations. We find quitegood accuracy, with an inaccuracy scaling inversely pro-portional to the number of photons in the system. Eventhough our approximation formally requires large num-ber of photons, it still works well when the resonator evo-lution starts from the ground state. In our simulationswith a few hundred photons in the system, the typical in-fidelity compared with the master equation simulations is about 10 − − − , while being faster by a factor of over10 (fractions of a second instead of hours). Comparedwith the coherent-state approximation, our method forthe simulated cases is more accurate by about a factorof 10 , which indicates the importance of taking into ac-count self-developing squeezing and heating.Thus, our main result in this paper is the derivation ofrelatively simple and computationally efficient equations,which describe the quantum evolution of a driven anddamped weakly nonlinear resonator in the case of largephoton numbers. As an example of using these equations,we derive the 3 dB squeezing limit discussed above forthe steady state and numerically show that this limitcan be exceeded during transients. Note that we analyzeonly the state of the resonator (intracavity field), whilethe analysis of the reflected field is left for future studies(the problem is that for a non-pulsed propagating field,the standard definition of squeezing is applicable only inthe steady state, so for transients we will need to modifythe definition; analysis will probably require the use ofeither the input-output theory [66–69] for the linearizedsystem or the approach of weak measurements [70, 71]).The range of validity for our approach seems to beessentially the same as for validity of the Gaussian ap-proximation. Note that for small number of photons inthe resonator, the resonator is practically linear, whilefor large number of photons, the resonator is practicallysemi-classical, and in both cases the Gaussian approxi-mation is applicable. This is why our approach workswell in a rather wide range, except the vicinity of thebifurcation point, where unsqueezing becomes too large;also, within the bistability region our approach cannotdescribe gradual mixing of quasistable states, which cor-responds to classical switching between them. We ana-lyze the accuracy of our approach numerically, by com-paring its results with results of simulations based on themaster equation.The paper is organized as follows. In Sec. II we de-scribe the system and pose the problem. In Sec. IIIwe review the Gaussian states and corresponding phase-space evolution equations for a driven and damped lin-ear resonator. Then in Sec. IV A we introduce Fock-spaceGaussian states and discuss their equivalence to the usual(phase-space) Gaussian states in the case of large pho-ton numbers, with explicit conversion relations betweenparameters of the phase-space and Fock-space represen-tations. Using these conversion relations, in Sec. IV Bwe combine the Fock-space evolution due to nonlinearitywith the phase-space evolution due to drive and damp-ing, thus deriving the hybrid phase-Fock-space evolutionequations, which are the main result of this paper. Sec-tion V is devoted to analysis of the numerical accuracyof our approach. We start with calculating the fidelityof the conversion between the Gaussian and Fock-spaceGaussian states in Sec. V A, and then in Sec. V B wecompare results of the hybrid evolution equations withthe master equation simulations. In Sec. VI the hybridevolution equations are used to show that steady-statesqueezing of the resonator state is limited by 3 dB, andit is also shown numerically that squeezing during tran-sients can exceed the 3 dB limit. We conclude in Sec.VII. In Appendix A we discuss derivation of the Gaussianstate evolution equations for a linear resonator under co-herent drive and damping. In Appendix B we show thatat large photon numbers, a Fock-space Gaussian statecan be approximated by a phase-space Gaussian state,and derive the corresponding conversion relations. Ap-pendix C discusses analytical results for squeezing in thesteady state. II. SYSTEM AND PROBLEM
We analyze the quantum state evolution of a weaklynonlinear resonator, which is coherently (classically)driven at frequency ω d and damped due to energy relax-ation with rate κ at bath temperature T b . The goal is tofind a reasonably simple approximate description of thisevolution, suitable for large number of photons in the res-onator (we will use the terminology of photons, thoughfor a mechanical resonator the terminology of phononswould be more appropriate).Without damping, the laboratory-frame Hamiltonianof the considered system is ( (cid:126) = 1) H lf = H lfr + H lfd , (1) H lfr = (cid:88) n E ( n ) | n (cid:105)(cid:104) n | , E ( n ) = n − (cid:88) k =0 ω r ( k ) , (2) H lfd = 2Re[ ε ( t ) e − iω d t ] ( a † + a ) , (3)where | n (cid:105) is n th eigenstate of the resonator, with cor-responding eigenenergy E ( n ) expressed via the resonatorfrequency ω r ( n ) = E ( n +1) − E ( n ), which slightly changeswith the level number [we use E (0) = 0], ε ( t ) is thecomplex amplitude of the drive at frequency ω d , and a = ˆ x + i ˆ p is the annihilation operator, while a † = ˆ x − i ˆ p isthe creation operator. Here ˆ x and ˆ p are normalized po-sition and momentum operators, ˆ x = ˆ X (cid:112) mω r0 / p = ˆ P / √ mω r0 , where ˆ X and ˆ P are actual positionand momentum operators, m is effective mass, and inthe normalization we use ω r0 ≡ ω r (0); however, thisparticular value is not important, since we assume aweak nonlinearity, | ω r ( n ) − ω r (0) | (cid:28) ω r (0). The as-sumption of weak nonlinearity also allows us to use thestandard matrix elements for the annihilation operators, (cid:104) k | a | n (cid:105) = √ n δ n − ,k . Note that for a linear resonator, ω r ( n ) = ω r0 , the Hamiltonian (2) reduces to the stan-dard form H lfr = ω r0 a † a . Within the rotating wave ap-proximation (RWA), the drive Hamiltonian (3) becomes H lfd = ε ( t ) e − iω d t a † + ε ∗ ( t ) e iω d t a . The RWA is natu-ral for a weakly nonlinear resonator and near-resonantdrive, | ω d − ω r ( n ) | (cid:28) ω d . In some cases RWA missesexperimentally important effects [72]; however, it shouldbe sufficient for the simple system we consider here.In the rotating frame based on the drive frequency ω d , the RWA Hamiltonian becomes H rf = H rfr + H rfd with H rfr = (cid:88) n E rf ( n ) | n (cid:105)(cid:104) n | , E rf ( n ) = n − (cid:88) k =0 [ ω r ( k ) − ω d ] , (4) H rfd = ε ( t ) a † + ε ∗ ( t ) a. (5)In this paper we will mostly use the rotating frame.The evolution of the system density matrix ρ due toHamiltonian H (in either laboratory or rotating frame)and energy relaxation with rate κ is described by thestandard master equation in the Lindblad form [30, 73,74], ˙ ρ = i [ ρ, H ] + κ ( n b + 1)( aρa † − a † aρ/ − ρa † a/ κ n b ( a † ρa − aa † ρ/ − ρaa † / , (6)where n b = 1 e ω r0 /T b − ω r0 / T b ) −
12 (7)is the average number of thermal photons for the bathtemperature T b . Note that the evolution equation (6)is generally not correct for a nonlinear resonator (e.g.,Appendix B4 of [75]); however, we use it, assuming aweak nonlinearity. The problem with applicability of theLindblad equation (6) stems from the fact that it requiresindistinguishability of the emitted and/or absorbed pho-tons [75]. However, for a weakly nonlinear resonator, thephotons emitted from (absorbed by) different levels haveslightly different frequencies and can be distinguishedspectroscopically if the frequency difference exceeds thelevel width. To estimate the effect, let us assume thatunsqueezing is not too large, so the typical number ofphotons is ¯ n ± √ ¯ n , where ¯ n is the average photon num-ber. Then the frequency difference is about √ ¯ n ( dω r /dn ),while the level width is approximately κ ¯ n . Therefore, in-distinguishability requires ¯ n (cid:29) κ − ( dω r /dn ) . For ourtypical parameters used in Sec. V, the nonlinearity isquite small, so that κ − ( dω r /dn ) ∼ − ; therefore,the indistinguishability condition is well satisfied and theLindblad equation (6) is accurate.Solving Eq. (6) numerically in the Fock space, we canfind the resonator state evolution. However, for over ∼
100 average photons in the resonator the numerical so-lution becomes slow, and for over ∼
500 photons it be-comes computationally intractable on a personal com-puter because of too large Hilbert space. Note that over500 photons in the resonator can be used for a dispersivemeasurement of a superconducting qubit [72, 76].In this paper, we develop an approach which permitsa simple analysis of evolution at this large number ofphotons. To a significant extent, the approach is basedon the observation that evolution of a linear resonatorcan be described by Gaussian states in many situations[77]. Using the fact that a weak nonlinearity keeps theevolving state Gaussian (in the leading order), we willfind the corresponding evolution equations. This greatlysimplifies analysis, since a Gaussian state is characterizedby only 5 real parameters, instead of N parameters fora density matrix involving up to N Fock states.We will first review Gaussian states and evolution of adriven linear resonator, and then will show how a Gaus-sian state can be approximately converted into a Fock-space state, for which it is easy to introduce evolutiondue to nonlinearity.
III. EVOLUTION OF A LINEAR RESONATOR
Without nonlinearity, a Gaussian initial state remainsGaussian during evolution, while initially non-Gaussianstate gradually becomes Gaussian [77, 78]. In this sectionwe briefly review properties of the Gaussian states anddiscuss evolution of a linear resonator state due to applieddrive and damping.
A. Brief review of Gaussian states
Gaussian states [66, 79–81] are defined as states forwhich the Wigner function [66, 82] has a Gaussian form(generally with an arbitrary number of dimensions). Fora one-dimensional (single-mode) system with positionoperator ˆ X and conjugate momentum operator ˆ P , theWigner function of a Gaussian state is W ( X, P ) = exp (cid:16) − (cid:126)V T DDD − (cid:126)V (cid:17) π (cid:112) Det(
DDD ) (8)where (cid:126)V = ( X − X c , P − P c ) T , X c = (cid:104) ˆ X (cid:105) , P c = (cid:104) ˆ P (cid:105) , andelements of the covariance matrix DDD are D = D X = (cid:104) ˆ X (cid:105) − (cid:104) ˆ X (cid:105) , D = D P = (cid:104) ˆ P (cid:105) − (cid:104) ˆ P (cid:105) , and D = D = D XP = (cid:104) ˆ X ˆ P + ˆ P ˆ X (cid:105) / − (cid:104) ˆ X (cid:105)(cid:104) ˆ P (cid:105) . The Husimi Q -function, Glauber-Sudarshan P -function and densitymatrix (in X or P space) of a Gaussian state have aGaussian form as well [66, 83].For a linear resonator with Hamiltonian H lfr = ω r a † a (constant frequency ω r ), we can introduce the dimen-sionless (normalized) operators of position and momen-tum in the standard way as ˆ x = ˆ X/ (2 σ x, gr ) and ˆ p =ˆ P / (2 σ p, gr ), where σ x, gr and σ p, gr are the standard de-viations of the ground state in the position and mo-mentum representations, so that ˆ x = ( a + a † ) / p = ( a − a † ) / i . For the normalized operators, theWigner function W ( x, p ) has exactly the same form asEq. (8), except now (cid:126)V = ( x − x c , p − p c ) T , x c = (cid:104) ˆ x (cid:105) , p c = (cid:104) ˆ p (cid:105) , and elements of the covariance matrix are now D = D x = (cid:104) ˆ x (cid:105) − (cid:104) ˆ x (cid:105) , D = D p = (cid:104) ˆ p (cid:105) − (cid:104) ˆ p (cid:105) , and D = D = D xp = (cid:104) ˆ x ˆ p + ˆ p ˆ x (cid:105) / − (cid:104) ˆ x (cid:105)(cid:104) ˆ p (cid:105) . Explicit formof the Wigner function for a Gaussian state is W ( x, p ) = (cid:16) π (cid:113) D x D p − D xp (cid:17) − × exp (cid:20) − D p (∆ x ) + D x (∆ p ) − D xp ∆ x ∆ p )2( D x D p − D xp ) (cid:21) , (9) FIG. 1. Phase-space illustration of a Gaussian state. Theellipse corresponds to one standard deviation for the quadra-ture operators along any direction. It is also the contour linefor the Wigner function being a factor √ e less than its maxi-mum value. The ellipse center has coordinates ( x c , p c ), whichon the complex plane correspond to (cid:104) a (cid:105) = x c + ip c . Theminimum and maximum quadrature variances are D − b and D + b , respectively. The minimum-variance-direction angleis Θ /
2. In the rotating frame we use notation θ instead of Θ. where ∆ x = x − x c and ∆ p = p − p c . The Wigner func-tions (8) and (9) are normalized as (cid:82) W ( X, P ) dX dP = (cid:82) W ( x, p ) dx dp = 1.With the quadrature operator along direction ϕ de-fined as ˆ x ϕ ≡ ae − iϕ + a † e iϕ x cos ϕ + ˆ p sin ϕ, (10)the variance σ x ϕ ≡ (cid:104) ˆ x ϕ (cid:105) − (cid:104) ˆ x ϕ (cid:105) of this quadrature forthe Gaussian state is σ x ϕ = D x cos ϕ + D p sin ϕ + 2 D xp cos ϕ sin ϕ. (11)Let us introduce real variables D > b ≥ D ≡ D x + D p , b ≡ ( D x − D p ) D xp , (12)then the quadrature variance (11) can be rewritten as σ x ϕ = D − b cos(2 ϕ − Θ) , (13)Θ = arctan (cid:18) D xp D x − D p (cid:19) + π D x − D p )] . (14)Equation (13) shows that D − b and D + b are the min-imum and maximum quadrature variances respectively,and the direction of the minimum quadrature makes theangle Θ / x -axis (see Fig. 1). Note that( D + b )( D − b ) = D x D p − D xp . (15)The Wigner function in the rotated “diagonal basis”with x d being the coordinate along the minimum quadra-ture is W ( x d , p d ) = (cid:16) π (cid:112) ( D − b )( D + b ) (cid:17) − × exp (cid:20) − ( x d − x dc ) D − b ) − ( p d − p dc ) D + b ) (cid:21) , (16)where x d + ip d = ( x + ip ) e − i Θ / and similarly x dc + ip dc =( x c + ip c ) e − i Θ / . This formula shows that the contourlines for the Wigner function in the phase space of x and p are ellipses (Fig. 1).The Husimi Q -function [30] for the Gaussian statecan be obtained using the standard relation Q ( x, p ) = π (cid:82) W ( x (cid:48) , p (cid:48) ) e − x − x (cid:48) ) +( p − p (cid:48) ) ] dx (cid:48) dp (cid:48) . In particular, inthe diagonal basis we find Q ( x d , p d ) = (cid:16) π (cid:112) ( D − b + 1 / D + b + 1 / (cid:17) − × exp (cid:20) − ( x d − x dc ) D − b + 1 / − ( p d − p dc ) D + b + 1 / (cid:21) . (17)We see that the Q -function (17) has the same Gaussianform as the Wigner function (16), but variances for theboth axes are increased by 1 / a , a , and a † a , D = 12 (cid:104) (cid:104) a † a (cid:105) + 12 − (Re (cid:104) a (cid:105) ) − (Im (cid:104) a (cid:105) ) (cid:105) , (18) b = 12 (cid:104) (cid:2) Re (cid:104) a (cid:105) − (Re (cid:104) a (cid:105) ) + (Im (cid:104) a (cid:105) ) (cid:3) + (cid:0) Im (cid:104) a (cid:105) − (cid:104) a (cid:105) Im (cid:104) a (cid:105) (cid:1) (cid:105) / (19)Θ = arctan (cid:18) Im (cid:104) a (cid:105) − (cid:104) a (cid:105) Im (cid:104) a (cid:105) Re (cid:104) a (cid:105) − (Re (cid:104) a (cid:105) ) + (Im (cid:104) a (cid:105) ) (cid:19) + π { (cid:104) a (cid:105) − (Re (cid:104) a (cid:105) ) + (Im (cid:104) a (cid:105) ) ] } , (20) x c + ip c = (cid:104) a (cid:105) . (21)Besides introducing the Gaussian states via the Wignerfunction, it is also possible to introduce them as displacedsqueezed thermal states (DSTS) [83–85], so that the den-sity matrix is ρ DSTS = D ( α ) S ( ξ ) ν n th S ( ξ ) † D ( α ) † , (22)where α = (cid:104) a (cid:105) = x c + ip c is the phase-plane state center, D ( α ) = exp( αa † − α ∗ a ) is the displacement operator, S ( ξ ) = exp[ ξ ∗ a − ξ ( a † ) ] is the squeezing operatorwith squeezing parameter ξ = re i Θ (the angle Θ / ν n th is the thermal state,defined as ν n th = 11 + n th ∞ (cid:88) k =0 (cid:18) n th n th (cid:19) k | k (cid:105)(cid:104) k | , (23)where | k (cid:105) is k th Fock state and n th = Tr( a † a ν n th ) isthe average number of thermal photons. Note that Eq.(23) describes an equilibrium state of a linear resonatorat finite temperature without drive, and in that case n th is equal to the thermal photon number for the bath, n b ,given by Eq. (7). However, in the non-equilibrium caseconsidered in this paper, n th is not equal to n b . It is still possible to define an effective temperature T eff for aGaussian state (22) via the same relation,coth( ω r / T eff ) = 1 + 2 n th . (24)Note that the average photon number ¯ n for a Gaussianstate has a contribution proportional (but not equal) to n th ,¯ n = Tr( a † a ρ DSTS ) = | α | +(1+2 n th ) sinh r + n th , (25)while from Eq. (18) we find a simple expression¯ n = | α | + 2 D − / . (26)To relate parameters r and n th of the DSTS state tothe parameters of the Gaussian state (9), we can calculateaverages (cid:104) a (cid:105) , (cid:104) a (cid:105) , and (cid:104) a † a (cid:105) for the state (22), and usethese results to find the variances D x = (1 / n th / r − sinh 2 r cos Θ) , (27) D p = (1 / n th / r + sinh 2 r cos Θ) , (28) D xp = − (1 / n th /
2) sinh 2 r sin Θ . (29)Comparing Eqs. (27)–(29) with Eqs. (12)–(14), we findthe equivalence for n th = 2 (cid:112) ( D + b )( D − b ) − , tanh 2 r = bD , (30)and the same angle Θ.As follows from the discussion above, a Gaussian stateis determined by five real parameters. Two parameters, x c and p c , define the state center on the phase plane; it isconvenient to use their complex combination α = x c + ip c .Three real parameters define the “shape” (see Fig. 1),which can be characterized either by D x , D p , and D xp or by D , b , and Θ or by r , Θ, and n th . A Gaussianstate is in general a mixed state. A pure Gaussian stateis a minimum-uncertainty squeezed state, characterizedby 4 real parameters; for such a state D x D p − D xp =( D − b )( D + b ) = 1 /
16 and n th = 0. A coherent state ischaracterized by only 2 real parameters, which define thecenter; then D x = D p = D = 1 / D xp = b = n th = 0,and Θ is not important.Note that our discussion in this section used the labo-ratory frame. In this frame, the evolution due to Hamil-tonian H lfr = ω r a † a (in the absence of drive and damping)rotates the state center in Fig. 1 clockwise with angularvelocity ω r . Moreover, the whole phase-space picture inFig. 1 rotates clockwise with ω r . This means that pa-rameters D and b do not change with time, while theangle Θ / d (Θ / /dt = − ω r , and therefore˙Θ = − ω r . Since D and b do not change, the param-eters r and n th are also constant – see Eq. (30). In therotating frame based on the frequency ω d , the picture inFig. 1 additionally rotates counterclockwise with angularvelocity ω d , so that the net evolution is clockwise rota-tion with angular velocity ω r − ω d . Thus, in the rotatingframe, the parameters D , b , r , and n th are the same asin the laboratory frame, while the rotating-frame angleparameter θ is related to Θ as θ = Θ + 2 ω d t, (31)and it evolves as ˙ θ = − ω r − ω d ). Descriptions of theGaussian states in the rotating and laboratory frames arepractically the same, except Θ is replaced with θ and ω r isreplaced with ω r − ω d , as expected for the rotating-frameHamiltonian H rfr = ( ω r − ω d ) a † a . Note, however, thatthe conversion between the actual position and momen-tum operators ( ˆ X , ˆ P ) and the corresponding normalizedoperators (ˆ x , ˆ p ) should still be based on the actual fre-quency ω r and not on ω r − ω d . The relation between thelaboratory frame and the rotating frame is discussed inmore detail in the Appendix A. Evolution in the presenceof drive and damping is discussed next. B. Evolution equations
For a linear harmonic oscillator with H lfr = ω r a † a , theevolution (6) due to drive (3) and damping κ at bathtemperature T b , preserves state Gaussianity and leads tothe following evolution equations in the laboratory frame [43, 77, 86, 87],˙ x c = ω r p c , (32)˙ p c = − ω r x c − κp c − εe − iω d t ) , (33)˙ D x = 2 ω r D xp , (34)˙ D p = − ω r D xp − κD p + ( κ/
2) coth( ω r / T b ) , (35)˙ D xp = − ω r ( D x − D p ) − κD xp . (36)Note that the evolution of the state center ( x c and p c )is decoupled from the evolution of the variances, and thedrive ε contributes only to ˙ p c (as a classical force). Thestate center oscillates with the resonator frequency ω r (in-trinsically, neglecting effects of κ and ε ), while the vari-ances oscillate with doubled frequency, 2 ω r . Also notethat Eqs. (32)–(36) do not rely on the RWA assumption.Using the RWA (which symmetrizes coordinates x and p ) and going into the rotating frame based on the drivefrequency ω d , so that the Gaussian state center is char-acterized by a slowly changing complex number β in thestandard phase space, β = ( x c + ip c ) e iω d t , (37)from Eqs. (32)–(36) we can derive (see Appendix A) thefollowing evolution equations [73, 88] (see also [83, 89])for the parameters β , D , b , and θ ,˙ β = − i ( ω r − ω d ) β − ( κ/ β − iε, (38)˙ D = − κD + ( κ/
4) coth( ω r / T b ) , (39)˙ b = − κb, (40)˙ θ = − ω r − ω d ) . (41) Note that the drive does not affect evolution of thediagonal-basis variances D ± b ; however, the short-axisdirection θ/ ω r − ω d , similar to the rotation of the state center.Equations (38)–(41) are the starting point of our anal-ysis. They describe evolution of a linear resonator usingthe phase-space language. However, to include nonlin-earity, we will need to approximately convert them intothe Fock-space representation. From now on, we will useonly the rotating frame. IV. EVOLUTION OF A WEAKLY NONLINEARRESONATORA. Fock-space Gaussian state
Generalizing the idea of Ref. [18], let us introduce astate, for which the density matrix in the basis of eigen-states | n (cid:105) (Fock space) has the following form, ρ mn = 1 (cid:112) πW | β | exp (cid:20) − ( n + m − | β | ) W | β | − ( n − m ) W | β | (cid:21) × exp (cid:20) iφ β ( n − m ) − i K | β | (cid:16) n + m − | β | (cid:17) ( n − m ) (cid:21) . (42)We call it a Fock-space Gaussian state (because ofquadratic dependence on n and m inside exponents) or,following the terminology of Ref. [18], a sheared Gaus-sian state (because of a shearing effect produced by the K -term in the phase space). The state (42) is character-ized by five real parameters: | β | , φ β , W , W , and K .Note that a physical ρ mn requires0 < W ≤ W . (43)As shown in the Appendix B, in the case | β | (cid:29) W , W , and K are on the order of unity) this stateis approximately equal to the standard Gaussian statediscussed in Sec. III, so that β = e iφ β | β | (44)is (approximately) the state center, while the (approxi-mate) conversion relations for the parameters D , b , and θ are D = 18 (cid:20) W + W (1 + 16 K ) (cid:21) , (45) b = (cid:113) D − W / (16 W ) , (46) θ = 2 φ β + arctan (cid:16) KW D − W / (cid:17) + ( π/
2) [1 − sign( D − W / . (47)The conversion becomes exact for | β | → ∞ .While in the leading order (cid:104) a (cid:105) = e iφ β | β | for the Fock-space Gaussian state (42), more accurate calculationsshow the next-order correction proportional to | β | − , (cid:104) a (cid:105) = e iφ β (cid:20) | β | − W + 1 /W − | β | − K W | β | − i KW | β | (cid:21) . (48)The overlap fidelity between the Gaussian and Fock-space Gaussian states becomes somewhat better if thiscorrection is taken into account, so that a slightly shiftedcenter corresponds to the same (cid:104) a (cid:105) for the Gaussian andFock-space Gaussian states (see numerical results in Sec.V A). However, for simplicity we will not use the centercorrection (48) unless specifically mentioned.Note that the trace of the state (42) is not exactly 1;however, the difference is negligible (exponentially small)for | β | (cid:29)
1. The Fock-space Gaussian state (42) is in gen-eral mixed; it becomes pure if W = W , and in this caseit reduces to the sheared Gaussian state introduced inRef. [18]. [Note a misprint in Eq. (33) of Ref. [18], wherethe last exponent should actually be − iK ( n −| β | ) / | β | .]Comparing Eqs. (45) and (46) with Eq. (30), we find auseful relation for the thermal photon number, n th = ( (cid:112) W /W − / , (49)which is equivalent to the relation W /W = coth ( ω r / T eff ) = 16( D + b )( D − b ) . (50)Note that the ratio of the variances, ( D + b ) / ( D − b ),and the angle θ/ − φ β are both functions of only twoparameters: K and W W .The quadrature variance σ x ϕ along a direction ϕ forthe state (42) can be calculated as σ x ϕ = D − b cos(2 ϕ − θ ) from Eqs. (45)–(47). In particular, for the directionalong β ( ϕ = φ β ) we find the variance σ x ϕ = W /
4, whilefor the orthogonal direction ( ϕ = φ β + π/
2) we find thevariance σ x ϕ = 1 / (4 W ) + 4 K W .As follows from Eq. (47), in the case K = 0, the shortaxis (minimum variance) is either along the direction of β ( θ/ φ β ) or orthogonal to it ( θ/ φ β + π/ β is W / ϕ = φ β + π/
2] thevariance is 1 / W , the short axis is along β if W W < β if W W > | β | →∞ ) from the Fock-space parameters W , W , and K tothe phase-space parameters D , b , and θ , the inverse con-version is given by equations W = 4[ D − b cos( θ − φ β )] , (51) W = D − b cos( θ − φ β )4( D − b ) , (52) K = b sin( θ − φ β )4[ D − b cos( θ − φ β )] . (53)The main idea of introducing the Fock-space Gaussianstate (42) is that it has a simple evolution due to res-onator nonlinearity. Let us consider the evolution only due to Hamiltonian (4), i.e., with ε = κ = 0. Then ρ nm ( t ) = ρ nm (0) exp {− i [ E rf ( n ) − E rf ( m )] t } . Compar-ing this phase evolution with the second line of Eq.(42) and expanding the resonator frequency ω r ( n ) inEq. (4) up to first order around n ≈ | β | (assumingthat nonlinearity is practically constant within the range | n − | β | | (cid:46) √ W | β | ), we find evolution equations˙ φ β = − [ ω r ( | β | ) − ω d ] , (54)˙ K = 12 | β | dω r ( n ) dn (cid:12)(cid:12)(cid:12)(cid:12) | β | , (55)where we neglected discreteness of ω r ( n ). We see that β rotates due to detuning of the resonator frequency ω r ( | β | ) at the state center from the rotating-frame fre-quency ω d (as should be expected), while nonlinearitychanges K , leading to accumulation of the quadraticphase factor in Eq. (42).We emphasize that a weak nonlinearity approximatelypreserves the Fock-space Gaussian form (42), and there-fore approximately preserves the Gaussian-state form inthe phase space, assuming a large photon number | β | .Since the evolution due to the drive and damping alsopreserves the Gaussian-state form, as discussed in Sec.III (for weak nonlinearity we can use approximately thesame matrix elements of operator a in the Fock space asfor a linear oscillator), the state remains approximatelyGaussian in both phase and Fock spaces during the com-bined evolution. B. Hybrid phase-Fock-space evolution equations
We have separately described the evolution due to non-linearity, Eqs. (54)–(55), and due to drive and damping,Eqs. (38)–(41). The combined evolution is simply thesum of the corresponding terms. However, Eqs. (38)–(41)assume the phase-space representation of Fig. 1, whileEq. (55) is based on the Fock-state representation (42).Thus, we need to convert the equations into a commonrepresentation using the conversion formulas (44)–(47).We will characterize the evolving state by four param-eters: β ( t ), W ( t ), W ( t ), and K ( t ). We call it a hybridrepresentation, since β is a phase-space parameter, while W , W , and K originate from the Fock-space descrip-tion.As discussed in Sec. IV A, evolution due to nonlinearityproduces Eq. (55) for ˙ K , the center β evolves as˙ β = − i [ ω r ( | β | ) − ω d ] β, (56)while W and W do not evolve, ˙ W = ˙ W = 0. Notethat Eq. (56) essentially implies that the average num-ber of photons in the resonator is ¯ n ≈ | β | , neglectingcorrections in Eq. (25).To find evolution of parameters W , W , and K due todrive and damping, we write Eqs. (39)–(41) expressingthe time derivatives ˙ D , ˙ b , and ˙ θ via the partial deriva-tives over the parameters of the conversion equations(45)–(47), ∂D ∂W ˙ W + ∂D ∂W ˙ W + ∂D ∂K ˙ K = − κD + ( κ/
4) coth( ω r0 / T b ) , (57) ∂b∂W ˙ W + ∂b∂W ˙ W + ∂b∂K ˙ K = − κb, (58) ∂θ∂W ˙ W + ∂θ∂W ˙ W + ∂θ∂K ˙ K + 2 d [arg( β )] dt = 0 , (59)where in the last term of Eq. (59) we need to use˙ β = − βκ/ − iε , not including the evolution (56) due todetuning. This is because the evolution (56) compensatesthe right-hand-side term of Eq. (41), which we thereforedo not write in Eq. (59). Another justification of writingEq. (59) in this way is that we consider evolution onlydue to drive and damping [not due to detuning, whichis already considered in Eq. (56)]; then the angle θ doesnot change in time, and we should exclude the detuningterm from ˙ β .Equations (57)–(59) with the partial derivatives ob-tained from Eqs. (45)–(47), give us a system of three lin-ear equations for ˙ W , ˙ W , and ˙ K . Solving this system,we find˙ W = 8 KW Re( ε/β ) + κ [coth( ω r0 / T b ) − W ] , (60)˙ W = 8 KW Re( ε/β )+ κW [1 − W (1 + 16 K ) coth( ω r0 / T b )] , (61)˙ K = 14 [( W W ) − − (1 + 16 K )] Re( ε/β ) − κ ( K/W ) coth( ω r0 / T b ) . (62)Note that in the term coth( ω r0 / T b ) we neglect changingresonator frequency because of the weak nonlinearity as-sumption. In the special case when κ = 0, Eqs. (60)–(62)reduce to Eq. (47) of Ref. [18].Finally, combining the terms from Eqs. (55)–(56) (forevolution due to nonlinearity) and from Eqs. (60)–(62)(for evolution due to drive and damping), we obtain thehybrid phase-Fock-space evolution equations˙ β = − i [ ω r (¯ n ) − ω d ] β − κ β − iε, ¯ n ≈ | β | , (63)˙ W = 8 KW Re( ε/β ) + κ [coth( ω r0 / T b ) − W ] , (64)˙ W = 8 KW Re( ε/β )+ κW [1 − W (1 + 16 K ) coth( ω r0 / T b )] , (65)˙ K = (cid:18) W W − K (cid:19) Re( ε/β ) − κKW coth( ω r0 / T b ) + 12 | β | dω r ( n ) dn (cid:12)(cid:12)(cid:12)(cid:12) n = | β | . (66)Evolution equations (63)–(66) complemented with theconversion formulas (45)–(47) are the main result of this paper. To our knowledge, this approach to the quantumevolution of a weakly nonlinear resonator has never beenused previously.Equations (63)–(66) describe evolution of five real pa-rameters of a Gaussian state. Equation (63) describingevolution of the state center (2 real parameters) is de-coupled from the other three equations. The equationsare approximate and assume | β | (cid:29) | β | (cid:29)
1, Eqs. (63)–(66) canbe used to numerically analyze evolution starting evenfrom β = 0 with a good accuracy (the numerical resultsare discussed later). There is no divergence of Re( ε/β )in Eqs. (64)–(66) at β = 0 because if β ( t ) = 0, thenclose to this time moment β = − iε ( t − t ) and thereforeRe( ε/β ) = Re[ i/ ( t − t )] = 0. A numerical divergence canbe easily avoided by shifting the denominator of Re( ε/β )by a negligible amount.Equation (63) has a simple physical meaning; it takesinto account that the resonator frequency ω r ( n ) changeswith the photon number n , and approximates n with theaverage photon number ¯ n ≈ | β | . One may think that asimple generalization of Eq. (63) is to use a more accuratevalue for ¯ n from Eq. (25) in ω r (¯ n ) [it would also requireconversion equations (30) and (45)–(47)]. However, nu-merical simulations show that this correction does notalways give a better agreement with full master equationsimulations using Eq. (6). Because of that, we do not usethis correction in the numerical analysis in Secs. V andVI.Note that Eqs. (63)–(66) permit three natural rescal-ings. First, by rescaling the time axis, it is possi-ble to use κ = 1. Second, since discreteness of n is not important in our approach, we can rescale the n axis and normalize nonlinearity, for example setting dω r ( n ) /dn | n =0 = ±
1. Third, non-zero bath temperature T b is equivalent to rescaling W → W coth( ω r0 / T b )and W → W / coth( ω r0 / T b ), while using T b = 0 inEqs. (63)–(66); this leads to D → D coth( ω r0 / T b ) and b → b coth( ω r0 / T b ), with unchanged β and θ .Equations (64)–(66) describe evolution of the Fock-space parameters W , W , and K . It is also possibleto write evolution equations for the phase-space param-eters D , b , and θ . Note that without the last term inEq. (66), Eqs. (64)–(66) exactly correspond to Eqs. (39)–(41). Therefore, we only need to convert the last term in(66) into the phase space, that can be done by using par-tial derivatives from the conversion relations (51)–(53).In this way we obtain the following evolution equations,˙ D = − κD + ( κ/
4) coth( ω r0 / T b ) + 2 η β | β | b sin(∆ θ ) , (67)˙ b = − κb + 2 η β | β | D sin(∆ θ ) , (68) d (∆ θ ) dt = 2 Re( ε/β ) − η β | β | b − D cos(∆ θ ) b , (69)where ∆ θ ≡ θ − β ), η β ≡ dω r ( n ) /dn (cid:12)(cid:12) n = | β | , andevolution of β is still given by Eq. (63). Note that diver-gence in Eq. (69) at β = 0 can be avoided numericallyin the same way as discussed above: by a negligible shiftof β . The divergence in Eq. (69) at b = 0 can also beavoided numerically by a negligible increase of b (physi-cally, this divergence is because ∆ θ is undefined at b = 0).Equations (67)–(69) are equivalent to Eqs. (64)–(66). Wehave checked this equivalence numerically. However, inthe simulations discussed below we used Eqs. (64)–(66)rather than Eqs. (67)–(69). One of the reasons for ourpreference is that evolution of W , W , and K is alwayssmooth, while ∆ θ evolves very fast when b approacheszero, thus potentially creating a problem with numericalsolution of differential equations (even though our simu-lations never suffered from this potential problem).Note that from Eqs. (67) and (68) we can obtain ddt ( D ± b ) = − [ κ ∓ η β | β | sin(∆ θ )] ( D ± b )+ ( κ/
4) coth( ω r0 / T b ) , (70)which shows that for the maximum-variance andminimum-variance quadratures, the effective dampingrate is different, κ eff = κ ∓ η β | β | sin(∆ θ ), and changeswith time. Similarly, the effective bath temperature isalso different, coth( ω r0 / T b ) → ( κ/κ eff ) coth( ω r0 / T b ).Discussion in terms of different effective damping ratesfor the two quadratures makes an obvious connection tothe case of a parametric drive with doubled frequency.We have checked that Eqs. (67)–(69) are consistentwith the results of Ref. [29] for Gaussian variances ofclassical fluctuations around the trajectory (63), causedby classical (complex) white noise √ κ ζ ( t ) applied to theresonator, with the correlation function (cid:104) ζ ∗ ( t ) ζ ( t (cid:48) ) (cid:105) =(1 /
2) coth( ω r / T b ) δ ( t − t (cid:48) ), (cid:104) ζ ( t ) ζ ( t (cid:48) ) (cid:105) = 0 (as in, e.g.,[71]). Note, however, that in order to get correct equa-tions, we had to exchange B with B † in Eq. (3.2.4) of Ref.[29]. The correspondence between Eqs. (67)–(69) and re-sults of Ref. [29] confirms that the quantum squeezing issimilar to squeezing of classical fluctuations, and it alsoshows that our approach is physically similar to lineariza-tion of fluctuations around the classical trajectory withinthe Gaussian approximation.In Appendix C we derive analytical results for D , b ,and ∆ θ in the steady state and discuss their equivalenceto the results of Refs. [2] and [26] for a Duffing oscillator(Kerr nonlinearity). V. NUMERICAL ACCURACY
In this section we discuss numerical accuracy of ourapproach. We start with analyzing fidelity of the con-version between the Gaussian and Fock-space Gaussianstates, and then discuss numerical accuracy of the hy-brid phase-Fock-space evolution equations by comparingresults with full simulation.
FIG. 2. Infidelity 1 − F between the Gaussian and Fock-spaceGaussian states as a function of (real) β for several values ofthe parameters 4( D + b ) and θ/ n th = 0 (solid lines) and n th = 1 / | β | all lines show the scaling | β | − , illustratedby the long-dashed line. A. Fidelity of the conversion
As was discussed in Sec. IV A, the Gaussian state (9)is approximately equal to the Fock-space Gaussian state(42) with the conversion relations (44)–(47), in the caseof large photon numbers, | β | (cid:29)
1. Let us check theaccuracy of this conversion numerically. For that we cal-culate the overlap fidelity F between the states (9) and(42) using the standard definition [90] F = (cid:16) Tr (cid:112) √ ρ ρ √ ρ (cid:17) Tr( ρ ) Tr( ρ ) , (71)where ρ and ρ are the density matrices of the comparedstates. Note that for normalized states the denominatorin Eq. (71) is not needed, but we use the more generalversion (71) because the Fock-space Gaussian state (42)is not exactly normalized. When at least one of the statesis pure, Eq. (71) reduces to the usual state overlap, e.g., F = (cid:104) ψ | ρ | ψ (cid:105) = Tr( ρ ρ ) if ρ = | ψ (cid:105)(cid:104) ψ | and bothstates are normalized.To find the conversion fidelity for a Gaussian state withparameters β , D , b , and θ , we use conversion relations(45)–(47) to find corresponding parameters W , W , and K [ β is the same unless we use the correction (48)], whichgives us the Fock-space Gaussian state (42). Then wecalculate exact Fock-space representation of the Gaussianstate of (9) using Eq. (22) with parameters | ξ | and n th obtained from the relations (30) (using α = β and Θ = θ ). Finally, we use Eq. (71) in the Fock space to find thefidelity F between the Gaussian and Fock-space Gaussianstates. Note that F does not depend on the phase arg( β )for a fixed value of θ/ − arg( β ), so it is sufficient toconsider arg( β ) = 0, i.e., β = | β | ; this is what we assumebelow in the numerical analysis of the conversion fidelity;in this case θ/ − arg( β ) → θ/ . . . . . . θ/ (2 π ) − − ( − F ) | β | FIG. 3. Scaled infidelity (1 − F ) | β | as a function of the short-axis angle θ/
2. Solid lines are for n th = 0 (pure states) and4( D + b ) = 8, 4, 2, and 1 (top to bottom); dashed lines are for n th = 1 / D + b ) = 8, 4, and 2 (top to bottom). Weused β = 40, which is sufficiently large so that the presentedresults do not depend on | β | . Figure 2 shows infidelity 1 − F as a function of | β | on a log-log scale for several values of other parameters:4( D + b ) = 1, 2, and 4 (this parameter is the long-axisvariance compared with the coherent state; we call it“unsqueezing factor”), θ/ π/
2, and π/ n th = 0 and 1 / D + b ) and θ/
2; solid and dashed lines correspondto n th = 0 and 1 / θ when 4( D + b ) = 1 + 2 n th [see Eq.(30)], then we show only the line θ = 0; also note thatfor n th = 1 / D + b ) ≥ − F ∝ | β | − at large | β | (this scaling is il-lustrated by the long-dashed line). The deviation fromthis dependence at small | β | is mainly caused by two rea-sons. First, the “shoulder” feature may develop when | β | < (cid:112) D − b cos θ ) ≤ (cid:112) D + b ) because then | β | < √ W in Eq. (42) and thus the Gaussian approx-imation near n = 0 becomes inaccurate (less than 3standard deviations). Second, deviation from the scal-ing | β | − starts to develop when 1 − F (cid:38) .
05 because F cannot exceed 1; actually, a natural metric for distancebetween the states is arccos( √ F ) [90], which is approx-imately √ − F when 1 − F (cid:28)
1; for this metric theabove condition is √ − F (cid:38) .
22. From Fig. 2 we con-clude that the scaling 1 − F ∝ | β | − is almost perfect if | β | > (cid:112) D + b ) and 1 − F < . − F ) | β | forsufficiently large | β | (here we used β = 40), as a func-tion of the short-axis angle θ/
2. We used parameters4( D + b ) = 1, 2, 4, and 8, while n th = 0 (solid lines) and1 / θ/ D + b ) = 1+2 n th , since in this case the long-axis and short-axis variances coincide, D + b = D − b .When 4( D + b ) > n th , the local minima of the infi- D + b ) − ( − F ) | β | n th = 0 n th = 1 / FIG. 4. Solid lines: scaled infidelity (1 − F ) | β | maximizedover the angle θ/ β = 40), as a function of the quadru-pled long-axis variance 4( D + b ). The upper (blue) solid lineis for n th = 0, the lower (orange) solid line is for n th = 1 / delity are reached at θ/ θ/ π/
2; both thesecases correspond to K = 0 in Eq. (42) [note that K = 0minimizes the state center shift in Eq. (48), which affectsinfidelity, as discussed below]. For relatively small valuesof 4( D + b ), the minimum is reached at θ/ D + b ), the mini-mum infidelity is at θ/ π/ θ/ ± π/
4. Note that the infidelity dependence on θ/ π , and the dependence is symmetric about thepoints θ/ θ/ π/ − F ) | β | maximized over the angle θ/ D + b ) (long-axis variance in units of the coherentstate variance) for the case n th = 0 (zero effective tem-perature). We see that this line can be approximatelyfitted by the formula1 − F ≈ .
04 [4( D + b )] | β | , (72)which is drawn as the dotted black line.The infidelity scaling 1 − F ∝ ( D + b ) can be crudelyunderstood as a consequence of the Fock-space Gaussianstate center shift described by Eq. (48). Considering forsimplicity the case K = 0 and W = W (cid:28) n th = 0, θ/ | β | ≈ − (8 W | β | ) − ≈− ( D + b ) / (2 | β | ) along the short axis. The relative shiftcompared with the “width” of the state along the shortaxis is then ∆ | β | / √ D − b ≈ − [4( D + b )] / / (4 | β | ).Since the infidelity scales quadratically with this relativeshift, 1 − F ∝ (∆ | β | / √ D − b ) , we obtain the scaling1 − F ∝ [4( D + b )] / | β | .1The same numerical scaling of the infidelity in Eq. (72)indicates that the state center shift may play a signifi-cant role in fidelity reduction. To check this hypothesis,we used the correction from Eq. (48) to produce Gaus-sian and Fock-space Gaussian states with the same (cid:104) a (cid:105) by making a small compensating shift of β . The corre-sponding result for the infidelity 1 − F is shown by theupper (blue) dashed line in Fig. 4. As we see, the cor-rection has really decreased the infidelity; however, theimprovement is only by a factor of about 2, so the scalingis approximately the same as in Eq. (72), with the factor0.04 replaced by 0.02. We have also checked that numer-ical optimization of the infidelity over the center shift ofthe Fock-space Gaussian state [instead of using Eq. (48)]produces practically the same result. The infidelity de-crease by a factor of about 2 can be crudely understoodin the following way. The Fock-space Gaussian state hasa slightly crescent (non-elliptical) shape of the Wignerfunction in the phase plane. Slightly shifting its center,it is possible to improve the state fidelity compared withthe Gaussian state (which has a perfect elliptical shape);however, this improvement cannot be very significant.Now let us discuss the lower (orange) lines in Fig.4, for which n th = 1 / T eff = 0 . ω r ); as above, the dashed line takes into ac-count the center correction (48), while the solid line iswithout the correction. We see that non-zero n th im-proves the fidelity compared with the case n th = 0 for thesame long-axis variance D + b (the short-axis variance inthis case is increased by a factor of 4). The improvementcan be qualitatively understood using the above deriva-tion based on the state center shift: since the short-axis“width” is now larger, the relative inaccuracy is smaller,thus decreasing the infidelity. Note, however, that suchderivation would predict infidelity reduction by a factorof 4, while numerically the distance between the upperand lower solid lines in Fig. 4 is less than a factor of 2.5.Comparing the solid and dashed orange lines, we see thatthe state center correction decreases the infidelity; how-ever, the improvement is only by crudely a factor of 1.5,even less than in the zero-temperature case.We can make the following conclusions from the nu-merical results discussed in this section. First, the infi-delity of the conversion between the Gaussian and Fock-space Gaussian states is not larger than in Eq. (72), so theconversion becomes almost perfect for sufficiently large | β | . Second, correction (48) to the state center improvesthe fidelity; however, the improvement is not very sig-nificant (we will not use this correction in analyzing theevolution). Let us also note that the change of effectivetemperature from zero to 0 . ω r ( n th = 1 /
2) did not pro-duce a very significant change in the infidelity.
B. Accuracy of the hybrid phase-Fock-spaceevolution equations
The main result of this paper is the hybrid phase-Fock-space evolution equations (63)–(66), which permit a veryefficient approximate simulation of the state dynamicsfor a slightly nonlinear resonator in the large-photon-number regime. In contrast, full simulation using themaster equation (6) is highly resource-consuming in thisregime because of large Hilbert space. In this section wenumerically analyze the accuracy of our hybrid equationsby comparing the results with the full master equationsimulation.For the numerical analysis let us consider a constantdrive, ε ( t ) = ε , and a constant (Kerr) nonlinearity, ω r ( n ) = ω r0 + nη, (73)which corresponds to the rotating-frame resonator energylevels E rf ( n ) = ( ω r0 − ω d ) n + n ( n − η/
2. We also assumethat initial state is vacuum, ρ (0) = | (cid:105)(cid:104) | . Note that thehybrid evolution equations still work well when initialstate is vacuum, because for sufficiently weak nonlinear-ity, the photon number becomes large before the effectsdue to nonlinearity (e.g., squeezing) become important.Also note that in RWA the considered resonator Hamil-tonian is equivalent to H lfr = P / (2 m ) + ( m/
2) ˜ ω X +( η/ m ˜ ω X , where ˜ ω r0 = ω r0 − η . The difference be-tween the first-excitation frequency ω r0 and the “plasmafrequency” ˜ ω r0 for a Duffing oscillator is negligible be-cause we focus on the regime of large n .In the considered case, the RWA dynamics describedby the master equation (6) depends on five parame-ters: nonlinearity η , drive amplitude ε , initial detuning ω r0 − ω d , damping rate κ , and bath temperature T b char-acterized by the bath photon number n b via Eq. (7).Rescaling the time axis (using κ − as the time unit), itis easy to see that the dynamics depends on four dimen-sionless parameters: η/κ , ε/κ , ( ω r0 − ω d ) /κ , and n b .For simulation using the hybrid evolution equations(63)–(66), it is possible to further reduce the number offree parameters from four to only two (this is not possiblefor full master equation simulation). Since discretenessof n is not used in Eqs. (63)–(66), it is possible to rescale n -axis using κ/ | η | as the unit of n ; this eliminates non-linearity as a free parameter, d ( ω r /κ ) /d [ n/ ( κ/ | η | )] = ± η ). This rescaling renor-malizes the drive amplitude as ( ε/κ ) / (cid:112) κ/ | η | , while notaffecting dimensionless detuning. Furthermore, it is pos-sible to rescale W and W using 1 / coth( ω r0 / T b ) andcoth( ω r0 / T b ) respectively; this eliminates bath temper-ature as a free parameter, such that it can always beassumed zero. Then the rescaled dynamics is deter-mined by only two free parameters: ε (cid:112) | η | /κ / and( ω r0 − ω d ) /κ , and we can use Eqs. (63)–(66) with thefollowing parameters: κ → dω r /dn → ± η ), ε → ε (cid:112) | η | /κ / , ω r0 − ω d → ( ω r0 − ω d ) /κ , and T b →
0; this automatically rescales β as2 tκ − − − − − − F Naive (coherent state)Hybrid phase-FockGaussian-state fit
FIG. 5. Blue (lower) solid line: time dependence of infidelity1 − F ( t ) between the exact solution ρ m ( t ) obtained from themaster equation (6) and state ρ h ( t ) obtained from the hybridevolution equations (63)–(66). Parameters are close to typicalcircuit QED parameters (see text), time t is normalized bythe resonator decay time κ − . Dashed green line: infidelitybetween ρ m ( t ) and its Gaussian-state fit. Red (upper) solidline: infidelity of the conventional approach based on coherentstates. β → β/ (cid:112) κ/ | η | , time as t → κt , variables W and W as W → W coth( ω r0 / T b ) and W → W / coth( ω r0 / T b ),while K does not change.To check accuracy of the hybrid phase-Fock-space evo-lution equations, let us calculate the time-dependent fi-delity F ( t ) [Eq. (71)] between the exact solution ρ m ( t )of the master equation (6) and the state ρ h ( t ) obtainedfrom our approximate hybrid equations (63)–(66). Notethat in the hybrid method we evolve variables β , W , W ,and K , but the resulting state is always converted into aGaussian state using Eqs. (45)–(47), so the fidelity F ( t )is calculated between this Gaussian state and Fock-spacesolution of the master equation [for that the Gaussianstate is represented in the Fock space using Eq. (22)]. Insimulations we will use parameters somewhat close totypical parameters in circuit QED experiments for mea-surement of superconducting transmon qubits; a weaknonlinearity of the resonator in this case is induced bythe qubit nonlinearity; the resonator nonlinearity is muchmore significant when the transmon is in the ground state[18].The lower (blue) solid line in Fig. 5 shows the time-dependent infidelity 1 − F of the calculation based on thehybrid phase-Fock-space evolution equations (63)–(66).Here we used parameters κ/ π = 5 MHz, ω r0 − ω d = 0, η/ π = − .
02 MHz, ε/ π = 32 MHz (this corresponds to100 photons in the steady state), and n b = 3 . × − (this corresponds to T b = 50 mK for ω r0 / π = 6 GHz;we start with the vacuum state instead of the thermalstate, but the difference is negligible). We see a verygood accuracy provided by our approach, with infidelitybelow 10 − . For comparison, the upper (red) solid line shows the infidelity for the conventional naive approach,in which we assume a coherent state of the resonator,with the same center β ( t ) given by Eq. (63). We see thatthe conventional approach fails to describe the evolutionwith a good accuracy, thus emphasizing importance ofconsidering Gaussian states in our approach.For the dashed green line in Fig. 5, at each time t wefitted ρ m ( t ) by a Gaussian state having the same values of (cid:104) a (cid:105) , (cid:104) a (cid:105) , and (cid:104) a † a (cid:105) , and then calculated fidelity betweenthis Gaussian state and ρ m ( t ). Therefore, the dashed lineessentially shows the non-Gaussianity of the actual state ρ m ( t ) (we have checked that numerical optimization overthe state center β does not provide a noticeable furtherimprovement of the infidelity). Comparing the dashedgreen line with the blue solid line, we see that our hy-brid evolution equations (63)–(66) describe the resonatorstate almost as good as this Gaussian-state fit. We havefound numerically that almost all difference between thesolid blue and dashed green lines in Fig. 5 comes from asmall inaccuracy in calculation of the state center usingEq. (63) [see Fig. 6(b)]. We tried to improve this accuracyby using ¯ n from Eq. (25) for the center evolution (63) andalso by using the center correction (48). While this de-creased infidelity for some parameters, it increased it forsome other parameters, so we decided to use the simplestequation (63) for the state center evolution. As followsfrom Fig. 5, this already gives a very good accuracy.To clarify the origin of the “bump” on the lower lines inFig. 5, in Fig. 6(a) we show the corresponding evolutionof “squeezing parameter” 1 / [4( D − b )] (lower lines) and“unsqueezing parameter” 4( D + b ) (upper lines). We seethat the maximum infidelity in Fig. 5 occurs at approx-imately the same time as the maximum unsqueezing inFig. 6(a), thus hinting that the infidelity during evolu-tion originates from a mechanism similar to the infidelitybetween the Gaussian and Fock-space Gaussian states es-timated by Eq. (72). The quantitative comparison showsthat the maximum of the lower solid line in Fig. 5 isabout a factor of 4 smaller than the estimate given byEq. (72), while the steady-state infidelity is smaller thanthis estimate by a factor of 9.The solid lines in Fig. 6(a) are calculated using thehybrid evolution equations (63)–(66), while dashed linesare obtained from the Gaussian-state fit of the master-equation result ρ m ( t ). We see that the dashed and solidlines are very close to each other, indicating that our hy-brid approach is quite accurate in calculating the quadra-ture variances.Note that for a minimum-uncertainty (pure) state, thelower and upper lines (squeezing and unsqueezing) in Fig.6(a) should coincide; the ratio between these parame-ters is coth ( ω r0 / T eff ) – see Eq. (50). From Fig. 6(a)we see that the resonator state is considerably mixed,with the effective temperature T eff significantly exceed-ing [26] the bath temperature T b ; for example, in thesteady state T eff = 98 mK, in contrast to T b = 50 mK.A large corresponding ratio of thermal photon numbers, n th /n b = 17 .
3, indicates that the effective temperature3 tκ . . . . . . . ( D + b ) , / ( D − b ) (a) 4( D + b )1 / D − b )0 2 4 6 8 10 Re( β ) − − − − − I m ( β ) (b) Master equationEq. (62) FIG. 6. Panel (a): “Squeezing factor” [4( D − b )] − (lowerlines) and “unsqueezing factor” 4( D + b ) (upper lines) asfunctions of time, for parameters of Fig. 5. Solid lines are ob-tained from the hybrid evolution equations (63)–(66), dashedlines are obtained from the Gaussian-state fit to the master-equation result ρ m ( t ). Panel (b): Corresponding evolution ofthe state center β ( t ) on the phase plane, with points spacedin time by 0 . /κ . Solid blue line with dots is calculated usingEq. (63), almost coinciding red dashed line with squares show (cid:104) a (cid:105) for ρ m ( t ). T eff in this case is practically independent of the bathtemperature. Indeed, the same simulations with T b = 0showed a very close effective temperature, T eff = 96 mK.In Fig. 6(b) we show evolution of the state center β ( t )on the phase plane for the same parameters as in Figs. 5and 6(a). The dots (and squares) are separated by timeintervals 0 . /κ (which is 15.9 ns); the solid blue line withdots is for calculation using Eq. (63), while the dashedred line with squares shows (cid:104) a (cid:105) for the master-equationsimulation result ρ m ( t ). We see that Eq. (63) is quiteaccurate for calculating the state center. However, thereis a tiny (almost unnoticeable) difference between posi-tions of the dots and squares in Fig. 6(b); as mentionedabove, this tiny shift is mainly responsible for the dif-ference between the lower solid and dashed lines in Fig. . . . . . . Re( α ) − . − . − . − . − . − . I m ( α ) Master equationHybrid phase-Fock
FIG. 7. Contour plot for the Wigner function W ( α ) of theresonator state. The black solid lines are calculated using thehybrid evolution equations (63)–(66), the red dashed lines arecalculated using the master equation (6). The parameters arethe same as in Figs. 5 and 6, the snapshot is taken at time t = 15 /κ . The contours are drawn at the levels of 1 / π , 2 / π ,... 7 / π . The centers are indicated by black and red dots.
5. As another observation, the maximum photon number | β | is achieved at almost the same time as the maximumof 4( D + b ); however, we think that the infidelity bumpin Fig. 5 is caused by the maximum of 4( D + b ) and notby the almost simultaneous maximum of | β | .The main advantage of our method is a simple calcu-lation of the resonator state deviation from a coherentstate. For illustration, Fig. 7 shows the contour plot ofthe Wigner function W ( α ) of the resonator state at timemoment t = 15 /κ (practically the steady state) for thesame parameters as in Figs. 5 and 6. The solid blacklines are calculated for our approximate hybrid-evolutionstate ρ h , while the dashed red lines correspond to the ex-act state ρ m (at this snapshot 1 − F = 2 . × − ). Wesee that our approach gives a quite good approximationfor the Wigner function; the difference is mainly because W ( α ) contour plot for the actual state ρ m has a slightlycrescent shape, while in our Gaussian-state approxima-tion the contours are strictly elliptical. We used Eq. (16)to calculate W ( α ) for the Gaussian state ρ h , while for ρ m we used the formula [91, 92] W ( α ) = 2 π Tr (cid:104) D ( − α ) ρ D ( α ) e iπa † a (cid:105) , (74)in which the displacement operator D ( α ) was appliednumerically in the Fock space.Now let us check numerically the expectation that ourapproach should become more accurate with more pho-tons in the resonator. The solid lines in Fig. 8 show the4 tκ . . . . . . . . . − F n st = 50 n st = 100 n st = 200 FIG. 8. Solid lines: time dependence of infidelity 1 − F ( t )between the simulations based on the master equation andon our hybrid evolution equations, for the stationary-statephoton numbers n st ≈
50, 100, and 200 from top to bot-tom. The corresponding (color-matched, the same order)dashed lines show infidelity of the Gaussian-state fit to themaster-equation simulations. The dimensionless parameters, ε (cid:112) | η | /κ / = 0 .
40 and ( ω r0 − ω d ) /κ = 0, are the same as inFigs. 5–7, while ε and η change from line to line (see text). time-dependent infidelity 1 − F ( t ) for the calculationsusing Eqs. (63)–(66) (compared with the master equa-tion results) for different number of photons. All solidlines correspond to the same normalized drive amplitudeand detuning as in Figs. 5–7: ε (cid:112) | η | /κ / = 0 .
40 and( ω r0 − ω d ) /κ = 0; however, nonlinearity η varies: from topto bottom η/ π = − . − .
02, and − .
01 MHz; corre-spondingly, the drive amplitude ε also varies (with decayrate κ/ π = 5 MHz kept constant): ε/ π = 32 / √ , √ n st ≈ | β st | approximately equalto 50, 100, and 200 from top to bottom (note that thescaled evolution is the same as in Fig. 6). As expected,the solid lines in Fig. 8 show that the infidelity becomessmaller with more photons in the resonator. The scalingis crudely 1 − F ∝ | β st | − , as expected from Fig. 2 andEq. (72).In addition to better accuracy, for larger | β st | ourapproach becomes much more preferable computation-ally in comparison with the master-equation calculations.As an example, for our codes (which are rather simple,Mathematica-based) the calculation of the hybrid evolu-tion ρ h ( t ) for the solid lines in Fig. 8 took about 0.02 sec-onds, while obtaining the numerical master-equation so-lution ρ m ( t ) took 0.2, 1, and 4 hours on a high-end desk-top computer (longer time for larger | β st | ). The master-equation simulation duration scales crudely quadraticallywith the size of the Fock space, while for our hybrid equa-tions there is no scaling with the system size. For thelower solid line in Fig. 8, our method was faster by afactor exceeding 10 .Dashed lines in Fig. 8 show infidelity of the Gaussian- state fit of ρ m ( t ) for the same parameters. Comparing thesolid and dashed lines, we see that most of the infidelity inour approach comes from non-Gaussianity of the actualstate, thus making unimportant any possible improve-ments in the state center calculation by improving Eq.(63). We also see that the fraction of the infidelity com-ing from non-Gaussianity does not change significantlywith changing number of photons.Note that with zero initial detuning, ω d = ω r0 , as-sumed in Figs. 5–8, we automatically avoid the bista-bility region [93, 94] for the steady state of a classicalresonator with Kerr nonlinearity (73). Our method isgenerally not intended to work inside or close to thisbistability region. In particular, quantum treatment for-mally removes the bistability [2] because of transitionsdue to quantum fluctuations (tunneling or quantum ac-tivation [32]), even though the rate of these transitionscan be exponentially small. In contrast, our approachuses the classical equation (63) for the state center evo-lution, showing full bistability. The critical point [93, 94](start of the bistability) occurs at | ˜ ε | = 3 − / ≈ .
44 and∆˜ ω d = √ / ε ≡ ε (cid:112) | η | κ / , ∆˜ ω d ≡ − sign( η ) ω r0 − ω d κ . (75)For larger | ˜ ε | , the bistability range for ∆˜ ω d becomesnon-zero and grows. For a given ∆˜ ω d above √ /
2, thebistability region for the dimensionless drive amplitudeis | ˜ ε − | ≤ | ˜ ε | ≤ | ˜ ε + | , where | ˜ ε ∓ | = ˜ n ± [˜ n ± − ∆˜ ω d ] + ˜ n ± / n ± = [2∆˜ ω d ± (cid:112) ∆˜ ω − / / n is relatedto the photon number n st = | β st | as ˜ n = n st | η | /κ ). Asmentioned above, we should avoid this bistability regionwhen using our approach (63)–(66). We have checkednumerically that in the vicinity of the critical point aswell as near the bistability region, the unsqueezing pa-rameter 4( D + b ) may become large, indicating that ourapproach could become accurate only at very large num-ber of photons.The numerical results presented in this section showthat our approach based on the hybrid evolution equa-tions (63)–(66) typically provides a good accuracy, whichis orders of magnitude better than using the conventionalapproximation based on the coherent-state assumption.On the other hand, our approach is orders of magnitudefaster than the full simulation based on the master equa-tion. VI. dB SQUEEZING LIMIT AND ITSVIOLATION IN TRANSIENTS Squeezing of a resonator state due to Kerr nonlin-earity (73) has been discussed long ago [3, 23–25] (seealso [18]). A somewhat similar squeezing of the vacuumstate can be produced by a parametric drive at the dou-bled frequency [30, 31], and in this case the steady-statesqueezing of the resonator state is always less than 3 dB,5i.e., [4( D − b )] − ≤ W = ˙ W = 0 into Eqs. (64) and (65), wefind that in the steady state1 + 16 K = 2 W / coth( ω r0 / T b ) − W W . (76)Therefore, from Eq. (45) we obtain D = W / [4 W coth( ω r0 / T b )]. Now using Eq. (46) forthe parameter b , we obtain the scaled minimum quadra-ture variance 4( D − b ) = W / [ W coth( ω r0 / T b )] − (cid:112) [ W /W coth( ω r0 / T b )] − W /W . Representing thisresult as4( D − b ) = coth( ω r0 / T b )1 + (cid:113) − coth ( ω r0 / T b ) W /W , (77)we obtain [4( D − b )] − < ω r0 / T b ) ≥ W /W is positive. Thus, squeezing is less than 3 dB inthe steady state.Note that the 3 dB squeezing limit can be approachedonly when the bath temperature T b is zero [so thatcoth( ω r0 / T b ) = 1] and when W /W → ∞ . Corre-spondingly, effective temperature T eff becomes infinitelylarge because n th → ∞ , as follows from Eqs. (49) and(50). We also see that in this case the maximum quadra-ture variance becomes infinitely large, 4( D + b ) → ∞ ,which indicates instability (similar to the case of reach-ing the 3 dB limit for parametric doubled-frequency drive[30, 51]). Using Eqs. (63)–(66), we have checked numer-ically that 3 dB squeezing can be approached near thecritical point and also near the switching point on the up-per branch in the bistability region. As discussed above,our formalism is not actually intended to work in thisparameter range. The hybrid equations do not have anymathematical problems in this range; however, there canbe a problem with accuracy compared to the exact (mas-ter equation) evolution. In particular, when 4( D + b )becomes large near the critical point, the accuracy of theformalism requires a very large number of photons [see estimate (72)]. In addition, within the bistability regionour formalism neglects switching between the quasistablestates caused by fluctuations, so it can be reasonably ac-curate only when the switching rate is very small (thatalso requires a large number of photons). In spite of theseissues, we can still formally use our equations, keeping inmind the potential problems.Even simpler derivation of the 3 dB limit can be ob-tained using Eq. (70). This derivation follows very closelythe underlying physical idea of the derivation [30, 52] forthe case of a parametric drive. From Eq. (70) we find thatin the steady state the unsqueezing and inverse squeezingfactors are4( D ± b ) = coth( ω r0 / T b )1 ∓ η β | β | sin(∆ θ ) /κ . (78)Since D + b >
0, there is a limitation 2 η β | β | sin(∆ θ ) < κ (which is similar to the constraint of the parametric insta-bility). Therefore, for 4( D − b ) the denominator in Eq.(78) is less than 2 (and obviously positive), thus leadingto the inequality 4( D − b ) > (1 /
2) coth( ω r0 / T b ) ≥ / D − b )] − as afunction of time for the dimensionless drive amplitude˜ ε = ε (cid:112) | η | /κ / = 10, no initial detuning, ω d = ω r0 , andzero temperature of the bath. We see that the 3 dB limit(horizontal line) is exceeded repeatedly, even though inthe stationary state the squeezing is below 3 dB. Thisnumerical result was obtained using Eqs. (63)–(66). Tocheck it, we also performed the simulations using themaster equation (6). The dashed red line in Fig. 9(a)shows the corresponding result for the same parametersand η/κ = − .
03 (as discussed above, master equationrequires more dimensionless parameters than the hybridevolution equations); for example, this case can be re-alized with ω d / π = ω r0 / π = 6 GHz, κ/ π = 5 MHz, η/ π = − .
15 MHz, T b = 0, and ε/ π ≈
290 MHz (theseparameters can in principle be realized with a circuitQED setup by increasing the effective resonator nonlin-earity | η | using an increased qubit-resonator coupling).The maximum average number of photons in this case isapproximately 350 (at κt ≈ .
4) – see Fig. 9(c). Com-paring the solid blue and dashed red lines in Fig. 9(a),we see that the master equation gives a slightly smallersqueezing than the hybrid equations, but it still signifi-cantly exceeds the 3 dB value at the peaks. Note that thehybrid-equation calculation took about 0.02 seconds on adesktop computer, while the master-equation simulationtook over 15 hours (the ratio of over 10 ).A noticeable inaccuracy of the squeezing calculationin Fig. 9(a) using the hybrid equations is related to largevalues of the unsqueezing parameter 4( D + b ) shown inFig. 9(b). At the first peak ( κt ≈ .
4) the infidelityestimate using Eq. (72) for | β | ≈
350 gives 0.05, sowe would expect a noticeable inaccuracy. We checkedthat the inaccuracy decreases with decreasing nonlinear-6 . . . . . . . tκ / ( D − b ) (a) Hybrid phase-FockMaster equationWigner width . . . . . . . tκ ( D + b ) (b) Hybrid phase-FockMaster equationWigner width Re( β ) − − − I m ( β ) (c) MasterequationEq. (62) FIG. 9. Panel (a): Squeezing factor [4( D − b )] − as a functionof time t for ω d / π = ω r0 / π = 6 GHz, κ/ π = 5 MHz, η/ π = − .
15 MHz, T b = 0, and ε/ π = 290 MHz (so that ε (cid:112) | η | /κ / = 10). The solid blue line is calculated using thehybrid evolution equations, the dashed red line is obtainedfrom the master equation simulation, and the dotted blackline is the variance of the master-equation Wigner functionalong the short axis. Panel (b): Unsqueezing factor 4( D + b )for the same parameters (solid blue and dashed red lines).Dotted black line is the Wigner function variance along thelong axis. Panel (c): the corresponding evolution of the statecenter β ( t ) on the phase plane. The dots are separated intime by 0 . /κ , larger dots are separated by 0 . /κ . ity | η | /κ while keeping ε (cid:112) | η | /κ / fixed; this increasesthe number of photons, which scales as κ/ | η | . (Sincefurther increase of the photon number is very difficultfor the master-equation simulations, we actually checkedthat the inaccuracy in Fig. 9(a) increases with decreasingnumber of photons by increasing | η | /κ .) Note that theunsqueezing parameters calculated by the hybrid equa- tions and by the master equation [solid blue and dashedred lines in Fig. 9(b)] practically coincide with each other.Figure 9(c) shows the evolution of the state center β ( t )on the phase plane, with dots separated in time by 0 . /κ (larger dots are separated by 0 . /κ ); the results from Eq.(63) and master equation practically coincide with eachother. Comparing Fig. 9(c) with Figs. 9(a) and 9(b),we see that peaks in squeezing and unsqueezing approxi-mately correspond to maxima of the photon number | β | .The minima of the photon number correspond to smallbumps on the lines in Figs. 9(a) and 9(b).We expect that the difference between the solid blueand dashed red lines for the squeezing factor in Fig. 9(a)can be mostly explained by a non-Gaussian shape of theactual states produced by the master equation. Thisnon-Gaussianity can be seen as a slightly crescent shapeof the Wigner function in the phase plane (see Fig. 7),with slightly curved “arms” along the long axis, insteadof the perfect elliptical shape. However, the bendingof the “arms” produces a smaller effect along the shortaxis. To check this hypothesis, we have calculated theWigner function variance along the short axis by numer-ically fitting the master-equation Wigner function alongthe short axis (passing through the state center) with aone-dimensional Gaussian model. The result is shown bythe dotted black line in Fig. 9(a). It is almost indistin-guishable from the blue solid line, thus confirming thatsqueezing calculated by our hybrid-evolution method isessentially the squeezing of the Wigner function along theshort axis (which is slightly different from the usual “in-tegrated” definition based on the quadrature variance,which is affected by bending of the “arms”). In con-trast, the Wigner function variance along the long axis,shown by black dotted line in Fig. 9(b), noticeably differsfrom the quadrature variance shown by the solid blue (ordashed red) line. This is expected because the Wignerfunction along the long axis is significantly more affectedby bending of the “arms”.Figure 10 shows time-dependence of the squeezing fac-tor [4( D − b )] − for various parameters; these resultsare obtained using the hybrid equations (63)–(66). InFig. 10(a) we assume zero initial detuning and zero bathtemperature, ω d = ω r0 , T b = 0, while varying the dimen-sionless drive amplitude, ˜ ε ≡ ε (cid:112) | η | /κ / = 5, 10, and 15.In Fig. 10(b) we keep the amplitude fixed, ˜ ε = 10, andvary the detuning, ∆˜ ω d ≡ sign( η )( ω d − ω r0 ) /κ = − ε = 10 the bistability region starts at ∆˜ ω d = 8 . D + b ); for example, the maximum squeezing factor of5.6 in Fig. 10(a) for ˜ ε = 15 corresponds to 4( D + b ) = 7 . | β | = 13 . κ/η ). Similarly, the maximumsqueezing factor of 7.6 in Fig. 10(b) for ∆˜ ω d = 3 corre-sponds to 4( D + b ) = 16 . | β | = 14 . κ/η ).This means that to observe these large values of squeez-7 . . . . . . . tκ / ( D − b ) (a) ∆˜ ω d = 0 ˜ ε = 15˜ ε = 10˜ ε = 5 . . . . . . . tκ / ( D − b ) (b) ˜ ε = 10 ∆˜ ω d = +3∆˜ ω d = 0∆˜ ω d = − FIG. 10. Time dependence of the squeezing factor [4( D − b )] − , calculated using the hybrid evolution equations (63)–(66). The lines in panel (a) are for zero initial detuning, ω d = ω r0 , zero bath temperature, T b = 0, and dimensionless driveamplitudes ε (cid:112) | η | /κ / = 15, 10, and 5 (from top to bottom).The lines in panel (b) are for ε (cid:112) | η | /κ / = 10, T b = 0, anddimensionless initial detunings ( ω d − ω r0 ) /κ sign( η ) = 3, 0,and − ing, we would need very many photons in the resonator.From Eq. (72) and numerical results in Sec. V B, we ex-pect that validity of our formalism requires¯ n ≈ | β | (cid:29) [4( D + b )] . (79)Therefore, we estimate that for the upper (green) lines inFigs. 10(a) and 10(b) to be reasonably accurate, we needover 500 and 4,000 photons, respectively. Therefore, wecannot check results of Fig. 10 against the master equa-tion. However, since the results of the hybrid equationsand the master equation agree well with each other inthe range where the master equation requires reasonablecomputational resources, we believe that our Eqs. (63)–(66) can still be reliably used for parameters when themaster equation already cannot be used because of toolarge Hilbert space. VII. CONCLUSION
In this paper we have introduced a new approximatemethod for numerical calculation of quantum evolutionof a weakly nonlinear resonator due to drive and dissipa-tion. This method is most accurate for large number ofphotons in the resonator (hundreds, thousands or more).This is exactly the regime where the conventional methodbased on the master equation becomes inapplicable be-cause of too large Hilbert space. For a few hundred pho-tons in the resonator (when the master equation can stillbe used), our method is faster by a factor of over 10 ,while providing a very good accuracy.The method is based on a hybrid description of a quan-tum state, which uses both phase-space and Fock-spaceparameters. The advantage is that evolution due to driveand dissipation can be naturally described in the phasespace, while evolution due to nonlinearity has a simpledescription in the Fock space. We combined both de-scriptions by proving that a phase-space Gaussian statewith many photons has a simple approximate representa-tion in the Fock space, Eq. (42), which is also Gaussian.Thus, our method essentially uses the Gaussian-state ap-proximation for an evolving quantum state. It is not ap-plicable for quantum dynamics involving cat-states, butis well-applicable for analyzing squeezing, unsqueezing,and effective heating of the resonator state due to weaknonlinearity.The method describes the quantum evolution via solv-ing four ordinary differential equations, Eqs. (63)–(66).One of them, Eq. (63), is decoupled from other equa-tions and describes the evolution of the state center β ( t )on the (complex) phase plane. This is the usual classicalequation, which takes into account resonator nonlinear-ity. (This equation can be generalized by coupling itwith other equations; however, in our numerical analy-sis we did not find a significant improvement of accuracyby doing this.) Other three equations, Eqs. (64)–(66),essentially describe evolution of the three quantum pa-rameters of a Gaussian state (maximum and minimumquadrature variances D ± b and the short-axis angle θ/ W , W , and K ). For conversion of theresults into the phase-space description we use Eqs. (45)–(47). It is also possible to use Eqs. (67)–(69) to simulateevolution of the parameters D , b , and θ directly, thoughin this paper we have not focused on this way of analy-sis. Physically, our approach is related to linearization offluctuations around a classical trajectory [29]; however,formally it is based on a different framework.Numerical accuracy of our method has been studied inSec. V. Somewhat surprisingly, it works well not only fora very large number of photons (as expected), but mayalso provide a reasonable accuracy when there are only afew dozen photons in the resonator. It is important thatthe method accurately describes the evolution startingwith vacuum (where it formally should not work); thisis because during the evolution, effects of nonlinearity8become important at larger number of photons wherethe method already works well.The method becomes inaccurate when a quantum statecannot be reasonably represented as a Gaussian state. Inour simulations this has been usually the case when thelong-axis quadrature variance D + b is large, while thenumber of photons | β | is not sufficiently large, so thatthe Wigner function of the state has a noticeable cres-cent shape in the phase plane. We have found numeri-cally that Eq. (79) can be used for a crude estimate ofthe applicability range of the method; a weaker condi-tion, | β | > [4( D + b )] , still provides a reasonably goodaccuracy. Because of a growing inaccuracy, the methodis not intended to be used close to the critical point ofthe resonator bistability, where the long-axis quadraturevariance D + b becomes large. Similarly, the methodis not intended to be used within the bistability region,since it neglects switchings between the quasistable statescaused by fluctuations. Nevertheless, the equations ofthe method can be formally used in any regime, keep-ing in mind these reasons for potential inaccuracy of theresults compared with full master-equation simulations.We have checked (Appendix C) that our analytical re-sults for the steady state agree with the results of Refs.[2] and [26].As an example, In Sec. VI the equations of our methodhave been used to derive the 3 dB limit for the steady-state squeezing of a pumped and damped weakly non-linear resonator. We have also shown numerically thatsqueezing during transients can significantly exceed this3 dB limit (Fig. 10). We emphasize that such an analysisis very difficult using the master equation because a largesqueezing typically requires large number of photons inthe resonator and therefore large Hilbert space. In con-trast, our calculations take only a fraction of a second,independently of the photon number.We hope that our method can be useful in various fieldsof research involving squeezing of weakly nonlinear res-onators with large number of quantum excitations. Inparticular, it can be useful for circuit QED systems, inwhich a weak resonator nonlinearity is induced by in-teraction with a qubit. Note that our method describessqueezing of the resonator state, but it is not directly ap-plicable to a transmitted/reflected microwave field out-side of the resonator (such generalization can be a sub-ject of future research). Our method can also be useful inanalysis of nanomechanical systems at low temperatures. ACKNOWLEDGMENTS
We would like to thank Juan Atalaya, Aashish Clerk,Justin Dressel, and Mark Dykman for useful discussions.The work was supported by ARO grant W911NF-15-1-0496.
Appendix A: Rotating-frame evolution ofa linear-resonator state
In this appendix, we discuss derivation of the rotating-frame equations (38)–(41) for evolution of the Gaussian-state parameters β , D , b , and θ from the laboratory-frame equations (32)–(36), using the rotating wave ap-proximation (RWA).Let us start with introducing the rotating frame basedon the drive frequency ω d , by defining the dimensionlessrotating-frame position and momentum operators ˆ˜ x andˆ˜ p as ˆ˜ x + i ˆ˜ p = (ˆ x + i ˆ p ) e iω d t . (A1)This is equivalent to introducing a new lowering opera-tor ˆ˜ a = ˆ a e iω d t . From Eq. (A1) we obtain the canonicaltransformation ˆ x = ˆ˜ x cos ω d t + ˆ˜ p sin ω d t, (A2)ˆ p = ˆ˜ p cos ω d t − ˆ˜ x sin ω d t. (A3)To find the rotating-frame evolution equation for theGaussian state center, we use Eqs. (32) and (33) for theevolution of x c = (cid:104) ˆ x (cid:105) and p c = (cid:104) ˆ p (cid:105) , and convert theminto equations for ˜ x c = (cid:104) ˆ˜ x (cid:105) and ˜ p c = (cid:104) ˆ˜ p (cid:105) , thus obtaining ddt (˜ x c + i ˜ p c ) = − i ( ω r − ω d )(˜ x c + i ˜ p c ) − iε − iε ∗ e i ω d t − iκ (cid:18) ˜ p c e i ω d t x c − e i ω d t i (cid:19) . (A4)This equation is still exact. Now using RWA, we neglectthe terms oscillating with frequency 2 ω d , thus obtainingslow evolution of the Gaussian state center,˙ β = − i ( ω r − ω d ) β − κ β − iε, β ≡ ˜ x c + i ˜ p c , (A5)which is Eq. (38).To derive Eqs. (39) and (40) for ˙ D and ˙ b , let us startwith expressing D x , D p , and D xp via the correspondingrotating-frame quantities D ˜ x , D ˜ p , and D ˜ x ˜ p (with obviousdefinitions) D x = D ˜ x cos ( ω d t ) + D ˜ p sin ( ω d t ) + D ˜ x ˜ p sin(2 ω d t ) , (A6) D p = D ˜ x sin ( ω d t ) + D ˜ p cos ( ω d t ) − D ˜ x ˜ p sin(2 ω d t ) , (A7) D xp = D ˜ x ˜ p cos(2 ω d t ) + (1 / D ˜ p − D ˜ x ) sin(2 ω d t ) . (A8)Note that D ≡ ( D x + D p ) / D = ( D ˜ x + D ˜ p ) /
2; similarly, b ≡ ( D p − D x ) / D xp can also be expressed as b = ( D ˜ p − D ˜ x ) / D x ˜ p .For the evolution of D , from Eqs. (34) and (35) wefind ˙ D = − κD p + ( κ/
4) coth( ω r / T b ). Then using Eq.(A7), we obtain˙ D = − κ [ D ˜ x sin ( ω d t ) + D ˜ p cos ( ω d t ) − D ˜ x ˜ p sin(2 ω d t )]+ ( κ/
4) coth( ω r / T b ) . (A9)9Now using RWA, we neglect the terms oscillating withfrequency 2 ω d , so that sin ( ω d t ) → /
2, cos ( ω d t ) → /
2, and sin(2 ω d t ) →
0. This gives us˙ D = − κD + ( κ/
4) coth( ω r / T b ) , (A10)which is Eq. (39).For the evolution of b , from Eqs. (34)–(36) we obtain d ( b ) /dt = ( κ/ D p − D x ) coth( ω r / T b ) − κ ( D p − D x ) D p − κD xp . (A11)Within RWA, the first term on the right-hand side is zerobecause D p − D x oscillates with frequency 2 ω d [see Eqs.(A6) and (A7)]. The second term is not zero because D p has also a part oscillating with 2 ω d ; averaging over theseoscillations we obtain − κ [ D x ˜ p + ( D ˜ p − D ˜ x ) / − κb . Similarly, for the third term we use Eq.(A8) and averaging over the oscillations obtain − κ [ D x ˜ p +( D ˜ p − D ˜ x ) / − κb . Thus, within RWA d ( b ) /dt = − κb . (A12)Equivalently, ˙ b = − κb , which is Eq. (40).To derive Eq. (41) for ˙ θ , we start with Eqs. (14) and(31), which give θ = arctan (cid:18) D xp D x − D p (cid:19) + 2 ω d t + ( π/ D x − D p )] . (A13)Neglecting the last term, the time derivative is˙ θ = ˙ D xp ( D x − D p ) − D xp ( ˙ D x − ˙ D p )2 b + 2 ω d . (A14)Using Eqs. (34)–(36), we find that the numerator here is − ω r b − κD xp D +( κ/ D xp coth( ω r / T b ), in which theonly non-oscillating term is − ω r b . Dividing it by 2 b and adding 2 ω d , from Eq. (A14) we obtain ˙ θ = − ω r − ω d ), which is Eq. (41). Appendix B: Equivalence between Gaussian andFock-space Gaussian states
In this appendix, we show that the Fock-space Gaus-sian state introduced in Eq. (42) is approximately thesame as the standard Gaussian state [Eq. (9)] in thelimit of large photon number, | β | (cid:29)
1, and derive theconversion relations (44)–(47). This is done by compar-ing the Husimi Q -functions of the Gaussian and Fock-space Gaussian states. We use the rotating frame andcharacterize the Gaussian state by the complex param-eter β (center) and three real parameters: D , b , and θ – see Eqs. (12)–(14). The Fock-space Gaussian state ischaracterized by the complex parameter e iφ β | β | (whichis chosen to be the same as β ) and three real parameters: W , W , and K – see Eq. (42). The Husimi Q -function Q ( α ) of a state with densitymatrix ρ is defined via its overlap with the coherent state | α (cid:105) , Q ( α ) = 1 π (cid:104) α | ρ | α (cid:105) , | α (cid:105) = e − | α | ∞ (cid:88) n =0 α n √ n ! | n (cid:105) , (B1)where α = ˜ x + i ˜ p assumes the rotating frame, in contrastto the notation α used in Sec. III A. The function Q ( α )can be calculated from the Wigner function W ( α ) (herein the rotating frame; note a slightly different notationused in Sec. III A), Q ( α ) = 2 π (cid:90) W ( α (cid:48) ) e − | α − α (cid:48) | d Re( α (cid:48) ) d Im( α (cid:48) ) . (B2)For the Gaussian state (9) it is equal Q ( α ) = π − [4( D − b + 1 / D + b + 1 / − / × exp (cid:26) − ( D + b cos θ + 1 /
4) [Re( α − β )] D − b + 1 / D + b + 1 / − ( D − b cos θ + 1 /
4) [Im( α − β )] D − b + 1 / D + b + 1 / − (2 b sin θ ) Re( α − β ) Im( α − β )2( D − b + 1 / D + b + 1 / (cid:27) . (B3)Recall that β is the Gaussian state center, D + b is themaximum quadrature variance, D − b is the minimumquadrature variance, and θ/ x -axis (see Fig. 1).Note that in the diagonal basis, Eq. (B3) reduces to Eq.(17), up to a slight change of notations.Now let us calculate the Q -function for the Fock-spaceGaussian state, Eq. (42), and compare it with Eq. (B3).We will use a series of approximations to calculate Q ( α ).First, for | β | (cid:29) | α | (cid:29)
1; then thecoherent state | α (cid:105) in Eq. (B1) can be approximated as | α (cid:105) ≈ (2 π | α | ) − / (cid:80) n exp[ − ( n − | α | ) / | α | ] exp[ inφ α ],where φ α = arg( α ), so that the Q -function is approxi-mately Q ( α ) = 1 π (cid:112) π | α | ∞ (cid:88) n,m =0 ρ nm exp (cid:20) − ( n − | α | ) | α | − ( m − | α | ) | α | − iφ α ( n − m ) (cid:21) . (B4)0Substituting ρ nm from Eq. (42), we obtain Q ( α ) = N (cid:88) n,m exp[ − An − ˜ Am − B ( m ) n − ˜ Bm − C ] , (B5) N = π − (4 π W | β | | α | ) − / , (B6) A = 14 | α | + 18 W | β | + 18 W | β | + i K | β | , (B7)˜ A = 14 | α | + 18 W | β | + 18 W | β | − i K | β | , (B8) B ( m ) = −
12 + i ( φ α − φ β ) + m W | β | − W − m W | β | − iK, (B9)˜ B = − − i ( φ α − φ β ) − W + 2 iK, (B10) C = | α | | β | W . (B11)Then replacing summation over n and m by integrationwithin infinite limits (assuming | β | (cid:29)
1) and calculatingthe integral over n , we find Q ( α ) = N √ πe − C √ A ∞ (cid:90) −∞ exp (cid:20) [ B ( m )] A − ˜ Am − ˜ Bm (cid:21) dm. (B12)Using Eq. (B9), we then represent [ B ( m )] / A as[ B ( m )] / A = ¯ Am + ¯ Bm + ¯ C, (B13)¯ A = 14 A (cid:16) W | β | (cid:17) (cid:16) − W W (cid:17) , (B14)¯ B = 14 A W | β | (cid:16) − W W (cid:17)(cid:104) − − W + i ( φ α − φ β ) − iK (cid:105) , (B15)¯ C = 14 A (cid:16) − − W + i ( φ α − φ β ) − iK (cid:17) . (B16)Then the exponent in Eq. (B12) is exp[ − ( ˜ A − ¯ A ) m − ( ˜ B − ¯ B ) m ] and its integral over dm can be easily calculated, Q ( α ) = (2 π (cid:112) W | β || α | ) − [ A ( ˜ A − ¯ A )] − / × exp { ( ˜ B − ¯ B ) / [4( ˜ A − ¯ A )] − C − ¯ C } . (B17)Since we want to compare this result with Eq. (B3),we need to find its dependence on the difference α − β .Assuming | β | (cid:29)
1, we expand Eq. (B17) up to secondorder in Re( α − β ) and Im( α − β ). Let us consider first thespecial case when β is real ( β > φ β = 0. Thenexpansion of Eq. (B17) produces (after some algebra) the result Q ( α ) ≈ π (cid:16)
14 + W W + 14 W + W K W (cid:17) − / × exp (cid:26) − W + 16 W W K ) [Re( α − β )] W + W + W W (1 + 16 K ) − W (1 + W ) [Im( α − β )] W + W + W W (1 + 16 K ) − W W K Re( α − β ) Im( α − β )1 + W + W + W W (1 + 16 K ) (cid:27) . (B18)Comparing this formula with Eq. (B3) for the Gaussianstate, we see that the formulas coincide if D = 18 W + W K W , (B19) b = 14 (cid:114)(cid:16) W + W K W (cid:17) − W W , (B20) θ = arctan (cid:16) KW W − W W + 16 K W W (cid:17) + ( π/
2) [1 − sign(1 − W W + 16 K W W )] , (B21)where we use notation θ instead of θ to remind that weconsider the special case of a real positive β . Note thatEqs. (B19)–(B21) coincide with Eqs. (45)–(47) in the caseof a real positive β .For a complex β , it is also possible to use the second-order expansion of Eq. (B17); however, it is easier to usethe fact that dependence of Q ( α ) on the complex phase φ β in Eqs. (B5)–(B11) comes only from the combination φ α − φ β . Therefore, the Q -function of the Fock-spaceGaussian state does not change in the transformation β → | β | , α → e − iφ β α , so for a complex β we can still useEq. (B18) with the substitution ( α − β ) → e − iφ β ( α − β ).Using this substitution in the equivalent Eq. (B3), weeasily find that it results in replacing the angle θ (forreal β ) with θ = θ + 2 φ β , (B22)while the parameters D and b do not change. Anotherway to obtain Eq. (B22) is to note that the parameters W , W , and K of the Fock-space Gaussian state do notchange when the phase space is rotated (i.e., β → e i ∆ φ β , α → e i ∆ φ α ), while for the Gaussian state this results inthe change θ → θ + 2∆ φ with unchanged parameters D and b (see Fig. 1). Therefore, we can first rotate the phasespace clockwise by the angle φ β (to make β real), thenconvert parameters W , W , and K , into D , b , and θ using Eqs. (B19)–(B21), and then move the phase spaceback by counterclockwise rotation with the same angle φ β , which results in θ change (B22).Thus we have derived the conversion relations (45)–(47) between the Gaussian and Fock-space Gaussianstates ( β does not change). Note that our derivationrelied on the fact that the Husimi Q -function uniquelydefines a quantum state [82]. Since Eq. (B18) is only an1approximation, a Fock-space Gaussian state is not ex-actly equal to a Gaussian state. However, the accuracyof the conversion improves at larger | β | , approaching ex-act equivalence in the limit | β | → ∞ . Numerical resultsin Sec. V A show that infidelity of the conversion scalesas | β | − . Appendix C: Steady-state squeezing and heating
In this appendix, we derive results for D , b , and θ in the steady state. The parameters r and n th can bethen calculated using Eq. (30). The squeezing factoris [4( D − b )] − , the effective temperature T eff is givenby coth( ω r0 / T eff ) = 4 (cid:112) ( D + b )( D − b ). All variablesdiscussed in this appendix are only for the steady state.The steady-state value of β can be calculated from Eq.(63); in general it does not have an analytical expression.Note that ε/β = ω d − ω r ( | β | ) + iκ/ , (C1)so Re( ε/β ) can be positive or negative, depending ondetuning.From Eqs. (67) and (68) in the steady state we find D = coth( ω r0 / T b )4 11 − [2 η β | β | sin(∆ θ ) /κ ] , (C2) b = coth( ω r0 / T b )4 2 η β | β | sin(∆ θ ) /κ − [2 η β | β | sin(∆ θ ) /κ ] , (C3)where η β = dω r ( n ) /dn | n = | β | is the steady-state nonlin-earity. To obtain explicit analytics for D and b , we stillneed to find sin(∆ θ ). For that we can substitute the ra-tio b/D = 2 η β | β | sin(∆ θ ) /κ into Eq. (69) in the steadystate, thus obtainingtan(∆ θ ) = κ/ η β | β | − Re( ε/β ) . (C4)Since η β sin(∆ θ ) ≥ b ≥ θ ) = sign( η β ) (cid:115) ( κ/ ( κ/ + [ η β | β | − Re( ε/β )] (C5)in Eqs. (C2) and (C3).The angle θ can be calculated as θ = 2 arg( β ) + arctan (cid:18) κ/ η β | β | − Re( ε/β ) (cid:19) + ( π/ { − sign[ | β | − η − β Re( ε/β )] } . (C6)These results can be compared with results of Ref.[2] in the case of Kerr nonlinearity (Duffing oscillator), H lfr = ω r0 a † a + ( η/ a † ) a , which is equivalent to ourHamiltonian when ω r = ω r0 + nη [see Eq. (73)], so that η β = η = const. In this case Eq. (4.4) of Ref. [2] (con-verted into our notations) gives (cid:104) a (cid:105) − β = − ηβ ( ω r0 + 2 η | β | − ω d + iκ/ n b )2 λ , (C7) (cid:104) a † a (cid:105) − | β | = η | β | (1 + 2 n b )2 λ + n b , (C8) λ = ( ω r0 + 2 η | β | − ω d ) + κ / − η | β | . (C9)From these values, D , b , and θ can be obtained usingEqs. (18)–(20) [also, Eq. (26) gives (cid:104) a † a (cid:105) − | β | = 2 D − / κ →
0) can also be directlycompared with the analytical results presented in Secs.2.1 and 2.5 of Ref. [26]. In this case the squeezing andheating are determined only by the parameter combina-tion ε η/ ( ω d − ω r0 ) (which was called β in Ref. [26]).Results of Ref. [26] show that the squeezing parameter ξ = re iθ is real and equals ξ = 14 ln 3 Q − Q − , (C10)where Q satisfies equation Q ( Q −
1) = (cid:112) ε η/ ( ω d − ω r0 ) . (C11)Here in the case ε η/ ( ω d − ω r0 ) > /
27, there is onlyone real solution for Q . The range 0 < ε η/ ( ω d − ω r0 ) < /
27 corresponds to bistability, and there are three realsolutions for Q , with the largest value corresponding tothe upper bistability branch and the middle value for thelower branch. In the case ε η/ ( ω d − ω r0 ) <
0, we needto use the purely imaginary solution for Q .The angle θ in this limit is zero (squeezing is in phasewith the drive), except θ = π for the lower bistabilitybranch (then ξ < n th = n b + (2 n b + 1) sinh r. (C12)We have numerically compared Eqs. (C10) and (C12) for ξ and n th with our results following from Eqs. (C2), (C3),and (C6). As expected, we have found that they coincidein the limit κ → ε η/ ( ω d − ω r0 ) .Thus, our results agree with the results of Ref. [26].2 [1] M. I. Dykman and M. A. Krivoglaz, “Quantum theory ofnonlinear oscillators interacting with the medium,” Sov.Phys. JETP , 506 (1973).[2] P. D. Drummond and D. F. Walls, “Quantum theory ofoptical bistability. I. Nonlinear polarisability model,” J.Phys. A: Math. Gen. , 725 (1980).[3] R. Tana´s, “Squeezed states of an anharmonic oscillator,”in Coherence and Quantum Optics V , Proceedings of thefifth Rochester conference on coherence and quantum op-tics, edited by L. Mandel and E. Wolf (Plenum Press,New York, 1984) p. 645.[4] L. A. Lugiato, “Theory of optical bistability,” Prog. Opt. , 69 (1984).[5] Fluctuating nonlinear oscillators , edited by M. I. Dyk-man (Oxford University Press, Oxford, 2012).[6] A. Wallraff, D. I. Schuster, A. Blais, L. Frunzio, R.-S.Huang, J. Majer, S. Kumar, S. M. Girvin, and R. J.Schoelkopf, “Strong coupling of a single photon to asuperconducting qubit using circuit quantum electrody-namics,” Nature (London) , 162 (2004).[7] I. Siddiqi, R. Vijay, F. Pierre, C. M. Wilson, L. Frunzio,M. Metcalfe, C. Rigetti, R. J. Schoelkopf, M. H. Devoret,D. Vion, and D. Esteve, “Direct observation of dynami-cal bifurcation between two driven oscillation states of aJosephson junction,” Phys. Rev. Lett. , 027005 (2005).[8] M. Blencowe, “Quantum electromechanical systems,”Phys. Rep. , 159 (2004).[9] K. C. Schwab and M. L. Roukes, “Putting mechanics intoquantum mechanics,” Phys. Today , 36 (2005).[10] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bial-czak, M. Lenander, E. Lucero, M. Neeley, D. Sank,H. Wang, M. Weides, J. Wenner, J. M. Martinis, andA. N. Cleland, “Quantum ground state and single-phonon control of a mechanical resonator,” Nature (Lon-don) , 697 (2010).[11] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh,A. Kronwald, F. Marquardt, A. A. Clerk, and K. C.Schwab, “Quantum squeezing of motion in a mechanicalresonator,” Science , 952 (2015).[12] A. Blais, R.-S. Huang, A. Wallraff, S. M. Girvin, andR. J. Schoelkopf, “Cavity quantum electrodynamics forsuperconducting electrical circuits: An architecture forquantum computation,” Phys. Rev. A , 062320 (2004).[13] J. M. Chow, J. M. Gambetta, E. Magesan, D. W. Abra-ham, A. W. Cross, B. R. Johnson, N. A. Masluk, C. A.Ryan, J. A. Smolin, S. J. Srinivasan, and M. Steffen,“Implementing a strand of a scalable fault-tolerant quan-tum computing fabric,” Nat. Commun. , 4015 (2014).[14] E. Jeffrey, D. Sank, J. Y. Mutus, T. C. White, J. Kelly,R. Barends, Y. Chen, Z. Chen, B. Chiaro, A. Dunsworth,A. Megrant, P. J. J. O’Malley, C. Neill, P. Roushan,A. Vainsencher, J. Wenner, A. N. Cleland, and J. M.Martinis, “Fast accurate state measurement with super-conducting qubits,” Phys. Rev. Lett. , 190504 (2014).[15] O.-P. Saira, J. P. Groen, J. Cramer, M. Meretska,G. de Lange, and L. DiCarlo, “Entanglement genesis byancilla-based parity measurement in 2D circuit QED,”Phys. Rev. Lett. , 070502 (2014).[16] M. Boissonneault, J. M. Gambetta, and A. Blais, “Im-proved superconducting qubit readout by qubit-inducednonlinearities,” Phys. Rev. Lett. , 100504 (2010). [17] M. D. Reed, L. DiCarlo, B. R. Johnson, L. Sun, D. I.Schuster, L. Frunzio, and R. J. Schoelkopf, “High-fidelityreadout in circuit quantum electrodynamics using theJaynes-Cummings nonlinearity,” Phys. Rev. Lett. ,173601 (2010).[18] M. Khezri, E. Mlinar, J. Dressel, and A. N. Korotkov,“Measuring a transmon qubit in circuit QED: Dressedsqueezed states,” Phys. Rev. A , 012347 (2016).[19] R. Vijay, Josephson bifurcation amplifier: Amplifyingquantum signals using a dynamical bifurcation , Ph.D.thesis, Yale University (2008).[20] N. Bergeal, F. Schackert, M. Metcalfe, R. Vijay, V. E.Manucharyan, L. Frunzio, D. E. Prober, R. J. Schoelkopf,S. M. Girvin, and M. H. Devoret, “Phase-preserving am-plification near the quantum limit with a Josephson ringmodulator,” Nature (London) , 64 (2010).[21] F. Mallet, M. A. Castellanos-Beltran, H. S. Ku,S. Glancy, E. Knill, K. D. Irwin, G. C. Hilton, L. R.Vale, and K. W. Lehnert, “Quantum state tomogra-phy of an itinerant squeezed microwave field,” Phys. Rev.Lett. , 220502 (2011).[22] R. Vijay, D. H. Slichter, and I. Siddiqi, “Observationof quantum jumps in a superconducting artificial atom,”Phys. Rev. Lett. , 110502 (2011).[23] R. Tana´s, “Squeezing from an anharmonic oscillatormodel: ( a † ) a versus ( a † a ) interaction Hamiltonians,”Phys. Lett. A , 217 (1989).[24] G. J. Milburn, “Quantum and classical Liouville dynam-ics of the anharmonic oscillator,” Phys. Rev. A , 674(1986).[25] M. Kitagawa and Y. Yamamoto, “Number-phaseminimum-uncertainty state with reduced number uncer-tainty in a Kerr nonlinear interferometer,” Phys. Rev. A , 3974 (1986).[26] M. I. Dykman, “Periodically modulated quantum non-linear oscillators,” in Fluctuating nonlinear oscilla-tors (Oxford University Press, Oxford, 2012) p. 165,arXiv:1112.2407.[27] K. Tomita and H. Tomita, “Irreversible circulation offluctuation,” Prog. Theor. Phys. , 1731 (1974).[28] M. I. Dykman and M. A. Krivoglaz, “Theory of fluctua-tional transitions between the stable states of a non-linearoscillator,” Sov. Phys. JETP , 30 (1979).[29] D. Ludwig, “Persistence of dynamical systems under ran-dom perturbations,” SIAM Rev. , 605 (1975).[30] D. F. Walls and G. J. Milburn, Quantum optics (Springer, Berlin, 2008).[31]
Quantum squeezing , edited by P. D. Drummond and Z.Ficek (Springer-Verlag, Berlin, 2004).[32] M. I. Dykman, “Critical exponents in metastable de-cay via quantum activation,” Phys. Rev. E , 011101(2007).[33] C. Laflamme and A. A. Clerk, “Quantum-limited ampli-fication with a nonlinear cavity detector,” Phys. Rev. A , 033803 (2011).[34] C. M. Caves, “Quantum-mechanical noise in an interfer-ometer,” Phys. Rev. D , 1693 (1981).[35] V. Giovannetti, S. Lloyd, and L. Maccone, “Quantum-enhanced measurements: Beating the standard quantumlimit,” Science , 1330 (2004).[36] LIGO Scientific Collaboration, “Enhanced sensitivity of the LIGO gravitational wave detector by using squeezedstates of light,” Nat. Photon. , 613 (2013).[37] E. A. Sete, A. Galiautdinov, E. Mlinar, J. M. Martinis,and A. N. Korotkov, “Catch-Disperse-Release readout forsuperconducting qubits,” Phys. Rev. Lett. , 210501(2013).[38] S. Barzanjeh, D. P. DiVincenzo, and B. M. Terhal, “Dis-persive qubit measurement by interferometry with para-metric amplifiers,” Phys. Rev. B , 134515 (2014).[39] N. Didier, A. Kamal, W. D. Oliver, A. Blais, and A. A.Clerk, “Heisenberg-limited qubit read-out with two-modesqueezed light,” Phys. Rev. Lett. , 093604 (2015).[40] N. Didier, J. Bourassa, and A. Blais, “Fast quantumnondemolition readout by parametric modulation of lon-gitudinal qubit-oscillator interaction,” Phys. Rev. Lett. , 203601 (2015).[41] L. C. G. Govia and A. A. Clerk, “Enhanced qubit read-out using locally generated squeezing and inbuilt Purcell-decay suppression,” New J. Phys. , 023044 (2017).[42] A. Eddins, S. Schreppler, D. M. Toyli, L. S. Martin,S. Hacohen-Gourgy, L. C. G. Govia, H. Ribeiro, A. A.Clerk, and I. Siddiqi, “Stroboscopic qubit measurementwith squeezed illumination,” arXiv:1708.01674 .[43] R. Ruskov, K. Schwab, and A. N. Korotkov, “Squeez-ing of a nanomechanical resonator by quantum nonde-molition measurement and feedback,” Phys. Rev. B ,235407 (2005).[44] V. Peano, H. G. L. Schwefel, Ch. Marquardt,and F. Marquardt, “Intracavity squeezing can en-hance quantum-limited optomechanical position detec-tion through deamplification,” Phys. Rev. Lett. ,243603 (2015).[45] R. Movshovich, B. Yurke, P. G. Kaminsky, A. D. Smith,A. H. Silver, R. W. Simon, and M. V. Schneider, “Ob-servation of zero-point noise squeezing via a Josephson-parametric amplifier,” Phys. Rev. Lett. , 1419 (1990).[46] C. Eichler, D. Bozyigit, C. Lang, M. Baur, L. Steffen,J. M. Fink, S. Filipp, and A. Wallraff, “Observationof two-mode squeezing in the microwave frequency do-main,” Phys. Rev. Lett. , 113601 (2011).[47] E. Flurin, N. Roch, F. Mallet, M. H. Devoret, andB. Huard, “Generating entangled microwave radiationover two transmission lines,” Phys. Rev. Lett. ,183901 (2012).[48] K. W. Murch, S. J. Weber, K. M. Beck, E. Ginossar, andI. Siddiqi, “Reduction of the radiative decay of atomiccoherence in squeezed vacuum,” Nature (London) ,62 (2013).[49] D. M. Toyli, A. W. Eddins, S. Boutin, S. Puri, D. Hover,V. Bolkhovsky, W. D. Oliver, A. Blais, and I. Sid-diqi, “Resonance fluorescence from an artificial atom insqueezed vacuum,” Phys. Rev. X , 031004 (2016).[50] G. Kirchmair, B. Vlastakis, Z. Leghtas, S. E. Nigg,H. Paik, E. Ginossar, M. Mirrahimi, L. Frunzio, S. M.Girvin, and R. J. Schoelkopf, “Observation of quantumstate collapse and revival due to the single-photon Kerreffect,” Nature (London) , 205 (2013).[51] G. Milburn and D. F. Walls, “Production of squeezedstates in a degenerate parametric amplifier,” Opt. Com-mun. , 401 (1981).[52] M. J. Collett and C. W. Gardiner, “Squeezing of intracav-ity and traveling-wave light fields produced in parametricamplification,” Phys. Rev. A , 1386 (1984).[53] D. Rugar and P. Gr¨utter, “Mechanical parametric am- plification and thermomechanical noise squeezing,” Phys.Rev. Lett. , 699 (1991).[54] B. Yurke, “Use of cavities in squeezed-state generation,”Phys. Rev. A , 408 (1984).[55] P. Rabl, A. Shnirman, and P. Zoller, “Generation ofsqueezed states of nanomechanical resonators by reser-voir engineering,” Phys. Rev. B , 205304 (2004).[56] A. Kronwald, F. Marquardt, and A. A. Clerk, “Arbitrar-ily large steady-state bosonic squeezing via dissipation,”Phys. Rev. A , 063833 (2013).[57] A. Szorkovszky, A. C. Doherty, G. I. Harris, and W. P.Bowen, “Mechanical squeezing via parametric amplifi-cation and weak measurement,” Phys. Rev. Lett. ,213603 (2011).[58] K. J¨ahne, C. Genes, K. Hammerer, M. Wallquist, E. S.Polzik, and P. Zoller, “Cavity-assisted squeezing of amechanical oscillator,” Phys. Rev. A , 063819 (2009).[59] M. R. Vanner, I. Pikovski, G. D. Cole, M. S. Kim,ˇC Brukner, K. Hammerer, G. J. Milburn, and M. As-pelmeyer, “Pulsed quantum optomechanics,” Proc. Natl.Acad. Sci. , 16182 (2011).[60] C. U. Lei, A. J. Weinstein, J. Suh, E. E. Wollman,A. Kronwald, F. Marquardt, A. A. Clerk, and K. C.Schwab, “Quantum nondemolition measurement of aquantum squeezed state beyond the 3 dB limit,” Phys.Rev. Lett. , 100801 (2016).[61] M. I. Dykman, D. G. Luchinsky, R. Mannella, P. V. E.McClintock, N. D. Stein, and N. G. Stocks, “Supernar-row spectral peaks and high-frequency stochastic reso-nance in systems with coexisting periodic attractors,”Phys. Rev. E , 1198 (1994).[62] E. Buks and B. Yurke, “Mass detection with a nonlin-ear nanomechanical resonator,” Phys. Rev. E , 046619(2006).[63] R. Almog, S. Zaitsev, O. Shtempluck, and E. Buks,“Noise squeezing in a nanomechanical Duffing res-onator,” Phys. Rev. Lett. , 078103 (2007).[64] I. Serban, M. I. Dykman, and F. K. Wilhelm, “Relax-ation of a qubit measured by a driven Duffing oscillator,”Phys. Rev. A , 022305 (2010).[65] A. N. Cleland, Foundations of nanomechanics: fromsolid-state theory to device applications (Springer, Berlin,2003).[66] C. Gardiner and P. Zoller,
Quantum noise (Springer,Berlin, 2004).[67] C. W. Gardiner and M. J. Collett, “Input and output indamped quantum systems: Quantum stochastic differen-tial equations and the master equation,” Phys. Rev. A , 3761 (1985).[68] B. Yurke and J. S. Denker, “Quantum network theory,”Phys. Rev. A , 1419 (1984).[69] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt,and R. J. Schoelkopf, “Introduction to quantum noise,measurement, and amplification,” Rev. Mod. Phys. ,1155 (2010).[70] H. M. Wiseman and G. J. Milburn, Quantum measure-ment and control (Cambridge University Press, Cam-bridge, UK, 2010).[71] A. N. Korotkov, “Quantum Bayesian approach to circuitQED measurement with moderate bandwidth,” Phys.Rev. A , 042326 (2016).[72] D. Sank, Z. Chen, M. Khezri, J. Kelly, R. Barends,B. Campbell, Y. Chen, B. Chiaro, A. Dunsworth,A. Fowler, E. Jeffrey, E. Lucero, A. Megrant, J. Mu- tus, M. Neeley, C. Neill, P. J. J. O’Malley, C. Quintana,P. Roushan, A. Vainsencher, T. White, J. Wenner, A. N.Korotkov, and J. M. Martinis, “Measurement-inducedstate transitions in a superconducting qubit: Beyondthe rotating wave approximation,” Phys. Rev. Lett. ,190503 (2016).[73] B. Ya. Zel’dovich, A. M. Perelomov, and V. S. Popov,“Relaxation of a quantum oscillator,” Sov. Phys. JETP , 308 (1969).[74] U. Weiss, Quantum dissipative systems (World Scientific,London, 2012).[75] A. N. Korotkov, “Error matrices in quantum process to-mography,” arXiv:1309.6405 .[76] C. C. Bultink, M. A. Rol, T. E. O’Brien, X. Fu, B. C. S.Dikken, C. Dickel, R. F. L. Vermeulen, J. C. de Sterke,A. Bruno, R. N. Schouten, and L. DiCarlo, “Active res-onator reset in the nonlinear dispersive regime of circuitQED,” Phys. Rev. Applied , 034008 (2016).[77] J. Halliwell and A. Zoupas, “Quantum state diffusion,density matrix diagonalization, and decoherent histories:A model,” Phys. Rev. D , 7294 (1995).[78] W. H. Zurek, S. Habib, and J. P. Paz, “Coherent statesvia decoherence,” Phys. Rev. Lett. , 1187 (1993).[79] C. Weedbrook, S. Pirandola, R. Garc´ıa-Patr´on, N. J.Cerf, T. C. Ralph, J. H. Shapiro, and S. Lloyd, “Gaus-sian quantum information,” Rev. Mod. Phys. , 621(2012).[80] S. L. Braunstein and P. van Loock, “Quantum informa-tion with continuous variables,” Rev. of Mod. Phys. ,513 (2005).[81] A. Ferraro, S. Olivares, and M. G. A. Paris, Gaus-sian states in quantum information (Bibliopolis, Napoli,2005); arXiv:quant-ph/0503237 .[82] C. C. Gerry and P. Knight,
Introductory quantum optics (Cambridge University Press, Cambridge, 2005).[83] P. Marian and T. A. Marian, “Squeezed states with ther- mal noise. I. Photon-number statistics,” Phys. Rev. A ,4474 (1993).[84] H. Fearn and M. J. Collett, “Representations of squeezedstates with thermal noise,” J. Mod. Opt. , 553 (1988).[85] M. S. Kim, F. A. M. de Oliveira, and P. L. Knight,“Properties of squeezed number states and squeezed ther-mal states,” Phys. Rev. A , 2494 (1989).[86] A. C. Doherty and K. Jacobs, “Feedback control of quan-tum systems using continuous state estimation,” Phys.Rev. A , 2700 (1999).[87] A. Hopkins, K. Jacobs, S. Habib, and K. Schwab, “Feed-back cooling of a nanomechanical resonator,” Phys. Rev.B , 235328 (2003).[88] M. G. A. Paris, F. Illuminati, A. Serafini, andS. De Siena, “Purity of Gaussian states: Measurementschemes and time evolution in noisy channels,” Phys.Rev. A , 012314 (2003).[89] A. Serafini, M. G. A. Paris, F. Illuminati, andS. De Siena, “Quantifying decoherence in continuousvariable systems,” J. Opt. B: Quantum Semiclass. Opt. , R19 (2005).[90] M. A. Nielsen and I. L. Chuang, Quantum computationand quantum information (Cambridge University Press,Cambridge, UK, 2000).[91] K. E. Cahill and R. J. Glauber, “Density operators andquasiprobability distributions,” Phys. Rev. , 1882(1969).[92] S. Haroche and J.-M. Raimond,
Exploring the quantum:Atoms, cavities, and photons (Oxford University Press,New York, 2006).[93] L. D. Landau and E. M. Lifshitz,
Mechanics (Butterworth-Heinemann, Amsterdam, 1976), Sec. 29.[94] A. H. Nayfeh and D. T. Mook,