The dissipative quantum Duffing oscillator: a comparison of Floquet-based approaches

# The dissipative quantum Duffing oscillator: a comparison of Floquet-based approaches

Carmen Vierheilig and Milena Grifoni Institut für Theoretische Physik, Universität Regensburg
###### Abstract

We study the dissipative quantum Duffing oscillator in the deep quantum regime with two different approaches: The first is based on the exact Floquet states of the linear oscillator and the nonlinearity is treated perturbatively. It well describes the nonlinear oscillator dynamics away from resonance. The second, in contrast, is applicable at and in the vicinity of a -photon resonance and it exploits quasi-degenerate perturbation theory for the nonlinear oscillator in Floquet space. It is perturbative both in driving and nonlinearity. A combination of both approaches yields the possibility to cover the whole range of driving frequencies. As an example we discuss the dissipative dynamics of the Duffing oscillator near and at the one-photon resonance.

###### keywords:
Nonlinear quantum oscillator, Floquet theory
journal: Chemical Physics

## 1 Introduction

Classical nonlinear systems show interesting phenomena like bistability, frequency doubling and nonlinear response and cover a wide range of applicability and various physical realizations Nayfeh (); Jordan (). In the last years it has been possible to build nonlinear devices which can potentially reach the quantum regime. These are for example cavities incorporating a Josephson junction Boaknin (); Metcalfe (), SQUIDs used as bifurcation amplifiers to improve qubit read-out Lee (); Picot (); Siddiqi (); Siddiqi2 () or nanomechanical resonators Almog (); Alridge (). Recently, a novel class of devices combining SQUIDS and resonators has been demonstrated. For example, sensitive detection of the position of a micromechanical resonator embedded in a nonlinear, strongly damped DC-SQUID has been achieved Etaki (). In the deep quantum regime, where bistability is no longer observed, there has been to date no experimental operation of nonlinear driven oscillators to our knowledge.
From the theoretical side, semiclassical approaches have been used to describe the situation where the underlying classical bistability still plays a dominant role. A DC-SQUID embedded in a cavity allowing displacement detection and cooling was analyzed in Nation (). Composed qubit-Josephson bifurcation amplifier systems have been considered in SerbanDet (); Serbanqubit (). Dynamical tunneling in a Duffing oscillator was accounted for in SerbanMQDT () within a semiclassical WKB scheme, while in Guo (); Katz1 (); Katz2 () a Wigner function analysis near the bifurcation point is put forward.
The behaviour of the Duffing oscillator (DO) in the deep quantum regime has attracted recently lot of interest. In particular Rigo et al. Rigo () demonstrated, based on a quantum diffusion model, that in the steady state the quantum DO does not exhibit any bistability or hysteresis. It was also shown that the response of the Duffing oscillator displays antiresonant dips and resonant peaks Fistul (); Peano1 (); Peano2 (); Peano3 () depending on the frequency of the driving field, originating from special degeneracies of the eigenenergy spectrum of the nonlinear oscillator Fistul (). While the antiresonances persist in the presence of a weak Ohmic bath, for strong damping the nonlinear response turns to a resonant behaviour, namely the one of a linear oscillator at a shifted frequency Peano2 (). Finally, recently Nakano et al. Nakano () looked at the composed qubit-DO dynamics during read-out process.
In this work we investigate the deep quantum limit of the quantum Duffing oscillator and present two different approaches covering different parameter regimes. The first approach is based on the exact Floquet energies and states of the driven linear oscillator with the nonlinearity treated perturbatively. As there is no restriction on the driving amplitude, this scheme can also be applied to the regime where the driving amplitude is larger than the nonlinearity. The second approach treats both the driving and the nonlinearity perturbatively. It is applicable for driving frequencies which can resonantly excite two states of the nonlinear oscillator, requiring that the driving cannot overcome the nonlinearity. In general a combination of both approaches allows to cover the whole range of driving frequencies. Exemplarily we consider the dynamics of the Duffing oscillator near the one-photon resonance, where the oscillator dynamics is described analytically. As in Fistul (); Peano1 (); Peano2 (); Peano3 () we obtain that for weak dissipation the amplitude of the oscillations displays an antiresonance rather than a resonance. We find a characteristic asymmetry of the antiresonance lineshape. In contrast to Fistul (); Peano1 (); Peano2 (); Peano3 (), our analytic results are obtained without applying a rotating wave approximation (RWA) on the Duffing oscillator.
The paper is organized as follows: In section 2 we introduce the Hamiltonian of the non-dissipative Duffing oscillator and the two Floquet based approximation schemes to treat it. The energy spectrum and eigenstates of the non-dissipative system are calculated with the two different schemes in section 3 and section 4 and both approaches are compared in section 5. Afterwards dissipative effects are included within a Born-Markov-Floquet master equation in section 6. Section 7 addresses the special case of the one-photon resonance including dissipative effects. In section 8 conclusions are drawn.

