PHYSICAL REVIEW A 89, 063805 (2014)

Two-mode squeezed states in cavity optomechanics via engineering of a single reservoir
M. J. Woolley1 and A. A. Clerk2
1

School of Engineering and Information Technology, UNSW Canberra, ACT, 2600, Australia
2
Department of Physics, McGill University, Montréal, QC, Canada H3A 2T8(Received 9 April 2014; published 11 June 2014)

We study theoretically a three-mode optomechanical system where two mechanical oscillators are independently coupled to a single cavity mode. By optimized two-tone or four-tone driving of the cavity, one canprepare the mechanical oscillators in an entangled two-mode squeezed state, even if they start in a thermal state.
The highly pure, symmetric steady state achieved allows the optimal fidelity of standard continuous-variable
teleportation protocols to be achieved. In contrast to other reservoir engineering approaches to generating
mechanical entanglement, only a single reservoir is required to prepare the highly pure entangled steady state,
greatly simplifying experimental implementation. The entanglement may be verified via a bound on the Duaninequality obtained from the cavity output spectrum. A similar technique may be used for the preparation of ahighly pure two-mode squeezed state of two cavity modes, coupled to a common mechanical oscillator.
DOI: 10.1103/PhysRevA.89.063805

PACS number(s): 42.50.Dv, 03.67.Bg, 42.50.Lc, 85.85.+j

I. INTRODUCTION

The generation and detection of entangled states of macroscopic mechanical oscillators is an outstanding task in the
study of mechanical systems in the quantum regime [1].
There exist a number of proposals for the generation ofsuch states [2–6]. Perhaps most promising amongst these areapproaches based on reservoir engineering [7–10], whereby
the dissipation is engineered such that the steady state ofthe dissipative dynamics is the desired target state. These
proposals are highly attractive from an experimental pointof view, requiring relatively minor modifications of existing
experimental configurations [11–14].
Here, we propose an approach for generating highly pure,highly entangled two-mode squeezed states of two mechanicaloscillators via coupling to a driven cavity mode. The two-mode
squeezed state is the simultaneous vacuum of two nonlocal
bosonic operators (so-called Bogoliubov modes) [15]. Hence,
it is possible to prepare a mechanical two-mode squeezed state
by cooling these two modes. This can be achieved using two
independent reservoirs (see, e.g. Refs. [7–9]). Here, however,
we show that the same goal can be achieved using just a single
reservoir. By making the Bogoliubov modes nondegenerate,they will couple to different frequency components of a single
reservoir (that is, the damped cavity). Since the Bogoliubovtransformation preserves the difference in number operators,
this simply corresponds to a frequency difference of the two
mechanical oscillators. Viewed differently, one can say that
we are exploiting the coherent dynamics of the mechanicaloscillators to effectively cool both Bogoliubov modes. Ultimately, our protocol simply involves appropriately weighted
and detuned two-tone or four-tone driving of the coupled cavity
mode.
The approach we take here may be regarded as the coherent
feedback [16] analog of our measurement-based approach tothe same task [5], and is related to a recent proposal for the
preparation of a quantum squeezed state of a single mechanicaloscillator [17]. As compared with the measurement-based
approach for entanglement generation in Ref. [5], the purity ofthe steady state achieved here is greater and the implementation
is greatly simplified. Further, in contrast to the proposal
1050-2947/2014/89(6)/063805(17)of Ref. [8], the steady state is a (highly pure) two-mode
squeezed state rather than an entangled mixed state. Therefore,
using the steady state as an Einstein-Podolsky-Rosen (EPR)
channel, the optimal teleportation fidelity for a given amountof entanglement can be achieved via the standard continuousvariable teleportation protocol, without the need for additional
local operations [18,19]. In contrast to the proposal of Ref. [9](also described in the supplement to Ref. [8]), the two-mode
squeezed state is generated here using only one, rather thantwo, auxiliary cavity modes, simplifying the experimental
implementation.
Reservoir engineering has earlier been the subject ofsignificant theoretical study in the context of optical andatomic systems [7,20–25], culminating in the experimental
demonstration by Krauter and co-workers of the entanglement of atomic ensembles [26]. The utility of reservoir
engineering has also been shown with two-level systems,with demonstrations of superconducting qubit state control[27], as well as entanglement in both trapped ion [28] andsuperconducting [29] systems. In other work pertaining tomechanical entanglement, the entanglement of mechanicalmotion with a microwave field has been demonstrated [30].
Earlier, entanglement of phonons at the single-quantum level
was demonstrated [31], as well as the entanglement of motional
degrees of freedom of trapped ions [32].
Here, in Sec. II we introduce the multimode optomechanical system that we shall study. Section III describes
approaches that one may take to reservoir engineering inthis system, while in Sec. IV we describe how these strategies could be implemented in our system. In Sec. V, we
consider the adiabatic limit, in which the cavity responds
rapidly to the mechanical motion, and obtain analyticalexpressions for the entanglement, purity, and teleportation
fidelity possible with the steady state. Section VI gives ananalysis of the full linearized system, including the effects ofcounter-rotating terms and the possibility of instability. In Sec.
VII, we derive a bound on the entanglement based on the cavity
output spectrum. Section VIII provides an analysis of the threemode optomechanical system composed of two cavity modes
coupled to a single mechanical oscillator, and demonstrates
that the same physics can be realized in this system.

063805-1©2014 American Physical Society

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

precisely in Eqs. (26a) and (26b), are rotating in a nontrivial
manner. In particular,
1 +i(ωa −)t[âe+ â † e−i(ωa −)t2+ b̂e+i(ωb +)t + b̂† e−i(ωb +)t ],iP̂− = − [âe+i(ωa −)t − â † e−i(ωa −)t2+ b̂e+i(ωb +)t − b̂† e−i(ωb +)t ].

X̂+ =FIG. 1. (Color online) (a) The system studied consists of two
mechanical oscillators, each coupled to a common cavity (or circuit)mode. (b) Frequencies in this system defined with respect tothe cavity resonance frequency ωc . The blue lines indicate thestandard mechanical sidebands, at ±ωa and ±ωb . If the single-photon
optomechanical coupling rates (ga and gb ) are equal, the required
Hamiltonian (8) can be realized using only two cavity drive detunings(±ωm ), indicated by the vertical red lines.
II. SYSTEM AND HAMILTONIAN

The system [see Fig. 1(a)] is composed of two mechanicaloscillators, with resonance frequencies ωa and ωb , eachindependently, dispersively coupled (with strengths ga andgb , respectively) to a common cavity mode having resonancefrequency ωc . The Hamiltonian isĤ = ωa â † â + ωb b̂† b̂ + ωc ĉ† ĉ + ga (â + â † )ĉ† ĉ
+ gb (b̂ + b̂† )ĉ† ĉ + Ĥdrive + Ĥdiss ,(1)

where â and b̂ denote mechanical mode lowering operators,ĉ denotes the electromagnetic mode lowering operator, andĤdrive accounts for driving of the electromagnetic mode. The
term Ĥdiss accounts for dissipation, with the modes subject
to damping at rates γa , γb , and κ, respectively. The system
dynamics, within the usual approximations [33], is describedby the master equationρ̇ = −i[Ĥ ,ρ] + γa (n̄a + 1)D[â]ρ + γa n̄a D[â † ]ρ+ γb (n̄b + 1)D[b̂]ρ + γb n̄b D[b̂† ]ρ + κD[ĉ]ρ, (2)with the Hamiltonian Ĥ = Ĥ − Ĥdiss and the dissipative
superoperator D[ŝ]ρ = ŝρ ŝ † − 12 ŝ † ŝρ − 12 ρ ŝ † ŝ.
III. RESERVOIR ENGINEERING STRATEGIESGiven the two mechanical oscillator modes â and b̂, one
can introduce mechanical two-mode Bogoliubov operators insome rotating frame,β̂1 = â cosh r + b̂† sinh r,†β̂2 = b̂ cosh r + â sinh r.

(3a)
(3b)

Typically, the choice of rotating frame is set by the resonancefrequencies of the system, although here the rotating frame
shall be defined with respect toĤ0 = (ωa − )â † â + (ωb + )b̂† b̂ + ωc ĉ† ĉ,(4)with the choice of  to be specified in Sec. IV. Note that inthis frame the collective mechanical quadratures, to be defined

(5a)

(5b)

In any case, the two-mode squeezed state is defined as|r2 = Ŝ2 (r)|0,0, whereŜ2 (r) ≡ exp[r(â b̂ − â † b̂† )](6)is the two-mode squeezing operator with squeezing parameter r [34]. Starting from â |0,0 = 0 and b̂ |0,0 = 0 it†
is straightforward to show that [Ŝ2 (r)â Ŝ2 (r)] |r2 = 0 and††that [Ŝ2 (r)b̂Ŝ2 (r)] |r2 = 0. However, β̂1 = Ŝ2 (r)â Ŝ2 (r) and†β̂2 = Ŝ2 (r)b̂Ŝ2 (r), and therefore, the ground state of β̂1 andβ̂2 is the two-mode squeezed state with squeezing parameterr. Our goal then is to engineer the driving Hamiltonian of (1)such that the steady state of (2) results in the β̂i modes being
cooled to their ground state, implying two-mode squeezing of
the mechanical oscillators.
One method to achieve this is to use two cavity modes
to independently cool the Bogoliubov modes [7–9]. Both
Bogoliubov modes are independently coupled to a cavity†
mode with a beam-splitter-like interaction, i.e., β̂i ĉi + H.c.
(i = 1,2), an interaction which can be used to cool the
Bogoliubov modes [35]. While such an approach can beeffective, from a practical point of view it would be highly
advantageous if this could be achieved using only a single
cavity mode.
A seemingly simple way of using only a single reservoir
would be to couple the cavity to one of the Bogoliubov modes,
say β̂1 , and then couple β̂1 to the other Bogoliubov mode(β̂2 ) via an all-mechanical contribution to the Hamiltonian of†
the form β̂1 β̂2 + H.c., again a beam-splitter-like interaction.
This interaction will allow β̂2 to be cooled (by swapping
quanta into β̂1 ) even though it is not directly coupled to thecooling reservoir. While conceptually simple, the requisite
mechanical interaction would be difficult to realize, requiring
the direct coupling of the mechanical oscillators. Such a
Hamiltonian was considered in Ref. [36], albeit without aphysical implementation specified.
A third approach, which we will pursue here, is to couple
the cavity
to the sum of Bogoliubov modes β̂sum ≡ (β̂1 +√β̂2 )/ 2, and then couple the sum mode to the differencemode via an all-mechanical Hamiltonian contributionof the√†form β̂sum β̂diff + H.c., where β̂diff ≡ (β̂1 − β̂2 )/ 2. Again,
the swap interaction allows β̂diff to be cooled even though
it is not directly coupled to the cavity. Cooling both β̂sum
and β̂diff is equivalent to cooling both β̂1 and β̂2 since††††β̂sum β̂sum  + β̂diff β̂diff  = β̂1 β̂1  + β̂2 β̂2 . While this approach again seems to involve the realization of a challenging
interaction between Bogoliubov modes, this is not the case.
The beam-splitter interaction here takes the simple form††β̂1 β̂1 − β̂2 β̂2 = â † â − b̂† b̂. Thus, one does not require a direct