## 2 Quantum Duffing oscillator

A quantum Duffing oscillator is described by the Hamiltonian:

 ^HDO(t)=^P2y2M+MΩ22^y2+α4^y4+^yFcos(ωext), (1)

where and are the mass and frequency of the Duffing oscillator which is driven by a monochromatic field of amplitude and frequency . For later convenience we introduce the oscillator length . In the following we will consider the case of hard nonlinearities, , such that the undriven potential is monostable.
To treat the quantum Duffing oscillator problem we observe that the Hamiltonian can be rewritten as:

 ^HDO(t) = ^HLO(t)+α4^y4 (2a) = ^HNLO+^yFcos(ωext), (2b)

where describes a driven linear oscillator, while is the Hamiltonian of an undriven nonlinear oscillator. Due to the periodic driving Floquet theory Floquet (); Shirley (); Sambe (); Grifoni304 (), reviewed in A, can be applied. In particular Floquet theorem states that solutions of the time-dependent Schrödinger equation for are of the form

 |ψj(t)⟩=exp(−iϵjt/ℏ)|ϕj(t)⟩, (3)

where . The quasienergies and Floquet states , respectively, are eigenvalues and eigenfunctions of the Floquet Hamiltonian . As discussed in A, yields a physically equivalent solution but with shifted quasienergy . Eqs. (2a) and (2b) suggest two different approaches, shown in Figure 1, to solve the eigenvalue problem described by Eqs. (91) and (92):

In the first one, called App I, starting point are the exact Floquet states and eigenenergies of the driven linear oscillator , see Eq. (2a). The nonlinearity is treated as a perturbation. A similar problem was considered by Tittonen et al. Tittonen (). This approach is convenient if the Floquet states of the time-dependent Hamiltonian are known.
For the driven harmonic oscillator they have been derived by Husimi and Perelomov Husimi (); Popov () and are given in B.
In the second approach, which we call App II, one considers as unperturbed system the undriven nonlinear oscillator (NLO) and the driving is the perturbation, see Eq. (2b).
As we shall see, the different ways of treating the infinite dimensional Floquet Hamiltonian result in crucial differences when evaluating observables of the Duffing oscillator.

## 3 Perturbation theory for a time-periodic Hamiltonian with time-independent perturbation

The starting point of the perturbative treatment App I is the Floquet equation for the full Floquet Hamiltonian in the extended Hilbert space , see Eq. (91),

 ^H|ϕj,m⟩⟩ = ϵj,m|ϕj,m⟩⟩, (4)

where . Moreover, the Floquet states of the Floquet Hamiltonian satisfying the eigenvalue equation (91) are known, see e.g. Eq. (93a) and Eq. (96):

 ^H0|ϕj,m⟩⟩0 = ϵ(0)j,m|ϕj,m⟩⟩0. (5)

We look for an expression of and in first order in . Hence we introduce the first order corrections and as:

 ϵj,m = ϵ(0)j,m+ϵ(1)j,m, (6) |ϕj,m⟩⟩ = |ϕj,m⟩⟩0+|ϕj,m⟩⟩1.

Because the perturbation is time-independent, it is diagonal in the Hilbert space . Additionallly, we introduce the Fourier coefficients:

 0⟨⟨ϕk,n|^Vα|ϕj,m⟩⟩0 ≡ v(n−m)kj. (7)