063805-2

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)

interaction between the mechanical oscillators, but rather justa difference in their resonance frequencies. This is the key
insight that allows one to appropriately engineer the reservoir
via a single cavity mode.
As noted in Sec. I, we may take another perspective on
this third approach. By introducing a frequency difference
between the two mechanical oscillators, we are breaking
the degeneracy of the Bogoliubov modes 1 and 2 since
the Bogoliubov transformation preserves the number-operator
difference. Consequently, the Bogoliubov modes couple todifferent frequency components of the reservoir. Due to thefinite bandwidth of the cavity, it effectively functions as two
independent reservoirs such that both Bogoliubov modes arecooled.
We thus have that the desired Hamiltonian, in terms of the
Bogoliubov modes defined in Eqs. (3a) and (3b), is††††Ĥ = (β̂1 β̂1 − β̂2 β̂2 ) + G[(β̂1 + β̂2 )ĉ + H.c.] + Ĥdiss ,ωm = (ωa + ωb )/2 is the average of the two mechanicalfrequencies, i.e.,Ĥdrive = (E+∗ e+iωm t + E−∗ e−iωm t )e+iωc t ĉ + H.c.

This situation is depicted in Fig. 1(b). The driving tones
must be applied with a fixed relative phase. Working in aninteraction picture defined with respect to the Ĥ0 = ωm (â † â +b̂† b̂) + ωc ĉ† ĉ, one finds the effective Hamiltonian to be given
by Eq. (8) where = (ωa − ωb )/2,

Ĥ = (â † â − b̂† b̂) + G+ [(â + b̂)ĉ + H.c.]+ G− [(â + b̂)ĉ† + H.c.] + Ĥdiss .

(8)

The optomechanical couplings in Eqs. (7) and (8) are related
by

(9a)
G ≡ G2− − G2+ ,tanh r ≡ G+ /G− ,(9b)with r being the squeezing parameter entering in thedefinitions of the Bogoliubov modes in Eqs. (3a) and (3b).
Note that if G+ = G− in Eq. (8), then we recover the
Hamiltonian required for a two-mode backaction-evading
measurement of the mechanical oscillators [5], in which twocollective mechanical quadratures commute with the system
Hamiltonian. For G+ = G− , the backaction evasion is lost, but
now there is a backaction that may be regarded as a coherent
feedback process. This enables two-mode squeezing without
an explicit measurement.

(11)

and the many-photon optomechanical couplings are
G± = (ga + gb )c̄± /2,(12)with c̄± denoting the (assumed real) steady-state amplitudesof the fields at the driven sidebandsc̄± ≡ ĉ± ss =(7)
where  is an effective oscillation frequency and G is aneffective optomechanical coupling. In terms of the original
mechanical annihilation operators, Eq. (7) is

(10)

iE±
.
±iωm − κ/2(13)

The details of this derivation are given in Appendix A. It
relies on the assumptions that we are working in the resolvedsideband regime (ωa ,ωbκ) and that the driving strengths
E± are large. The former assumption allows us to discard
time-dependent contributions to the Hamiltonian (8), while
the latter assumption allows us to linearize the optomechanicalinteraction. Note that the ratio G+ /G− shall be referred to here
as the drive asymmetry since it is set by the ratio of the cavity
drives on either side of the cavity resonance frequency.
If the single-photon optomechanical couplings are unequal,
the two-tone cavity driving can not yield the complete
matching of oscillator a and b sideband processes requiredin the ideal Hamiltonian of Eq. (8). Instead, there will be
additional contributions to the Hamiltonian (8), given bym††Ĥm = Gm
+ [(â − b̂)ĉ + H.c.] − G− [(â − b̂ )ĉ + H.c.], (14)where the effective coupling imperfections areGm± = ±(ga − gb )c̄± /2.

(15)

In the two-tone driving case, the imperfection is due to themismatch in the single-photon optomechanical coupling rates.
B. Four-tone driving

IV. IMPLEMENTATIONThe Hamiltonian (8) is readily implemented in conventional
cavity optomechanics setups. We shall focus on the regime|G+ | < |G− | such that the dynamics corresponding to (8)are stable. If the single-photon optomechanical coupling ratesin (1) are equal (ga = gb ), then just two cavity drives are
required to realize (8). If they are unequal (ga = gb ), thenfour cavity drives are required. Of course, if ga ∼ gb , we canapproximately realize (8) with only two cavity drives and still
generate useful entanglement in the steady state. We consider
each of these cases in turn.

The Hamiltonian (8) involves four sideband processes; theup conversion and down conversion of drive photons via theabsorption (or emission) of quanta from (or to) the mechanicaloscillator a or b. The realization of (8) requires a balance of therates at which these processes take place. By using four cavity
driving tones, one tone associated with each sideband process,the balancing of the rates of these processes is possible even
if the single-photon optomechanical couplings are unequal(see Fig. 2). These driving tones are applied with a detuningof  from the mechanical sidebands, at ωc ± (ωa − ) andωc ± (ωb + ), as depicted in Fig. 2. The appropriate Hamiltonian contribution is

A. Two-tone driving

If the single-photon optomechanical coupling rates areequal, we require cavity driving tones at ωc ± ωm , where
063805-3∗ +i(ωa −)t
∗ +i(ωb +)tĤdrive = e+iωc t ĉ(E1+
e
+ E2+
e
∗ −i(ωa −)t
∗ −i(ωb +)t
+ E1−
e
+ E2−
e
) + H.c. (16)

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

condition (22b) ensures that the unwanted sideband processes
have a negligible effect on the system dynamics. It issometimes convenient to refer to the drive frequencies via their
detunings from the center of the two mechanical sidebands,
given byδ ≡ (ωa − ωb )/2 − .

FIG. 2. (Color online) (a) The three-mode optomechanical system, as in Fig. 1, under four-tone driving. (b) Frequencies in thissystem defined with respect to the cavity resonance frequency ωc .
The blue lines indicate the standard mechanical sidebands, at ±ωa and±ωb . If the single-photon optomechanical coupling rates are unequal,
the required Hamiltonian (8) can be realized using four cavity driving
frequencies ±(ωa − ), ± (ωb + ), indicated by vertical red lines.

The steady-state amplitudes at the driven sidebands aredenoted by c̄k± (k = 1,2), withc̄k± ≡ ĉk± ss =iEk±,±iωk − κ/2(17)where we have introduced the notation for the drive detunings
ω1 ≡ ωa − ,

(18a)

ω2 ≡ ωb + .

(18b)Then, we demand that the driving strengths are “matched,”
meaning that
gb
c̄1±
= .
(19)
c̄2±
ga
That is, we require that the two steady-state amplitudes (i.e.,
drives) on the same side of the cavity resonance frequency havean asymmetry set by the optomechanical coupling asymmetry.
With the condition (19) satisfied, and working in an interaction
picture defined with respect to the Hamiltonian (4), the system
can again be described by the Hamiltonian (8), now with the(assumed real) many-photon optomechanical coupling rates
G± = (ga c̄1± + gb c̄2± ) /2.

Gm± = ± (ga c̄1± − gb c̄2± ) /2.

It is interesting to note that with the driving conditionc̄1+ = c̄1− , in addition to (19), we could realize a two-mode
backaction-evading measurement of the mechanical oscillators
irrespective of the coupling asymmetry, and so generalize theresults of Ref. [5].
V. ADIABATIC LIMIT

We first consider the dynamics of the system governed by(8) in the adiabatic limit, where the cavity responds rapidly
to the mechanical motion; that is, where κ > ,G± (but stillin the regime where ωa ,ωbκ). In this limit, we eliminate
the cavity mode, obtaining an effective description for themechanical modes alone. This adiabatic limit will simplify the
analysis and thus provide insight into our mechanism; it will
also prove to be a useful regime for the task of mechanicalentanglement generation.
We stress that the Hamiltonian (8) applies both to thecase of equal single-photon optomechanical couplings andtwo-tone driving, and to the case of unequal single-photon
optomechanical couplings with matched four-tone driving(although  is determined differently in each case). Further,
the imperfections (asymmetry in couplings in the former caseand mismatch in driving conditions in the latter case) are both
described by the Hamiltonian (14).
In this adiabatic limit, the cavity annihilation operator
is given by ĉ = −2iG(β̂1 + β̂2 )/κ. Substituting this into thedissipative terms of the master equation (2), the adiabaticallyeliminated master equation is††ρ̇ = −i[β̂1 β̂1 − β̂2 β̂2 ,ρ] + γa (n̄a + 1)D[â]ρ + γa n̄a D[â † ]ρ+ γb (n̄b + 1) D[b̂]ρ + γb n̄b D[b̂† ]ρ + D[β̂1 + β̂2 ]ρ,
(24)

(20)

Again, the details of the derivation are left to Appendix A.
Imprecision in the matching condition (19) gives additional
contributions to the Hamiltonian (8), of the form of Eq. (14),but now with

with the optomechanical damping rate
≡ γC ≡(21)

In this case, the effective coupling imperfection arises fromthe drives not being weighted precisely according to thecondition (19).
The cavity drive frequencies should be set such that satisfies the following conditions:


γ,

(22a)



(ωa − ωb )/2 − γ .

(22b)Condition (22a) ensures that the sum and difference Bogoliubov modes are sufficiently coupled [cf. Eq. (7)]. The

(23)

4G 2
,
κ

(25)

where G is the effective optomechanical coupling introduced in
Eq. (9a), and C is the corresponding cooperativity parameter.
Now, the steady state of Eq. (24) is easily obtained, and its
entanglement and purity metrics readily calculated.
In view of Eq. (24), an alternative interpretation of thecooling of both Bogoliubov modes is possible. In the limit††γ /  → 0, the terms containing β̂1 β̂2 and β̂2 β̂1 will rapidly
average away, meaning that Eq. (24) will be equivalent tohaving independent dissipation of modes β̂1 and β̂2 . Physically,
this corresponds to the dissipation of β̂1 and β̂2 being due todistinct modes of the reservoir.

063805-4

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)

A. Entanglement

The case of symmetric mechanical damping (γa ,γb = γ )and symmetric thermal occupation of the mechanical baths(n̄a ,n̄b = n̄) allows simple analytical results for the steadystate second moments to be obtained. The assumption ofequal thermal occupations is reasonable for most experimental
situations, while it turns out that our results are not sensitive
to unequal mechanical damping rates provided that they
are both small. The simplest two-mode, continuous-variableentanglement criterion is provided by the Duan inequality [37].
It is expressed in terms of collective quadrature operators,
defined by√
X̂± = (X̂a ± X̂b )/ 2,(26a)√
P̂± = (P̂a ± P̂b )/ 2,(26b)where we have introduced the usual quadratures for each
oscillator mode√√(27)
X̂s = (ŝ + ŝ † )/ 2, P̂s = −i(ŝ − ŝ † )/ 2.
Then, the Duan criterion tells us that a Gaussian state for which2 + P̂−2  < 1
X̂+

(28)

is inseparable. Note that this could equally well be formulated
in terms of X̂− and P̂+ , although Eq. (28) shall be the suitableform here.
For our system we consider the limit γ /  → 0 since thisensures that the sum and difference Bogoliubov modes are
sufficiently coupled (or, equivalently, that the two individual
Bogoliubov modes see effectively independent reservoirs). For
G− = G+ , we find for the steady-state second moments2 = P̂∓2  =
X̂±

γ
γ+

(n̄ + 1/2) +

γ+

e∓2r
2

γκ
=
(n̄ + 1/2)γ κ + 4(G2− − G2+ )2(G− ∓ G+ )+.
γ κ + 4(G2− − G2+ )(29a)Fig. 3(a), it is seen that the Duan quantity takes a value below
one for experimentally reasonable parameters, achievable instate-of-the-art microwave cavity optomechanics experiments[12–14], indicating that the mechanical oscillators are entangled in the steady state. This continues to be the case even when
one accounts for large nonzero initial thermal occupations
and large imperfections in the effective couplings. Further,
the logarithmic negativity, shown in Fig. 3(b), takes a large
value for these parameters. For comparison, the logarithmic
negativity of a stationary two-mode squeezed state generated
via a Hamiltonian parametric amplifier interaction is bounded
above, due to a stability constraint, by EN = ln 2 ∼ 0.69.
As previously noted [8,17], the entanglement goes througha maximum as a function of the drive asymmetry. For themechanical steady state to be highly entangled, we require
both that the target steady state is highly squeezed (r large,
requiring G+ /G− → 1) and that the system is effectively
cooled towards this steady state ( large, requiring G+ /G− →0). Obviously, these limits are incompatible and the optimal
asymmetry is between these, leading to the observed maximum.
Given the simple analytical results we have obtained, we
may optimize the steady-state entanglement analytically byminimizing the Duan quantity over the drive asymmetry. Thisis most conveniently done using Eq. (29b). We find that the
optimal drive asymmetry is


G+ 
1 + 1/C−
1 + n̄= 1+−(31a)
G− opt.
C−C−
1∼ 1− √ .
C−

The result (31b) holds in the large-cooperativity limit,
provided that one is still within the adiabatic regime. It follows
that the Duan quantity (28), to first order in C−−1 , is
1 + n̄ n̄(1 + n̄)2 + P̂−2  = √
+
.
X̂+
C−
C−

2

(29b)

Equation (29a) takes a particularly simple form, describing
coupling to a squeezed reservoir with an optomechanicaldamping rate . The results (29a) and (29b) are easily checkedagainst the solution of the full system (i.e., without the
adiabatic elimination), as discussed in Sec. VI.
The G+ = G− limit is unclear from Eq. (29a), as itcorresponds to the limits → 0 and r → +∞. However,
the result is clear from Eq. (29b), and we recover the result2that X̂+ = P̂−2  = n̄ + 1/2 (evading the backaction) and2X̂−
 = P̂+2  = n̄ + 1/2 + C± /2 (heated by the backaction)[5], where the cooperativities associated with the blue and red
sideband drives alone (denoted by the subscripts “+” and “−”,respectively), are
C± ≡

±

γ

≡

4G2±
.
γκ

(30)

While the Duan inequality provides a simple entanglement
criterion, the entanglement may be quantified via the logarithmic negativity, defined in Appendix B 1. Both are shown,as functions of the drive asymmetry G+ /G− , in Fig. 3. In

(31b)

(32)

Clearly, the mechanical oscillators are entangled even for amodest cooperativity. We emphasize these results are only
valid in the adiabatic limit, corresponding to C−  4κ/γ . The
achievable entanglement beyond the adiabatic regime shall be
discussed in Sec. VI.
B. Purity

The purity of the steady state generated is relevant for
both experiments in quantum foundations and in quantum
information processing; its role in determining teleportation
fidelity shall be described in Sec. V C. Now, the fact that thesteady state is highly entangled does not necessarily imply thatthe steady state is highly pure. Indeed, if one cools only one
Bogoliubov mode, then the steady state is highly entangled, butalso highly impure [8]. The purity of the mechanical two-mode
steady state is defined asμ ≡ tr(ρ 2 ),(33)

where ρ is the density matrix of the two mechanical modes. It
can be directly evaluated from knowledge of its symmetricallyordered covariance matrix V. With quadratures as defined in

063805-5

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

FIG. 3. (Color online) Entanglement, expressed via (a) the Duanquantity (28) and (b) the logarithmic negativity, against the drive
asymmetry G+ /G− . We hold G− fixed (fixing the cooperativityC− ) and vary G+ . These results correspond to the time-independent(rotating-wave approximation) Hamiltonian (8), and apply to both thetwo-tone and four-tone driving cases, with the effective couplings G±
given by Eqs. (12) and (20), respectively. These plots are obtainedusing the adiabatic limit results of Sec. V, although they coincide withthe results for the full system obtained in Sec. VI A. The solid black
curve corresponds to a mechanical bath thermal occupation of n̄ = 0,m
and Gm± = 0 where G± are the effective coupling imperfections,introduced in Eqs. (14) and (21) for the two-tone and four-tone
driving cases, respectively. The blue curve (long dashes) corresponds
to n̄ = 25 and Gm± = 0, while the red curve (short dashes) corresponds
to n̄ = 25 and Gm± = 0.5G± . Parameters common to each curve are
C− = 1200, κ = 2π × 1.592 × 105 Hz = 106 s−1 , γ /κ = 4 × 10−5 ,and /κ = 0.1.

Eqs. (27) and the covariance matrix expressed in the orderedbasis (X̂a ,P̂a ,X̂b ,P̂b ), the purity is simply given by√
μ = 1/(4 det V).

(34)

The purity may also be assessed by calculating the thermal
occupations of the two mechanical Bogoliubov modes, definedin Eqs. (3a) and (3b).
Plots of both quantities, against the drive asymmetry
G+ /G− , are shown in Fig. 4. From Fig. 4(a) it is clear thatthe occupations of the two Bogoliubov modes are the sameprovided that the imperfections in the effective couplings,
Eqs. (15) and (21), are zero. Further, the occupations are close
to zero for reasonable parameters, verifying that our scheme
effectively cools both Bogoliubov modes. The purity of thestate is correspondingly high, being close to one for reasonable
experimental parameters [see Fig. 4(b)], and so outperforminga scheme in which only one Bogoliubov mode is cooled [8].
In the absence of effective coupling imperfections, we canobtain simple analytical results characterizing the purity of thesteady state. The occupations of the Bogoliubov modes in the

FIG. 4. (Color online) (a) Steady-state occupations of the mechanical Bogoliubov modes, defined in Eqs. (3a) and (3b), and (b)steady-state purity of the two mechanical modes, defined in Eq. (33),against drive asymmetry. The solid black curve corresponds to amechanical bath thermal occupation of n̄ = 0 and no imperfectionin the effective couplings (Gm± = 0), the blue curve (long dashes)corresponds to n̄ = 25 and Gm± = 0, and the red curves (short dashes)correspond to n̄ = 25 and Gm± = 0.5G± . Only one solid black curveand one dashed blue curve is shown in (a) since the occupations of
the two Bogoliubov modes are the same in each of these cases. The
black curves with long and short dashes correspond to bounds onthese quantities for the scheme in which only one Bogoliubov mode
is cooled [8], assuming a thermal occupation of n̄ = 0: the curve in(a) is a lower bound on the occupation of the uncooled Bogoliubov
mode, while the curve in (b) is an upper bound on the purity of thesteady state in this case. Remaining parameters for each curve are as
given in the caption of Fig. 3.

limit γ /  → 0 (and for G− = G+ ) areγ†[n̄ + (2n̄ + 1) sinh2 r](35a)β̂i β̂i  =γ+ 2
G+ + n̄(G2− + G2+ )γκ=, (35b)γ κ + 4(G2− − G2+ )G2− − G2+
for i = 1 and 2, consistent with the result for single-mode
squeezing in Ref. [17]. Again, it is clear that both Bogoliubov
modes are cooled equally. The purity itself, in the limit γ /  →0, is given byμ=(γ + )2[γ (1 + 2n̄) + ]2 + 4(1 + 2n̄)γsinh2 r.

(36)

The purity at the entanglement maximum, in the largecooperativity limit, is μ ∼ 1/[2(1 + n̄)]. However, a highly
pure (and still highly entangled) steady state can be achievedby choosing a drive asymmetry just below that corresponding
to the optimal entanglement (as this gives a large and hence
more cooling, at the expense of a smaller squeeze parameter r).

063805-6

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)

C. Teleportation fidelity