As in the case of conventional stationary perturbation theory, the perturbed states are written as a linear combination of the unperturbed states:

 |ϕj,m⟩⟩ = |ϕj,m⟩⟩0+∑(i,n)≠(j,m)cnmij|ϕi,n⟩⟩0, (8)

where denotes the couple of quantum numbers and . Inserting ansatz Eq. (6) in the Floquet equation (4) we obtain:

 (^H0+^Vα−ϵ(0)j,m−ϵ(1)j,m)(|ϕj,m⟩⟩0+|ϕj,m⟩⟩1) = 0 (9)

Because and solve the Floquet equation for the last equation reduces to:

 (^Vα−ϵ(1)j,m)|ϕj,m⟩⟩0+(^H0−ϵ(0)j,m)|ϕj,m⟩⟩1 = 0 (10)

This equation allows to determine the modification to the quasienergy and the actual form of the coefficients . To calculate we multiply it from the left with . It follows:

 ϵ(1)j,m = 0⟨⟨ϕj,m|^Vα|ϕj,m⟩⟩0. (11)

To determine the coefficients for the states we multiply Eq. (10) from the left with , where the couple . Moreover we exclude the case of degenerate quasienergies and impose: . In case of degenerate perturbation theory should be applied. We obtain:

 0 = 0⟨⟨ϕk,n|^Vα|ϕj,m⟩0+(ϵ(0)k,n−ϵ(0)j,m)cnmkj, (12)

yielding:

 cnmkj = 0⟨⟨ϕk,n|^Vα|ϕj,m⟩⟩0ϵ(0)j,m−ϵ(0)k,n≡c(n−m)kj. (13)

If we set the driving to zero the quasienergy and the states reduce to the ones obtained by applying conventional stationary perturbation theory on the unforced system.

### 3.1 Application to the quantum Duffing oscillator

We can now determine the actual form of the quasienergy spectrum and the corresponding expansion coefficients for the case of the quantum Duffing oscillator using as perturbation . The matrix elements Eq. (7) defined as:

 v(n−m)kj = 1Tωex∫Tωex0dtexp(i(n−m)ωext)0⟨ϕk(t)|^Vα|ϕj(t)⟩0, (14)

are given in C. The quasienergies of the quantum Duffing oscillator, exact in all orders of the driving strength and up to first order in the nonlinearity, read:

 ϵj,m = ℏΩ(j+12)+F24M(ω2ex−Ω2)+α4⎡⎣32(2j+1)y20(FM(ω2ex−Ω2))2 (15) +32(j(j+1)+12)y40+38(FM(ω2ex−Ω2))4⎤⎦−ℏωexm +O(α2).

In the limit of no driving Eqs. (8) and (15) yield:

 Ej = limF→0ϵj=ℏΩ(j+12)+38αy40(j(j+1)+12), (16a) |j⟩ = limF→0|ϕj(t)⟩=limF→0(t|ϕj⟩⟩ = |j⟩0+αy404⎡⎣√j(j−1)(j−12)2ℏΩ|j−2⟩0 +√(j+1)(j+2)(j+32)−2ℏΩ|j+2⟩0 +14√j(j−1)(j−2)(j−3)4ℏΩ|j−4⟩0+ 14√(j+1)(j+2)(j+3)(j+4)−4ℏΩ|j+4⟩0⎤⎦,

such that the modifications due to the nonlinearity are exactly those obtained by conventional stationary perturbation theory Carmen (), where are the eigenstates of the undriven harmonic oscillator. Expanding up to second order in the driving amplitude we obtain from Eq. (8) for the result:

 |ϕj(t)⟩ = |ϕj(t)⟩0+α4⎡⎣+[y40(j+32)+32y20A2ξ]√(j+1)(j+2)−2ℏΩ|ϕj+2(t)⟩0 +[y40(j−12)+32y20A2ξ]√j(j−1)2ℏΩ|ϕj−2(t)⟩0 +y404√(j+1)(j+2)(j+3)(j+4)−4ℏΩ|ϕj+4(t)⟩0 +y404√j(j−1)(j−2)(j−3)4ℏΩ|ϕj−4(t)⟩0 +34(2j+1)y20A2ξ[exp(−i2ωext)ℏ2ωex+exp(i2ωext)−2ℏωex]|ϕj(t)⟩0 +3!√24(j+1)√j+1Aξy30[exp(−iωext)ℏωex−ℏΩ−exp(iωext)ℏωex+ℏΩ]|ϕj+1(t)⟩0
 +3!√24j√jAξy30[exp(−iωext)ℏωex+ℏΩ+exp(iωext)−ℏωex+ℏΩ]|ϕj−1(t)⟩0 +√(j+3)(j+2)(j+1)23/24y30Aξ[exp(−iωext)ℏωex−3ℏΩ−exp(+iωext)ℏωex+3ℏΩ]|ϕj+3(t)⟩0 +√j(j−1)(j−2)23/24y30Aξ[exp(−iωext)ℏωex+3ℏΩ+exp(+iωext)−ℏωex+3ℏΩ]|ϕj−3(t)⟩0 +34y20A2ξ√(j+1)(j+2)[exp(−i2ωext)ℏ2ωex−2ℏΩ+exp(i2ωext)−ℏ2ωex−2ℏΩ]|ϕj+2(t)⟩0 +34y20A2ξ√j(j−1)[exp(−i2ωext)ℏ2ωex+2ℏΩ+exp(i2ωext)−ℏ2ωex+2ℏΩ]|ϕj−2(t)⟩0],