The two-mode squeezed state generated here may be
regarded as the entangled resource (“EPR channel”) in acontinuous-variable teleportation protocol [18]. If we write
the steady-state two-mode symmetrically ordered covariancematrix, in the ordered basis (X̂a ,P̂a ,X̂b ,P̂b ), in block form as
Va VabV=,(37)
VTab Vbthen the teleportation fidelity for a single-mode Gaussian input
state under the standard protocol is given by [38]
1
F= √
,
(38a)
det(2Vin + N)
N = σz Va σz + σz Vab + VTab σz + Vb ,(38b)where Vin is the covariance matrix of the state to be teleported.
For the teleportation of a coherent state, this is Vin = (1/2)I2 .
We find that the mechanical steady state in our system, inthe limit γ / → 0, is a thermal two-mode squeezed state.
Such a state is defined by†ρ ≡ Ŝ2 (ξ ) ρtha ⊗ ρthb Ŝ2 (ξ ),(39)

where ρtha(b) denotes the density matrix of a thermal state ofmode a(b) with occupation n̄a(b)th , and Ŝ2 (ξ ) is the two-mode
squeezing operator introduced in Eq. (6). Accordingly, wea(b)can assign effective occupations n̄th
and an effective twomode squeezing parameter ξ to our steady state, as detailed inAppendix B 2. Now, the purity of such a state is simply
1
μ=
,
(40)
1 + 2n̄ath 1 + 2n̄bth
and the teleportation fidelity, based on a thermal two-mode
squeezed state channel, is [38]1
F = −2ξ
.
(41)
1 + n̄ath + n̄bth + e2ξ
e
Clearly, larger effective occupations correspond to a lower
purity and a lower teleportation fidelity.
It is known that when the channel of a continuous-variable
teleportation protocol is a symmetric thermal two-mode
squeezed state (n̄ath ,n̄bth ≡ n̄th ), the teleportation fidelity issimply [19]
1
,
(42)
F=
1 + e−EN
where the logarithmic negativity is given by EN =Max[0,2ξ − ln(1 + 2n̄th )]. Equation (42) actually gives the
optimal teleportation fidelity achievable for a given amountof entanglement [19]. Asymmetry results in a teleportation
fidelity below that given by Eq. (42).
The effective thermal occupations for our steady state(n̄ath and n̄bth ), as a function of the drive asymmetry, areshown in Fig. 5(a). It is clear that n̄ath = n̄bth provided that theimperfections in the effective couplings are zero. Accordingly
in these cases, the teleportation fidelity is given by Eq. (42) andso attains its optimal value for a given amount of entanglement.
thConsider now the case of an asymmetric channel n̄tha = n̄b ;this is generically the kind of state produced using a scheme
which cools only a single Bogoliubov mode [8]. For such

FIG. 5. (Color online) (a) Effective thermal occupations for eachmechanical mode in the steady state, introduced in Eq. (39), against
the drive asymmetry. (b) The teleportation fidelity, using the generated
mechanical steady state as an EPR channel to teleport a coherent state
via the standard protocol, against the drive asymmetry. The solid black
curves correspond to a mechanical bath thermal occupation of n̄ = 0
and no imperfection in the effective couplings (Gm± = 0), the blue
curves (long dashes) correspond to n̄ = 25 and Gm± = 0, and the red
curves (short dashes) correspond to n̄ = 25 and Gm± = 0.5G± . The
black curve with long and short dashes on the left is the lower boundon the occupation of the uncoupled mode in a scheme that cools one
Bogoliubov mode [8]. The black curve with long and short dashes on
the right is the teleportation fidelity achievable with such a scheme.
The upper (lower) dotted purple curve is the upper (lower) bound on
the optimal teleportation fidelity achievable for a given amount ofentanglement. The given amount of entanglement corresponds to thatpossessed by a two-mode squeezed state with squeezing parameterr = tanh−1 (G+ /G− ) [cf. Eq. (9b)]. Remaining parameters for eachsolid curve are as given in the caption of Fig. 3.

states, the standard teleportation protocol does not achieve thefidelity in Eq. (42). This fidelity can be reached in principle
in the highly entangled regime, if one goes beyond thestandard protocol by allowing for additional local Gaussian
operations [39].
The teleportation fidelity, assuming that the mechanicalsteady state we have generated is used as an EPR channel
and assuming no other sources of imperfection, is plotted
as a function of drive asymmetry in Fig. 5(b). Crucially, bycooling both Bogoliubov modes, the optimal teleportation
fidelity tends to 1 rather than to 45 in the highly entangled
regime. Further, this upper bound is achievable with reasonable
parameters.
D. Coherent feedback

The reservoir engineering scheme that we have describedhere permits an alternative interpretation in terms of coherent
feedback [16], similar to that provided for the squeezing
scheme in Ref. [17]; we depict the process schematically in
Fig. 6. The Hamiltonian (8) can be rewritten in terms of the063805-7

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

positive and the collective quadrature X̂− will be antisqueezedin the steady state.
VI. FULL SYSTEM
A. Solution with time-independent Hamiltonian

FIG. 6. (Color online) Representation of the reservoir engineering scheme of Figs. 1 and 2 as a form of coherent feedback, with the
feedback (backaction) being applied autonomously via the cavity
mode. The description is in terms of the collective mechanicalquadrature operators introduced in Eqs. (26a) and (26b). (a) Thereis a measurement of X̂+ via the cavity at a rate ∼(G− + G+ ) anda feedback onto X̂+ via the cavity at a rate ∼(G− − G+ ). The
measurement rate being greater than the feedback rate leads tosqueezing. (b) There is a measurement of P̂+ via the cavity at a
rate ∼(G− − G+ ) and a feedback onto P̂+ via the cavity at a rate∼(G− + G+ ). The feedback rate being greater than the measurement
rate leads to amplification.

collective mechanical quadratures of Eqs. (26a) and (26b) as√
Ĥ = (X̂+ X̂− + P̂+ P̂− ) + 2(G− + G+ )X̂+ X̂c√+ 2(G− − G+ )P̂+ P̂c + Ĥdiss .
(43)
This Hamiltonian is a perturbation, via the third term, of a
Hamiltonian that we previously studied in the context of twomode backaction-evading measurement and feedback control[5]. From a coherent feedback point of view, the second term in
Eq. (43) may be regarded as the “measurement” interaction andthe third term in Eq. (43) may be regarded as the “feedback”
backaction, applied autonomously via the cavity mode.
After an adiabatic elimination of the cavity mode, theHeisenberg-Langevin equations corresponding to the Hamiltonian (43) may be written asd 
 + B1 · X
 in + B2 · Yin ,X = A0 · X
dt

(44)

 = (X̂+ ,P̂+ ,X̂− ,P̂− )T is the vector of mechanicalwhere X
 in is the corresponding
collective quadrature observables, X
vector of mechanical input noises, and Yin ≡ [Ŷ1 (t),Ŷ2 (t)]T
are the operators associated with the cavity input noise. The
remaining matrices in Eq. (44) are specified in Appendix C 1.
The new noise input operators have the correlation functions
1 G− − G+1
δ(t),δ(t) ≡ n̄1 +Ŷ1 (t)Ŷ1 (0) =2 G− + G+
2(45a)

1 G− + G+
1
δ(t),δ(t) ≡ n̄2 +2 G− − G+2(45b)

Ŷ2 (t)Ŷ2 (0) =where n̄1 and n̄2 denote effective thermal occupations of
the noise inputs. Clearly, the input noises seen by thecollective mechanical quadratures are weighted by the ratiosof the measurement and feedback rates. Since we have
G− > G+  0 here, the effective occupation n̄1 is negative[17]. Therefore, as far as the collective quadrature X̂+ isconcerned, the cavity behaves as a squeezed bath and thereforeX̄+ will be squeezed in the steady state [33]. Conversely, n̄2 is

In Sec. V, we calculated the steady state of our system in
the adiabatic limit, after mathematically removing the cavity
mode from the system. Even with the cavity mode retained, theeffective Hamiltonian (8) is quadratic and time independent,
and we may readily solve for the steady state. The HeisenbergLangevin equations corresponding to (8) may be written as thesystemd 
 + B0 · X
 in ,X = A0 · X
dt

(46)

 = (X̂a ,P̂a ,X̂b ,P̂b ,X̂c ,P̂c )T is defined in terms of
where X
individual oscillator quadratures, and the matrices are specified
in Appendix C 2. The steady-state, symmetrically ordered
covariance matrix V is obtained by solving the Lyapunov
equation
A0 V + VAT0 = −B0 BT0 .

(47)

Solving Eq. (47) allows us to assess the steady state even when
we are not in the adiabatic limit, and to validate results obtainedin the adiabatic limit. As before, knowledge of the covariancematrix allows the evaluation of entanglement, purity, andfidelity measures. The analytical results are easily obtained
but sufficiently complicated that we do not quote them here,
while the numerical results coincide with those previously
obtained in the adiabatic limit.
B. Solution with time-dependent Hamiltonian

The results of Sec. VI A are still only valid provided that
we are justified in making a rotating-wave approximation;that is, in discarding the time-dependent (“counter-rotating”)
contributions that arise in the derivation of the Hamiltonian (8).
Here, we account for these time-dependent contributions;the exact forms that they take are given in Appendix D 1.
The corresponding contributions to the Heisenberg-Langevinequations can be handled by making the replacement A0 →A(t) in Eq. (46), where the time-dependent drift matrix is given
byN(Ak+ e+2iδk t + Ak− e−2iδk t ).

A(t) = A0 +

(48)

k=1

For the case of two-tone driving (10), N = 1 and δ1 = ωm .
For four-tone driving (16), we have N = 4 and δ1 = δ, δ2 =
ωb + , δ3 = ωm , and δ4 = ωa −  [see Fig. 2 and recall thatδ was introduced in Eq. (23)]. The first of these contributions
is due to the second drive on the same side of the cavity
resonance frequency, while the remaining contributions are
due to the drives on the other side of the cavity resonancefrequency. The matrices Ak± are given in Appendix D 2.
Given that the drift matrix (48) is now time varying, thecovariance matrix V(t) is given by solution of the Lyapunovlike differential equation

063805-8

V̇ = AV + VA† + B0 BT0 .

(49)

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)

FIG. 7. (Color online) The effects of counter-rotating contributions to the Hamiltonian, described in detail in Appendix D, on theentanglement and purity of the mechanical two-mode steady state under four-tone driving. The plots show (a) the Duan quantity, (b) the
logarithmic negativity, (c) the occupations of the Bogoliubov modes, and (d) the purity, as functions of the drive asymmetry. The solid black
curve corresponds to the time-independent Hamiltonian (8), while the other curves correspond to cases where counter-rotating effects are
significant. The blue curve (long dashes) is for δ/κ = 0.5, the red curve (short dashes) is for ω1 /κ = 5, and the purple curve (long and shortdashes) is for ga /gb = 1.5. Although it is possible for counter-rotating contributions to have a significant effect on these measures, it should also
be possible to achieve sufficiently high sideband resolution that their effects may be neglected. Note that under two-tone driving we typicallyhave ga = gb and always have δ = 0; the effect of the parameter ωm /κ in that case is comparable to the effect of the parameter ω1 /κ in thefour-tone driving case. Parameters, unless otherwise specified, for each curve are δ/κ = 1, ω1 /κ = 100, and d = 1. The parameters common
5
6 −1
−5
to each curve are n̄ = 0, Gm± = 0, C− = 1200, κ = 2π × 1.592 × 10 Hz = 10 s , γ /κ = 4 × 10 , and /κ = 0.1.

The covariance matrix will be oscillatory in the long-time
limit; we seek the dc component of its solution. The direct
numerical solution of Eq. (49) is inefficient, so instead we usean ansatz to obtain an approximate numerical solution [40].
The procedure used is outlined in Appendix D 3.
C. Effects of counter-rotating terms

The effects of the counter-rotating Hamiltonian contributions on the entanglement and purity of the steady state areshown in Figs. 7 and 8. Figure 7 shows these as functions ofdrive asymmetry, while Fig. 8 shows the two-mode squeezing(entanglement), optimized over the drive asymmetry, as afunction of the cooperativity parameter C− .
From Fig. 7 it is clear that there is a degradation in theentanglement and purity measures as the frequency of thecounter-rotating terms is lowered. This is unsurprising since
the form of the time-dependent Hamiltonian contributions,detailed in Appendix D 1, depart from the ideal form of Eq. (7).
However, the same overall behavior in the entanglement andpurity is observed; that is, a maximum in the entanglement anda monotonic decrease in the purity.
The rotating-wave approximation results coincide withthe full time-dependent Hamiltonian results in the limit that
all counter-rotating frequencies greatly exceed the cavity
decay rate; that is, all |δk |κ. With the parameters chosen,counter-rotating effects become significant at δ/κ ∼ 0.5 andωa /κ ∼ 5; these correspond to modest sideband resolutions.

FIG. 8. (Color online) Mechanical two-mode squeezing as afunction of the cooperativity parameter C− , introduced in Eq. (30),for a range of sideband resolutions ωm /κ. The results are presented
for the case of two-tone cavity driving with the single-photon
optomechanical coupling rates being equal, ga = gb (that is, no imperfection in the effective coupling rates). The quantity plotted is thetwo-mode squeezing in dB, defined by TMS (dB) ≡ − log10 [(X̂2+  +P̂2− )/(X̂2+  + P̂2− )0 ], with both the numerator and denominatorbeing instances of the Duan quantity of Eq. (28). At each value ofthe cooperativity parameter C− , the Duan quantity is minimized overthe effective coupling asymmetry G+ /G− . The curves are shown fora rotating-wave approximation (solid black curve) meaning that the
sideband resolution is effectively infinite; ωm /κ = 103 (blue curvewith long dashes); ωm /κ = 102 (purple curve with short dashes);ωm /κ = 10 (red curve with long and short dashes). In the case offour-tone cavity driving, the behavior of the two-mode squeezing
as a function of ω1 /κ is similar to the behavior seen here as afunction of ωm /κ. The other parameters are as specified in the captionof Fig. 7.

063805-9

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

It should be possible to significantly exceed this resolution
and therefore largely avoid the effects of counter-rotating
terms. The ratio δ/κ can be reduced further than the ratioωa /κ without significant deleterious effects due to the distinct
manner with which the corresponding contributions enter thefull time-dependent Hamiltonian. Provided that the asymmetry
in the single-photon optomechanical coupling rates is small,the effective couplings associated with the terms rotating at±2δ are not exponentially enhanced in the large-r limit, while
those oscillating at ±2ωa , ±2ωm , and ±2ωb are exponentiallyenhanced [see Eq. (D3)].
Note that with counter-rotating terms included, the entanglement and purity of the mechanical steady state depend on
the asymmetry in the single-photon optomechanical coupling
rates, even if there are no imperfections in the effective
couplings. This is in contrast to the results obtained withthe time-independent Hamiltonian (8). As seen in Fig. 7, thisasymmetry leads to a significant degradation in entanglement
and purity measures if the imperfection is around ∼50% of
the effective coupling. Again, it should be possible to engineer
the optomechanical system such that the asymmetry is much
lower than this value, and the corresponding deleterious effects
are negligible.
The behavior of the optimized mechanical two-mode
squeezing as a function of the cooperativity parameter C− isshown in Fig. 8. Even neglecting the effects of counter-rotating
terms, the achievable two-mode squeezing (entanglement)
plateaus in the large-cooperativity limit. For a finite sideband
resolution, however, the squeezing goes through a maximum as
a function of the cooperativity, with the maximum occurring
at a lower value of the cooperativity parameter as the sideband resolution is decreased. Unsurprisingly, the discrepancy
between the rotating-wave approximation (RWA) result andthe full result increases as the sideband resolution is reduced.
Note, however, that these discrepancies become significant
only at very large values of the cooperativity parameter.
Also note that in the large-cooperativity limit, it is possible
for the dynamics associated with the full time-dependent
Hamiltonian to be unstable where the dynamics associatedwith the corresponding time-independent Hamiltonian arestable. However, this does not tend to be the case at the
optimal drive asymmetry. This potential for the onset of an
instability occurs at higher values of the cooperativity than
we have previously considered in this work. Our results showthat for high levels of steady-state entanglement, the reservoir
engineering scheme discussed here is robust against realistic
levels of counter-rotating corrections.
VII. EXPERIMENTAL OBSERVABILITY
A. Output spectrum

From an experimental point of view, reconstructing theentire covariance matrix would be extremely demanding. Even
performing direct measurements of both of the collective
quadatures required for testing the Duan criterion of Eq. (28)
would be difficult. However, one could perform a backactionevading measurement of the collective quadrature X̂+ , and takethis as some evidence for the existence of two-mode squeezingin the steady state. Alternatively, we can seek a signature of
the mechanical entanglement in the cavity output spectrum.

As usual, the output spectrum is calculated by first solving
the Heisenberg-Langevin equations in the frequency domain.
Taking the Fourier transform of Eq. (46) we findX[ω]= −(A0 + iωI6 )−1 · B0 · X[ω],(50)
where X[ω]= (â[ω],â † [ω],b̂[ω],b̂† [ω],ĉ[ω],ĉ† [ω])T and thematrices are given in Appendix C 3. The output spectrum is
calculated in the standard manner [33] as†S[ω] = dt eiωt δ ĉout (t)δ ĉout (0),(51)
where the output√cavity field is given by δ ĉout = ĉout − ĉout and ĉout = ĉin + κ ĉin .
We first calculate the spectrum assuming that there are
no imperfections in the effective couplings (Gm± = 0) andignoring time-dependent Hamiltonian contributions; that is,with the effective Hamiltonian (8). Then, the cavity output
spectrum, in the limit γ /  → 0, is given byS[ω] = κ32[G2− n̄ + G2+ (n̄ + 1)]γ [γ 2 + 4(ω2 + 2 )],|N (ω)|2

(52)

where
N (ω) = [8G 2 + (γ − 2iω)(κ − 2iω)](γ − 2iω) +42 (κ − 2iω). In the case  = 0, this reduces to the result
for a single mechanical oscillator [17]. As shown in Fig. 9, theoutput spectrum exhibits peaks at detunings around ± fromthe cavity resonance frequency. This corresponds to the drive
photons being scattered towards the cavity resonance, withan energy ωa or ωb being provided by or extracted from themechanical oscillators. As the optomechanical damping rateis increased, the widths of the spectral peaks increase and they
are shifted to larger detunings (for G− /G+ > 0), as expected.
If we now allow for the possibility of imperfections in theeffective couplings, but still ignore time-dependent contributions, the cavity output spectrum is again readily obtained. The
corresponding spectra are also shown in Fig. 9. The general
expression is complicated, but at detunings of ± we findS[±] = γ κ2m 2(G− ± Gm− ) n̄ + (G+ ± G+ ) (1 + n̄)
.
2
m 2 2
2
[G2− − (Gm− ) − G+ + (G+ ) ](53)

Clearly, the asymmetry in the spectral peaks is determinedby the imperfections Gm± in the effective optomechanicalcouplings.
B. Bogoliubov modes

Knowledge of the cavity output spectrum can be used toprovide us with information about the occupations of themechanical Bogoliubov modes. In particular, neglecting imperfections in the effective couplings and in the limit γ /  →0, we can show that the occupations of the Bogoliubov modes
are related to the integral of each peak in the output spectrum by 0 +∞S[ω]dω =S[ω]dω
−∞0= 8π κ

G2†
β̂ β̂i ,
4G 2 + κ(κ + γ ) i(54)

for i = 1 or 2 (in this, and subsequent, expressions). Therefore,
from the output spectrum and knowledge of the system parameters, one can determine the occupations of the Bogoliubov

063805-10

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)

FIG. 9. (Color online) Cavity output spectra, as defined in Eq. (51), centered around detunings from the cavity resonance frequency of (leftpanel) ω = − and (right panel) ω = +. The spectra are normalized to the peak in the output spectrum in the absence of imperfections in theeffective couplings, denoted by S0 [ω]|max . The spectra are shown for the case (solid black curve) without imperfections in the effective couplingsmm(Gm± = 0), and for the cases G± /G± = 0.3 (blue curve with long dashes) and G± /G± = −0.3 (red curve with short dashes). Imperfections in
the effective couplings lead to asymmetry in the observed output spectra [see Eq. (53)]. In the absence of these imperfections, the steady-statemechanical entanglement can be bounded based on a measurement of the output spectrum. The spectra are shown for a drive asymmetry
G+ /G− = 0.9 and a cooperativity C− = 1200, while other parameters are n̄ = 0, κ = 2π × 1.592 × 105 Hz = 106 s−1 , γ /κ = 4 × 10−5 , and/κ = 0.1.