where we used the abbreviation and are the Floquet states of the linear oscillator.

## 4 Perturbative approach for the one-photon resonance

When the nonlinearity becomes a relevant perturbation to the equidistant spectrum of the linear oscillator, it becomes preferable to use the second approximation scheme, App II, based on the decomposition Eq. (2b).
In this case it is convenient to express the Floquet Hamiltonian in the composite Hilbert space spanned by the vectors , where is an eigenstate of the nonlinear oscillator given in Eq. (16). Hence, in this basis the Floquet Hamiltonian of the nonlinear oscillator , see Eq. (19) below at vanishing driving amplitude, is diagonal. In contrast, the perturbation is time-dependent and thus non-diagonal also in the Hilbert space . From the relation:

 ⟨⟨j,n|^HDO|k,n′⟩⟩ = (^HDO)(n−n′)jk−ℏωexnδjkδnn′, (18)

where are the Fourier coefficients of the matrix , it follows

 ⟨⟨j,n|^HDO|k,n′⟩⟩ = Ej,nδkjδnn′+F2⟨j|^y|k⟩(δn,n′+1+δn,n′−1), (19)

with and the energies of the nonlinear oscillator Eq. (16a). From Eq. (19) is is thus apparent that two eigenstates , of become degenerate when , i.e. for a driving frequency satisfying

 ℏωex(n′−n) = Ek−Ej. (20)

Setting one speaks of a -photon resonance. From Eq. (16a) for the energies it follows, with ,

 ℏωexN = Ej+N−Ej=N[ℏΩ+38αy40(N+1+2j)]. (21)

In the following we restrict to the one-photon resonance , i.e., the quasienergies and are degenerate if .
Moreover, due to the arbitrariness in the choice of the Brillouin zone index , we fix it in the following to the zeroth Brillouin zone, i.e., . For our perturbative treatment we further require that the nonlinearity is large enough that if is resonant with the remaining quasi-energy levels are off resonance and not involved in the doublet spanned by the two degenerate levels. Having this in mind, we have to restrict ourselves to a certain range of possible values of , namely to the resonance region such that the chosen doublet remains degenerate or almost degenerate, i.e. for the one-photon resonance: . This results from the fact that if , the closest lying levels and are by away. Because of the manifold (doublet) structure of the quasi-energy spectrum, we apply in the following Van Vleck perturbation theory Shavitt1980 (); Cohen1992 () and treat the driving as a small perturbation, i.e. . Consequently, a consistent treatment in App II requires that either contributions are neglected if we consider the nonlinearity only up to first order, or that both driving and nonlinearity are treated up to second order. As the second order in both parameters is very involved, we restrict to the first order in the nonlinearity and neglect quadratic contributions in the driving strength, as long as their reliability cannot be verified within a different approach, i.e. App I.
Within Van Vleck perturbation theory we construct an effective Floquet Hamiltonian having the same eigenvalues as the original Hamiltonian and not containing matrix elements connecting states belonging to different manifolds. Therefore it is block-diagonal with all quasi-degenerate energy states in one common block. To determine the transformation and the effective Hamiltonian we write both as a power series in the driving:

 ^S = ^S(0)+^S(1)+^S(2)+… (22) ^Heff = ^H(0)eff+^H(1)eff+^H(2)eff+…. (23)