modes. The same information can be obtained from the heightsof the spectral peaks since (in the same limit) we also have
S[±] =γ κ + 4(G2− − G2+ ) †β̂i β̂i .
G2− − G2+

(55)

Note that Eqs. (54) and (55) only hold when the imperfectionsin the effective couplings are less than ∼ 1%. A similar result
to that of Eq. (54) was obtained for a single mechanicaloscillator [17], although in that case the integration is over all
frequencies.

from the perspective of quantum information processing, toconsider the opposite scenario in which there are two cavity
modes coupled to a single mechanical oscillator [42–45]. In
particular, the ubiquity of optomechanical couplings raises
the possibility of entangling cavity modes of vastly different
frequencies (e.g., microwave and optical modes). This threemode optomechanical system is depicted schematically in
Fig. 10(a).
Now, the Hamiltonian of the system is [cf. Eq. (1)]Ĥ = ωa â † â + ωb b̂† b̂ + ωc ĉ† ĉ + ga (ĉ + ĉ† )â † â
+ gb (ĉ + ĉ† )b̂† b̂ + Ĥdrive + Ĥdiss ,C. Entanglement criterion

Now, from Eqs. (54) and (55) it is clear that we can estimate
the occupations of the Bogoliubov modes using the cavity
output spectrum. Recall that the simplest means of verifying
the presence of mechanical entanglement is via the Duancriterion (28). The task then is to bound the Duan quantity
using our knowledge of the occupations of the Bogoliubov
modes.
Repeated application of the generalized Cauchy-Schwarz†
inequality [17,41] allows one to show that |β̂i2 |  β̂i β̂i  +††1. With the additional assumption that β̂1 β̂1  = β̂2 β̂2 ,
2known to be true in the absence of imperfections in the
couplings and in the limit γ /  → 0, we can bound the Duanquantity. Explicitly, we find that
†

2
X̂+
 + P̂−2   8e−2r (β̂i β̂i  + 1/2).

(57)

where â and b̂ describe two cavity modes, and ĉ describes themechanical oscillator. We shall assume a driving Hamiltonianof the formĤdrive = (Ea+ e+iωc t + Ea− e−iωc t )e+i(ωa −)t â
+ (Eb+ e+iωc t + Eb− e−iωc t )e+i(ωb +)t b̂ + H.c. (58)(56)

The parameter r is known from the drive asymmetry [cf.
Eq. (9b)]. The Duan quantity and its bound converge in thehighly entangled (large-r) regime, and therefore we expect it
to reliably indicate the existence of an entangled mechanicalsteady state.
VIII. TWO CAVITY MODES, ONE MECHANICALOSCILLATOR

Thus far, we have considered a three-mode optomechanicalsystem composed of two mechanical oscillators coupled to one
common cavity mode. It is of considerable interest, particularly

FIG. 10. (Color online) (a) Representation of a three-mode optomechanical system composed of two cavity modes (with resonancefrequencies ωa and ωb ), each independently coupled by radiationpressure to a single mechanical oscillator (with resonance frequency
ωc ). (b) The driving conditions, in terms of frequencies, that lead to thelinearized effective Hamiltonian (8) with the role of the cavity andmechanical modes interchanged. The cavity resonance frequencies
are indicated by the blue lines, while the driving tones are indicatedby red lines. The driving tones are placed symmetrically (at detunings±ωc ) about detunings of  from the cavity resonance frequencies.

063805-11

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

That is, driving at ωa −  ± ωc and ωb +  ± ωc for cavity
modes â and b̂, respectively, as depicted in Fig. 10(b). In
an interaction picture defined with respect to the Hamiltonian(4), the Hamiltonian (57) takes the form of the Hamiltonian (8)if we set ga ā+ = gb b̄+ ≡ G+ and ga ā− = gb b̄− ≡ G− , withā± and b̄± denoting the steady-state amplitudes at the driven
cavity sidebands. Thus, we can realize the same physics withthe cavity modes that we described previously for the mechanical modes, including the possibility of generating highly
pure, highly entangled electromagnetic modes. As before, thisderivation relies on the assumptions that we are operatingin the resolved-sideband regime, that the driving strengths at
the driven sidebands are large, and that the effects of counterrotating terms are negligible. The deviations arising from abandoning any of these assumptions may, of course, be calculated.

ACKNOWLEDGMENTS

This work was supported by NSERC, an ECR Grant fromUNSW Canberra, and the DARPA ORCHID program undera grant from the AFOSR. We thank J. Suh, C. U. Lei,
I. Petersen, and S. Armstrong for useful discussions.
APPENDIX A: DERIVATION OF HAMILTONIAN

We consider the three-mode optomechanical system depicted in Figs. 1(a) and 2(a), and start from the Hamiltonian (1)with the four-tone drive (16); the scenario with the two-tone
drive (10) shall be treated as a special case. Moving into arotating frame with respect to Ĥ0 = ω1 â † â + ω2 b̂† b̂ + ωc ĉ† ĉ,recalling Eqs. (18a) and (18b), we obtainĤ = (â † â − b̂† b̂) + ga (âe−iω1 t + â † e+iω1 t )ĉ† ĉ

IX. CONCLUSIONS

We have provided a detailed proposal for configuring athree-mode optomechanical system such that the steady stateincludes a highly pure, highly entangled two-mode squeezed
state. The generation of both mechanical and electromagnetic
two-mode squeezed states has been described. The symmetryof this steady state makes it an attractive platform for the
implementation of continuous-variable teleportation protocols. The proposal is efficient in the sense that it requires
only one driven auxiliary mode to be configured as theengineered reservoir. The problem of unequal single-photon
optomechanical couplings has been overcome by using afour-tone driving scheme, and potential instabilities arisingfrom counter-rotating terms have been accounted for. A
simple experimental signature for the presence of mechanicalentanglement, in terms of the cavity output spectrum, hasbeen provided. The proposal described is implementable withexisting technology.

+ gb (b̂e−iω2 t + b̂† e+iω2 t )ĉ† ĉ + Ĥdrive + Ĥdiss .

(A1)

The effective oscillation frequency  is chosen such that
Eqs. (22a) and (22b) are satisfied. The Heisenberg equations,neglecting noise terms, corresponding to Eq. (A1) are
γaâ˙ = −iâ − iga e+iω1 t ĉ† ĉ − â,
(A2a)
2
γ
b
b̂˙ = ib̂ − igb e+iω2 t ĉ† ĉ − b̂,(A2b)2ĉ˙ = −iga (âe−iω1 t + â † e+iω1 t )ĉ − iE1+ e−iω1 t− iE1− e+iω1 t − igb (b̂e−iω2 t + b̂† e+iω2 t )ĉκ
− iE2+ e−iω2 t − iE2− e+iω2 t − ĉ.
2Assuming resolved-sideband operation (ωa(b)the ansatz [35]ĉ(t) = ĉ0 (t) + ĉ1+ (t)e−iω1 t + ĉ1− (t)e+iω1 t + ĉ2+ (t)e−iω2 t + ĉ2− (t)e+iω2 t .

(A2c)κ), we take

(A3)

Then, substituting Eq. (A3) into the system (A2a)–(A2c) and separating out the Fourier coefficients of the cavity field we obtainthe system
γa
†
†
†
†
†
†
â˙ = −iâ − iga (ĉ1− ĉ0 + ĉ1+ ĉ0 ) − â − iga e+2iδt (ĉ2− ĉ0 + ĉ2+ ĉ0 ) − iga e+2iω1 t (ĉ1+ ĉ0 + ĉ1− ĉ0 )2††− iga e+2iωm t (ĉ2+ ĉ0 + ĉ2− ĉ0 ),

(A4a)
γ
b
†
†
†
†
†
†
†
†
b̂˙ = ib̂ − igb (ĉ2− ĉ0 + ĉ2+ ĉ0 ) − b̂ − igb e−2iδt (ĉ1− ĉ0 + ĉ1+ ĉ0 ) − igb e+2iωm t (ĉ1+ ĉ0 + ĉ1− ĉ0 ) − igb e+2iω2 t (ĉ2+ ĉ0 + ĉ2− ĉ0 ),
2
(A4b)
κ
ĉ˙0 = −iga â ĉ1− − iga â † ĉ1+ − igb b̂ĉ2− − igb b̂† ĉ2+ − ĉ0 − iga (â ĉ2− e−2iδt + â † ĉ2+ e+2iδt ) − igb (b̂ĉ1− e+2iδt + b̂† ĉ1+ e−2iδt )2−2iω1 t−2iωm t†+2iω1 t+ â ĉ2+ e+ â ĉ1− e+ â † ĉ2− e+2iωm t ) − igb (b̂ĉ1+ e−2iωm t + b̂ĉ2+ e−2iω2 t− iga (â ĉ1+ e+ b̂† ĉ1− e+2iωm t + b̂† ĉ2− e+2iω2 t ),

(A4c)

†

ĉ˙1− = (−iω1 − κ/2)ĉ1− − iga â ĉ0 − iE1− ,ĉ˙1+ = (iω1 − κ/2)ĉ1− − iga â ĉ0 − iE1+ ,(A4d)ĉ˙2− = (−iω2 − κ/2)ĉ1− − igb b̂† ĉ0 − iE2− ,ĉ˙2+ = (iω2 − κ/2)ĉ1− − igb b̂ĉ0 − iE2+ .

(A4f)

(A4e)
(A4g)
063805-12

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)

Note that in writing out Eqs. (A4a)–(A4g) we have retained fast-rotating terms only in Eqs. (A4a)–(A4c). Solving
Eqs. (A4d)–(A4g) for the steady-state amplitudes of the fieldat the driven sidebands, assuming that the single-photon
optomechanical couplings are relatively small, we obtain the
results of Eq. (17).
Replacing the operator Fourier components at the driven
sidebands by their classical steady-state values, we can write
an effective (quadratic) Hamiltonian for the system dynamics.
If we neglect fast-rotating terms in Eqs. (A4a)–(A4c), thisHamiltonian will be time independent; the time-dependent
contributions to the effective Hamiltonian are given in
Appendix D 1.
With four-tone driving, the effective Hamiltonian is (replacing ĉ0 → ĉ)Ĥ = (â † â − b̂† b̂) + ga [(c̄1− â + c̄1+ â † )ĉ† + H.c.]+ gb [(c̄2− b̂ + c̄2+ b̂† )ĉ† + H.c.] + Ĥdiss .

(A5)

Assuming that the drives are matched according to Eq. (19)and that the steady-state amplitudes in the driven sidebands
are real, we get the Hamiltonian (8) with the effective
couplings of Eq. (20). If there is an imperfection (drive mismatch), then we have the additional Hamiltonian contributions
given by Eq. (14) with the coupling imperfections given by
Eq. (15).
With two-tone driving we have ω1 = ω2 = ωm , and theappropriate effective Hamiltonian is now (A5) with c̄k±replaced by c̄± , the driving strengths at the frequencies
ωc ± ωm . If the single-photon optomechanical couplings areequal, then the effective Hamiltonian is given by Eq. (8)where the effective couplings are given by Eq. (12). Imperfections (unequal single-photon optomechanical couplings)
lead to additional contributions of the form of Eq. (14),where the effective coupling imperfections are given by
Eq. (12).
APPENDIX B: TWO-MODE GAUSSIAN STATES
1. Entanglement

The entanglement of a two-mode Gaussian state may bequantified from its symmetrically ordered covariance matrix V
via the logarithmic negativity [46,47]. Writing the covariancematrix as in Eq. (37), the logarithmic negativity is given by
EN = max {0, − ln 2η} ,whereη = 2−1/2 {(V) − [(V)2 − 4 det V]1/2 }1/2
(V) = det Vb + det Va − 2 detVab .

(B1)
and

where the coefficients are given bya(b)2+ 1/2 cosh2 ξ + n̄b(a)ca(b) = n̄thth + 1/2 sinh ξ, (B3a)cab = n̄ath + n̄bth + 1 sinh ξ cosh ξ.

APPENDIX C: HEISENBERG-LANGEVIN EQUATIONS
1. Collective quadratures

In the adiabatic limit and in terms of the collective
quadratures of Eqs. (26a) and (26b), the dynamics of thetwo-mode mechanical system is described by Eq. (44), withthe matrices


A+−
A+
,
(C1a)
A0 =
A+−
A−


B1+ B1−
B1 =
,
(C1b)
B1+ B1−
 
√ I2
B2 =
,
(C1c)
022
including the components

Va = ca I2 , Vb = cb I2 , Vab = −cab σz ,A+ = −(γ /2 + )I2 ,(C2a)A− = −(γ /2)I2 ,
B1± = γ (1 ± l)/2I2 ,−lγ /2A+− =,−
−lγ /2

(C2b)
(C2c)
(C2d)

where l = (γa − γb )/(2γ ) and γ = (γa + γb )/2.
2. Individual quadratures

The dynamics of the linearized, three-mode optomechanical system, with Hamiltonian (8), is described by Eq. (46).
The system matrix is⎤⎡
Aa 022 Ca(C3)A0 = ⎣022 Ab Cb ⎦,Ca Cb Ac
where 022 is the 2 × 2 zero matrix, and−γa /2,Aa =−
−γa /2−−γb /2Ab =,−γb /2
Ac = −(κ/2)I2 ,0Ca =−G− − G+ + Gm
d

0
Cb =
−G− − G+ − Gmd2. Thermal two-mode squeezed state

An alternative method for characterizing the mechanicaltwo-mode state is obtained by noting that the covariancematrix of our steady state takes the form of a thermal
two-mode squeezed state [48], as defined in Eq. (39). This
state is described by three parameters: a two-mode squeezing
parameter ξ , and the two initial thermal occupations n̄ath andn̄bth . The covariance matrix for the thermal two-mode squeezed
state, in the block form of Eq. (37), is

(B3b)


G− − G+ − Gms ,0
G− − G+ + Gm
s ,
0

(C4a)
(C4b)
(C4c)
(C4d)
(C4e)

m
m
with the shorthand notation Gm
ands = G− + G+mmm
Gd = G− − G+ . The noise matrix is given by
B0 = Diag( γa (n̄a + 1/2)I2 , γb (n̄b + 1/2)I2 , κ/2I2 ).

(B2)063805-13(C5)

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)



−G− − Gm−
Cb = iG+ − Gm+

3. Individual mode operators

For the purpose of evaluating the cavity output spectrum,
it is more convenient to work in terms of annihilation andcreation operators, as in Eq. (50). The corresponding matrices
are
√
√
√
B0 = Diag( γa I2 , γb I2 , κI2 ),(C6)
while A0 is given by the block matrix form of Eq. (C3), now
with


0
−i − γa /2
,
(C7a)
Aa =
0
i − γa /20i − γb /2
Ab =
,
(C7b)
0
−i − γb /2
−G+ − Gm
−G− + Gm
−
+
,
(C7c)
Ca = i
G+ + Gm
G− − Gm+−

−G+ + Gm
+ ,G− + Gm−(C7d)while Ac is still given by Eq. (C4c).
APPENDIX D: COUNTER-ROTATING CONTRIBUTIONS
1. Hamiltonians

In deriving the time-independent Hamiltonian (8) we
discarded fast-rotating terms; here, we include them. The fulltime-dependent Hamiltonian isĤ(t) = Ĥ + ĤCR ,(D1)

where Ĥ is the time-independent effective Hamiltonian (8) andĤCR is the time-dependent (“counter-rotating”) contribution.

a. Four-tone driving

With four driving tones (16), the time-dependent part of the Hamiltonian (D1) is
ĤCR = ga {â[c̄2− e−2iδt + c̄2+ e−2iωm t + c̄1+ e−2iω1 t ] + â † [c̄2+ e+2iδt + c̄1− e+2iω1 t + c̄2− e2iωm t ]}ĉ†+ gb {b̂[c̄1− e+2iδt + c̄1+ e−2iωm t + c̄2+ e−2iω2 t ] + b̂† [c̄1+ e−2iδt + c̄2− e+2iω2 t + c̄1− e+2iωm t ]}ĉ† + H.c.

(D2)

There are terms at four distinct oscillation frequency magnitudes: the ±2δ terms are associated with the two drives being on
the same side of the cavity resonance frequency, while the terms oscillating at ±2(ωa − ), ± 2ωm , ± 2(ωb + ) are associatedwith two drives on opposing sides of the cavity resonance frequency. In terms of Bogoliubov modes, the Hamiltonian (D2) may
be written as†

ĤCR = Ge−2iδt {β̂1 ĉ† [1/d + (d − 1/d) cosh2 r] − β̂2 ĉ† (d − 1/d) cosh r sinh r}†+ Ge+2iδt {β̂2 ĉ† [1/d − (d − 1/d) sinh2 r] + β̂1 ĉ† (d − 1/d) cosh r sinh r}
1 +i(ω1 +ω2 )t†+2iω2 t−2iω1 t−i(ω1 +ω2 )t− G cosh r sinh r β̂1 ĉ e
−e
+ e
− de
d


1
+ β̂2 ĉ† e+2iω1 t − e−2iω2 t + de+i(ω1 +ω2 )t − e−i(ω1 +ω2 )td
1 −i(ω1 +ω2 )t† †2
2iω1 ti(ω1 +ω2 )t2−2iω2 te
+ de) − sinh r+e+ G β̂1 ĉ cosh r(e
d


1
†
+ G β̂2 ĉ† cosh2 r e2iω2 t + ei(ω1 +ω2 )t − sinh2 r(de−i(ω1 +ω2 )t + e−2iω1 t ) + H.c.,dwhere the asymmetry in the single-photon optomechanicalcoupling rates is parametrized by
d≡

ga
.
gb

(D4)

b. Two-tone driving

2. Drift matrix

A time-dependent Hamiltonian (D1) leads to a timedependent drift matrix in the corresponding HeisenbergLangevin equations [see Eq. (46)]. The drift matrix takes theform given in Eq. (48), and we specify the coefficient matrices
here. In writing out these matrices it is useful to parametrize
the coupling imperfection by ε± , where

With two driving tones, as per Eq. (10), the time-dependent
contribution to the Hamiltonian (D1) is
ĤCR = G+ (â + b̂)e
†

c̄1±
gb
ε± ≡ .
c̄2±ga−2iωm t †ĉ†

+ G− (â + b̂ )e+2iωm t †ĉ + H.c.,

only containing fast-rotating terms oscillating at |2ωm |.

(D5)

(D3)

(D6)

Having both sets of drives matched according to Eq. (19),and therefore no imperfection in the effective coupling,corresponds to ε± = 1. The drift coefficient matrices, for the063805-14

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)


i(−G̃− ε̃+ ∓ G̃+ ε̃− )Q̃± =−G̃− ε̃+ − G̃+ ε̃−

general case of four-tone driving, are
⎤
⎡
022
d M̃+
1 ⎣ 022022
022
N+ /d ⎦ ,A1+ =2 d M̃
N
/d
022
−
−
⎤
⎡
022 022 022
1
A2+ = ⎣022 022 Q̃+ ⎦ ,2 0
Q̃− 022
22
⎡
⎤
022
022
d Q̃+
1⎢
⎥
022
Q+ /d ⎦ ,A3+ = ⎣ 022
2
d Q̃− Q− /d022⎡⎤022 022 Q̃+1⎢⎥
A4+ = ⎣ 022 022 022 ⎦ ,2Q̃− 022 022

with Q± given by Q̃± with the replacement G̃± → G± . The
scalar tilde quantities are defined as(D7b)

G̃± = G± ε± ,2(1 + ε± )ε̃± =.
(1 + ε+ )(1 + ε− )

(D7c)

(D9a)
(D9b)

3. Time-dependent drift matrix: Solution(D7d)

G− ε̃+ − G+ ε̃−i(±G− ε̃+ − G+ ε̃− )

(D8a)

,
(D8b)

N

(D8c)

(D7a)

where Ak− = A∗k+ , 022 is the 2 × 2 zero matrix, and we haveintroduced the notationG̃− ε̃+ − G̃+ ε̃−i(∓G̃− ε̃+ − G̃+ ε̃− ),M̃± =−G̃− ε̃+ − G̃+ ε̃−i(∓G̃− ε̃+ + G̃+ ε̃− )i(±G− ε̃+ + G+ ε̃− )N± =−G− ε̃+ − G+ ε̃−−G̃− ε̃+ + G̃+ ε̃−,i(G̃− ε̃+ ∓ G̃+ ε̃− )Given the form of the drift matrix (48), we expect thecovariance matrix V, given by the solution of (49), to be
oscillatory in the long-time limit. We approximate the solution
via the covariance matrix ansatz [40]N