In D the general formulas for the energies and the states up to second order are provided Shavitt1980 (); Cohen1992 (); Certain (); Kirtman ().
The zeroth order energies are and and the corresponding (quasi)-degenerate Floquet states are: and .
The quasi-degenerate block of the effective Hamiltonian in this basis up to second order in the driving strength acquires the form:

 ^Heff = ⎛⎝Ej,0+E(2)−−jE(1)jE(1)jEj+1,+1+E(2)++j⎞⎠, (24)

where

 E(1)j=⟨⟨j,0|^VF|j+1,+1⟩⟩=⟨⟨j+1,+1|VF|j,0⟩⟩=n1(j)y0F2√2, (25)

and

 E(2)−−j = y20F28[n21(j−1)Ej,0−Ej−1,−1+n21(j−1)Ej,0−Ej−1,+1+n21(j)Ej,0−Ej+1,−1], (26) E(2)++j = y20F28[n21(j+1)Ej+1,+1−Ej+2,+2+n21(j+1)Ej+1,+1−Ej+2,0+n21(j)Ej+1,+1−Ej,2],

with

 n1(j) = √j+1[1−38ℏΩα(j+1)]. (27)

Notice that the unperturbed quasienergies and are correct up to first order in the nonlinearity. For consistency also has to be treated up to first order in only.
As shown by Eq. (4) below it is essential to determine the eigenenergies of up to second order in . They are also the eigenenergies of and read:

 ϵ∓j = 12(Ej,0+Ej+1,+1+E(2)−−j+E(2)++j)±12[(Ej,0−Ej+1,+1)2 +2(Ej,0−Ej+1,+1)(E(2)−−j−E(2)++j)+4E(1)2j]1/2.

The convention is chosen such that for , whereas it jumps at resonance such that for . Because the first order correction in the driving enters Eq. (4) quadratically, a calculation of the quasienergies up to first order in merely yields (when ) the zeroth order results. Consequently to be consistent one has to take into account also the second order corrections and to the energies.
The eigenstates of the block Eq. (24) are determined by:

 |−j,0⟩⟩ eff := −sinηj2|j+1,+1⟩⟩+cosηj2|j,0⟩⟩, (29) |+j,1⟩⟩ eff := sinηj2|j,0⟩⟩+cosηj2|j+1,+1⟩⟩,

where

 tanηj = 2|E(1)j|−(Ej,0−Ej+1,+1+E(2)−−j−E(2)++j). (30)

In conventional Van Vleck perturbation theory the eigenstates of are obtained by applying a back transformation:

 |∓j,n⟩⟩ = exp(−i^S)|∓j,n⟩⟩eff. (31)

Expanding the exponential up to first order we obtain for the eigenstates:

 |∓j,n⟩⟩ = |∓j,n⟩⟩eff−i^S(1)|∓j,n⟩⟩eff, (32) = |∓j,n⟩⟩eff+^R^VF|∓j,n⟩⟩eff.

For reasons given in section 5 and the chosen parameter regime of App II we do not determine the second order correction for the states coming from the second order contribution to .
In the second line of Eq. (32) we used the fact that we can express the transformation by introducing the reduced resolvent , allowing a nice connection to conventional degenerate perturbation theory as shown in D. From Eq. (32) it follows:

 |−j,0⟩⟩ = −sinηj2(1+^R^VF)|j+1,+1⟩⟩+cosηj2(1+^R^VF)|j,0⟩⟩, (33) |+j,1⟩⟩ = sinηj2(1+^R^VF)|j,0⟩⟩+cosηj2(1+^R^VF)|j+1,+1⟩⟩,

where

 ^R^VF|j,0⟩⟩ = ∑(k,n)≠{(j,0),(j+1,+1)}|k,n⟩⟩⟨⟨k,n|^VF|j,0⟩⟩Ej,0−Ek,n = y0F2√2(n1(j−1)Ej,0−Ej−1,−1|j−1,−1⟩⟩+n1(j)Ej,0−Ej+1,−1|j+1,−1⟩⟩ +n1(j−1)Ej,0−Ej−1,+1|j−1,+1⟩⟩+n3(j,α)Ej,0−Ej+3,+1|j+3,+1⟩⟩ +n3(j−3,α)Ej,0−Ej−3,+1|j−3,+1⟩⟩+n3(j,α)Ej,0−Ej+3,−1|j+3,−1⟩⟩ +n3(j−3,α)Ej,0−Ej−3,−1|j−3,−1⟩⟩),
 ^R^VF|j+1,+1⟩⟩ = ∑(k,n)≠{(j,0),(j+1,+1)}|k,n⟩⟩⟨⟨k,n|^VF|j+1,+1⟩⟩Ej+1,+1−Ek,n = y0F2√2(n1(j)Ej+1,1−Ej,2|j,+2⟩⟩+n1(j+1)Ej+1,1−Ej+2,2|j+2,+2⟩⟩ +n1(j+1)Ej+1,1−Ej+2,0|j+2,0⟩⟩+n3(j−2,α)Ej+1,1−Ej−2,0|j−2,0⟩⟩ +n3(j−2,α)Ej+1,1−Ej−2,2|j−2,+2⟩⟩+n3(j+1,α)Ej+1,1−Ej+4,0|j+4,0⟩⟩ +n3(j+1,α)Ej+1,+1−Ej+4,2|j+4,+2⟩⟩),

and

 n3(j,α) = α16ℏΩ√(j+3)(j+2)(j+1).

The effect of the transformation is to yield a contribution from states outside the manifold. Notice that in order to obtain the states to first order in the trigonometric functions and should be expanded in powers of .
We conclude this section by mentioning that eigenenergies and eigenstates of the Duffing oscillator have been calculated near and at resonance also by Peano et al. Peano2 (). However in Peano2 () the nonlinear undriven Hamiltonian is approximated by , where is the occupation number operator of the undriven linear oscillator. This approximated Hamiltonian is diagonal in the linear oscillator basis and yields the result Eq. (16a) for the energies of . However, further corrections of order contained in the eigenstates (16) are neglected. The results of Peano2 () at finite driving can be retained from Eqs. (4) and (33) by treating the driving up to first order, by replacing by and by setting .

## 5 Comparison of the outcomes of the two approaches

The approximation scheme in Sec. 3, App I, is valid when the quasienergy spectrum of the linear oscillator is non-degenerate, i.e., away of a -photon resonance. In contrast, the perturbative approach of Sec. 4, denoted as App II, works at best near a -photon resonance in the quasienergy spectrum of the undriven nonlinear oscillator. Thus a comparison of the outcomes of the two approaches is possible in the frequency regime near resonance, i.e. within . Additionallly, as the Van Vleck-based approach is perturbative in the driving, remember , a comparison requires an expansion in of the results from App I.
This section is organized as follows: First the energies and then the matrix elements of the position operator are compared.

### 5.1 Comparison of the quasienergies

We start with the off resonant case and expand the result in Eq. (4) up to second order in the driving amplitude :

 ϵ−j = Ej,0+E(2)−−j+E(1)2j/(Ej,0−Ej+1,+1), (35) ϵ+j = Ej+1,+1+E(2)++j−E(1)2j/(Ej,0−Ej+1,+1).

Expanding further for consistency the eigenvalues up to first order in the nonlinearity we obtain:

 ϵ−j = Ej,0+y20F28[2Ωℏ(ω2ex−Ω2)+3αy40ℏ2(2j+1)Ω2(ω2ex−Ω2)2] +O(F3,α2) ϵ+j = Ej+1,+1+y20F28[2Ωℏ(ω2ex−Ω2)+3αy40ℏ2(2j+3)Ω2(ω2ex−Ω2)2]+O(F3,α2).

Inserting , these are exactly the results obtained from App I for <