(Vk+ e+2iδk t + Vk− e−2iδk t ).

V(t) = V0 +

(D10)

k=1

In general, the solution will contain harmonics of thebare frequencies that appear in Eq. (48), as well as theirsum and difference frequencies. However, the solution that
we really seek is the dc component of the covariancematrix V0 .
The equations of motion (46), with the drift matrix (48),
may be written in the frequency domain as in [ω] ≡ N [ω], − 2δk ] + Ak− X[ω + 2δk ]) + (A0 − iωI6 )X[ω](Ak+ X[ω= −B · X

(D11)

k=1

 +∞
where we have defined the Fourier transform as F[ω] = −∞ f (t) e−iωt dt. Now, we form the frequency-dependent state and noise
vectors. These are (2N + 1)-dimensional vectors where N is the number of positive-frequency counter-rotating terms in Eq. (48).
 − 2δ|N+1−n| sgn (N + 1 − n)] and N [ω − 2δ|N+1−n| sgn (N + 1 − n)], respectively. Subsequently, we
The nth elements are X[ωcan write the linear systemĀ[ω] · X[ω]= N[ω],
where⎡

A0 − i(ω − 2δN )I6

⎢
⎢
⎢
⎢
⎢
⎢
Ā[ω] = ⎢
⎢
⎢
⎢
⎢
⎢
⎣

..

⎤

AN−
..
.

.
..

AN+

(D12)

.
A1+

...

A1−
A0 − iωI6
A1+
..
.
AN+

A1−
..
.

...
..

AN−

.

⎥
⎥
⎥
⎥
⎥
⎥
⎥.
⎥
⎥
⎥
⎥
⎥
⎦

(D13)

A0 − i(ω + 2δN )I6

Introducing the noise correlation matrix i,j [ω,ω ] = Ni [ω]Nj∗ [ω ] + Nj∗ [ω ]Ni [ω]/2 leads toN[ω,ω ] = D0 δ[ω − ω ] +(Dk+ δ[ω − ω − 2δk ] + Dk− δ[ω − ω + 2δk ]),

(D14)

k=1

where the matrices are defined by (D0 )ii = BBT , (Dn+ )N+1−n,N+1 = (Dn+ )N+1,N+1+n = BBT for n ∈ {1, . . . ,N}, and Dn− =DTn+ . The indices refer to 6 × 6 blocks in the overall matrix.
063805-15

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 89, 063805 (2014)

Solving the linear system (D12) leads toV[ω,ω ] = Ā−1 [ω][ω,ω ](Ā−1 [ω ])† .

(D15)

We really want the 6 × 6 central block which we denote Ṽ[ω,ω ]. This is given byNṼ[ω,ω ] = Ṽ0 δ[ω − ω ] +(Ṽk+ [ω]δ[ω − ω − 2δk ] + Ṽk− [ω]δ[ω − ω + 2δk ]),

(D16)

k=1

where we have the coefficientsṼ0 [ω] = [Ā−1 [ω]D0 (Ā−1 [ω])† ]6 ,−1

−1

(D17a)†

Ṽk± [ω] = [Ā [ω]Dk± (Ā [ω ∓ 2δk ]) ]6 .
The coefficients in Eq. (D10) follow from
1
V0 =
2π
Vk± =

1
2π

(D17b)

 +∞−∞Ṽ0 [ω]dω,(D18a)

Ṽk± [ω]dω.

(D18b)

 +∞[1] G. J. Milburn and M. J. Woolley, Acta Phys. Slovaca 61, 483(2012).
[2] S. Mancini, V. Giovannetti, D. Vitali, and P. Tombesi, Phys. Rev.
Lett. 88, 120401 (2002).
[3] M. Schmidt, M. Ludwig, and F. Marquardt, New J. Phys. 14,125005 (2012).
[4] H. Tan, L. F. Buchmann, H. Seok, and G. Li, Phys. Rev. A 87,
022318 (2013).
[5] M. J. Woolley and A. A. Clerk, Phys. Rev. A 87, 063846 (2013).
[6] A. Szorkovszky, A. A. Clerk, A. C. Doherty, and W. P. Bowen,
arXiv:1402.6392.
[7] C. A. Muschik, E. S. Polzik, and J. I. Cirac, Phys. Rev. A 83,
052312 (2011).
[8] Y.-D. Wang and A. A. Clerk, Phys. Rev. Lett. 110, 253601(2013).
[9] H. Tan, G. Li, and P. Meystre, Phys. Rev. A 87, 033829 (2013).
[10] A. Tomadin, S. Diehl, M. D. Lukin, P. Rabl, and P. Zoller, Phys.
Rev. A 86, 033821 (2012).
[11] J. B. Hertzberg, T. Rocheleau, T. Ndukum, M. Savva, A. A.
Clerk, and K. C. Schwab, Nat. Phys. 6, 213 (2010).
[12] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman,
K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W.
Simmonds, Nature (London) 475, 359 (2011).
[13] F. Massel, S. U. Cho, J.-M. Pirkkalainen, P. J. Hakonen, T. T.
Heikkilä, and M. A. Sillanpää, Nat. Commun. 3, 987 (2012).
[14] J. Suh, A. J. Weinstein, C. U. Lei, E. E. Wollman, S. K. Steinke,
P. Meystre, A. A. Clerk, and K. C. Schwab, arXiv:1312.4084.
[15] H. Bruus and K. Flensberg, Many-Body Quantum Theory inCondensed Matter Physics (Oxford University Press, Oxford,
UK, 2004).
[16] J. Kerckhoff, R. W. Andrews, H. S. Ku, W. F. Kindel, K. Cicak,
R. W. Simmonds, and K. W. Lehnert, Phys. Rev. X 3, 021013(2013).
[17] A. Kronwald, F. Marquardt, and A. A. Clerk, Phys. Rev. A 88,
063833 (2013).

−∞[18] S. L. Braunstein and H. J. Kimble, Phys. Rev. Lett. 80, 869(1998).
[19] G. Adesso and F. Illuminati, Phys. Rev. Lett. 95, 150503(2005).
[20] J. I. Cirac, A. S. Parkins, R. Blatt, and P. Zoller, Phys. Rev. Lett.
70, 556 (1993).
[21] J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 77, 4728
(1996).
[22] M. B. Plenio and S. F. Huelga, Phys. Rev. Lett. 88, 197901(2002).
[23] A. S. Parkins, E. Solano, and J. I. Cirac, Phys. Rev. Lett. 96,
053602 (2006).
[24] E. G. Dalla Torre, J. Otterbach, E. Demler, V. Vuletic, andM. D. Lukin, Phys. Rev. Lett. 110, 120402 (2013).
[25] Q. He and Z. Ficek, Phys. Rev. A 89, 022332 (2014).
[26] H. Krauter, C. A. Muschik, K. Jensen, W. Wasilewski, J. M.
Petersen, J. I. Cirac, and E. S. Polzik, Phys. Rev. Lett. 107,
080503 (2011).
[27] K. W. Murch, U. Vool, D. Zhou, S. J. Weber, S. M. Girvin, and
I. Siddiqi, Phys. Rev. Lett. 109, 183602 (2012).
[28] Y. Lin, J. P. Gaebler, F. Reiter, T. R. Tan, R. Bowler, A. S.
Sorenson, D. Leibfried, and D. J. Wineland, Nature (London)
504, 415 (2013).
[29] S. Shankar, M. Hatridge, Z. Leghtas, K. M. Sliwa, A. Narla,
U. Vool, S. M. Girvin, L. Frunzio, M. Mirrahimi, and M. H.
Devoret, Nature (London) 504, 419 (2013).
[30] T. A. Palomaki, J. D. Teufel, R. W. Simmonds, and K. W.
Lehnert, Science 342, 710 (2013).
[31] K. C. Lee, M. R. Sprague, B. J. Sussman, J. Nunn, N. K.
Langford, X.-M. Jin, T. Champion, P. Michelberger, K. F. Reim,
D. England, D. Jaksch, and I. A. Walmsley, Science 334, 1253
(2011).
[32] J. D. Jost, J. P. Home, J. M. Amini, D. Hanneke, R. Ozeri,
C. Langer, J. J. Bollinger, D. Leibfried, and D. J. Wineland,
Nature (London) 459, 683 (2009).

063805-16

TWO-MODE SQUEEZED STATES IN CAVITY . . .

PHYSICAL REVIEW A 89, 063805 (2014)[33] D. F. Walls and G. J. Milburn, Quantum Optics (Springer, Berlin,2010).
[34] C. C. Gerry and P. L. Knight, Introductory Quantum Optics(Cambridge University Press, Cambridge, UK, 2008).
[35] M. J. Woolley, A. C. Doherty, G. J. Milburn, and K. C. Schwab,
Phys. Rev. A 78, 062303 (2008).
[36] K. Koga and N. Yamamoto, Phys. Rev. A 85, 022103 (2012).
[37] L. M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Phys. Rev.
Lett. 84, 2722 (2000).
[38] J. Fiurasek, Phys. Rev. A 66, 012304 (2002).
[39] A. Mari and D. Vitali, Phys. Rev. A 78, 062340 (2008).
[40] A. Mari and J. Eisert, Phys. Rev. Lett. 103, 213603 (2009).
[41] U. M. Titulaer and R. J. Glauber, Phys. Rev. 140, B676
(1965).

[42] Y.-D. Wang and A. A. Clerk, Phys. Rev. Lett. 108, 153603(2012).
[43] Sh. Barzanjeh, M. Abdi, G. J. Milburn, P. Tombesi, and D. Vitali,
Phys. Rev. Lett. 109, 130503 (2012).
[44] J. Bochmann, A. Vainsencher, D. D. Awschalom, and A. N.
Cleland, Nat. Phys. 9, 712 (2013).
[45] R. W. Andrews, R. W. Peterson, T. P. Purdy, K. Cicak, R. W.
Simmonds, C. A. Regal, and K. W. Lehnert, Nat. Phys. 10, 321(2014).
[46] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
[47] G. Adesso, A. Serafini, and F. Illuminati, Phys. Rev. A 70,
022318 (2004).
[48] P. Marian, T. A. Marian, and H. Scutaru, Phys. Rev. A 68, 062309(2003).

063805-17