PRL 110, 253601 (2013)

week ending
21 JUNE 2013

PHYSICAL REVIEW LETTERS

Reservoir-Engineered Entanglement in Optomechanical Systems
Ying-Dan Wang and Aashish A. Clerk
Department of Physics, McGill University, 3600 rue University, Montreal, Quebec H3A 2T8, Canada(Received 23 January 2013; published 18 June 2013)
We show how strong steady-state entanglement can be achieved in a three-mode optomechanical system(or other parametrically coupled bosonic system) by effectively laser cooling a delocalized Bogoliubov
mode. This approach allows one to surpass the bound on the maximum stationary intracavity entanglement
possible with a coherent two-mode squeezing interaction. In particular, we find that optimizing the relative
ratio of optomechanical couplings, rather than simply increasing their magnitudes, is essential forachieving strong entanglement. Unlike typical dissipative entanglement schemes, our results cannot be
described by treating the effects of the entangling reservoir via a Linblad master equation.
DOI: 10.1103/PhysRevLett.110.253601

PACS numbers: 42.50.Wk, 03.67.Bg, 42.50.Lc, 85.85.+jIntroduction.—The study of highly entangled quantum
states is of interest both for fundamental reasons and for a
myriad of applications in quantum information processing
and quantum communication. Of particular fundamental
interest is the possibility to entangle distinct macroscopic
objects, a task made difficult by the unavoidable decoherence and dissipation associated with such systems. Equallyinteresting would be the ability to entangle photons of very
different frequencies, e.g., microwave and optical photons.
A promising venue for the realization of both of thesekinds of entanglement is provided by quantum optomechanics, where macroscopic mechanical degrees of freedom can be controlled, measured, and coupled using themodes of an electromagnetic cavity. Recent milestones inthis field include the ability to cavity cool a mechanicalresonator to its ground state of motion [1,2] and the observation of many-photon strong coupling effects [3–6]. A
natural setting for entanglement generation is a three-mode
optomechanical system consisting of two ‘‘target’’ modes
to be entangled, which are each coupled to a third‘‘auxiliary’’ mode. One could either have two optical target
modes and a mechanical auxiliary mode, or vice versa; bothvariants have recently been achieved in experiment [7–9].
Several theoretical studies have described such schemes,
using the basic idea that the auxiliary mode mediates aneffective (coherent) two-mode squeezing interaction
between the two target modes (see, e.g., Refs. [10–13]).
However, such schemes typically yield at best a relatively
small amount of intramode entanglement (something which
we quantify more fully below).
In this Letter, we again consider generating steady-state
entanglement of two bosonic modes in a three-mode system; while we focus on an optomechanical realization, ourideas could also be realized using superconducting circuits
coupled via Josephson junctions [14,15] or other parametrically coupled three-bosonic-mode systems. Unlike previous works, we consider the possibility of entanglement
via reservoir engineering [16]: we wish to tailor the dissipative environment of the two target modes such that the0031-9007=13=110(25)=253601(5)dissipative dynamics relaxes the system into an entangledstate. Such dissipative entanglement has been discussed in
the context of atomic systems [17–21] and has even beenrealized experimentally [22].
The dissipative entanglement scheme we describe isrelated to optomechanical cavity-cooling schemes [23,24]
which have been used successfully to cool mechanicalresonators to the ground state. In our case, one is notcooling a simple mechanical mode to the ground state,
but rather a hybrid mode delocalized over both target
modes. In contrast to previous reservoir-engineering
approaches to entanglement generation, where the dynamics is reduced to a simple Markovian master equation for
the target degrees of freedom, our treatment is valid even in
the regime where a simple adiabatic elimination of theintermediate mode is not possible. As we show, this regime
turns out to be the most effective at generating entanglement. Our result shows that by optimizing the ratio of
optomechanical coupling strengths, rather than simply
increasing their magnitudes, this laser-cooling mechanism
can be used to yield large amounts of time-independent
intracavity entanglement. The amount of entanglement is
far greater than in previous studies and, in fact, far greaterthan the maximum possible entanglement allowed by acoherent parametric interaction. Note that reservoir engineering in optomechanics has previously been studied
theoretically, with the very different goal of generating
long-range coherence in arrays [25].
System and normal modes.—While our scheme applies
to a general bosonic three-mode system, we focus here on
an optomechanical system where two optical or microwave
cavity modes are coupled to a single mode of a mechanicalresonator (see Fig. 1); see the Supplemental Material [26]for a discussion of entangling two mechanical modes
coupled to a cavity mode. The Hamiltonian isX^ a^ y a^ i Þ þ H^ diss :H^ ¼ !M b^y b^ þ
ð!i a^ yi a^ i þ gi ðb^y þ bÞi

253601-1

i¼1;2

(1)
Ó 2013 American Physical Society

PRL 110, 253601 (2013)

PHYSICAL REVIEW LETTERS

cold mechanicalmode =
engineered reservoir

cavity 1

cavity 2

FIG. 1 (color online). Schematic of one realization of a threemode optomechanical system, where two cavity modes arecoupled to a single mode of a mechanical resonator. By driving
cavity 1 (2) at the red (blue) detuned mechanical sideband, adissipative entanglement mechanism is realized. To enhance thescheme, the mechanical resonator is cavity cooled and opticallydamped via a coupling to a third driven cavity mode, so thatits total damping rate  is greater than the damping rate  ofeach cavity.

a^ i is the annihilation operator for cavity i (frequency !i ,damping rate i ), b^ is the annihilation operator of themechanical mode (frequency !M , damping rate ), andgi are the optomechanical coupling strengths. H^ diss
describes the dissipation of each mode, as well as thedriving of the cavity modes. To achieve an entangling
interaction, cavity 1 (2) is driven at the red (blue) sideband
associated with the mechanical resonator: !d1 ¼ !1  !M
and !d2 ¼ !2 þ !M [12]. We work in an interaction
picture with respect to the cavity drives, and write a^ i ¼
a i þ d^i where a i is the classical cavity amplitude. We takeja 1;2 j  1, which allows us to linearize the optomechanical interaction in the usual way (i.e., drop interaction
terms not enhanced by the classical cavity amplitudes).
The linearized Hamiltonian in the rotating frame is thusH^ ¼ !M ðb^y b^ þ d^y1 d^1  d^y2 d^2 Þ þ H^ int þ H^ CR þ H^ diss with
^ þ G2 ðb^d^2 þ d^y b^y Þ;
H^ int ¼ G1 ðb^y d^1 þ d^y1 bÞ2(2)^ þ G2 ðb^y d^2 þ d^y bÞ:^
H^ CR ¼ G1 ðb^y d^y1 þ d^1 bÞ
2(3)Here Gi ¼gi a i (we take gi , a i >0 without loss of generality).
We further focus on the resolved-sideband regime !M 
1 , 2 , which suppresses the effects of the nonresonantinteractions in H^ CR . The remaining interaction H^ int in
Eq. (2) has the basic form suitable for entangling d^1 andd^2 : on a heuristic level, the parametric-amplifier interaction^ and then the beam-splitter(G2 term) first entangles d^2 and b,interaction (G1 term) swaps the b^ and d^1 states, thus yielding the desired entanglement.
Note that if one made the interactions in Eq. (2) nonresonant (e.g., by detuning the cavity drives from the
sideband resonances by ), one could adiabatically eliminate the mechanical mode, resulting in a two-mode squeezing interaction H^ TMS ’ ðd^1 d^2 þ H:c:Þ with   G1 G2 =
[27]. Such an interaction naturally leads to entanglement,
but the amount is severely limited by the requirement ofstability   1;2 =2. We quantify the entanglement using

week ending
21 JUNE 2013

the standard measure of the logarithmic negativity EN(see the Supplemental Material [26]). One finds that the
maximum stationary intracavity entanglement due to thetwo-mode squeezing coupling is (for 1 ¼ 2 and zerotemperature) EN ¼ lnð1 þ 2=Þ  ln2  0:7. Many suggested schemes for entanglement generation in optomechanical systems are limited by this stability requirement.
In contrast, the resonant case we consider allows for analternative dissipative entanglement mechanism capable ofmuch larger EN . We will focus attention on the regime
G2 < G1 where (for 1 ¼ 2 ) our linear system is always
stable (see the Supplemental Material [26]). Defining
the effective two-mode squeezing parameter r ¼
arctanhðG2 =G1 Þ, we introduce delocalized (canonical) cavity Bogoliubov mode operators:^ d^1 S^y ðrÞ;^ A ¼ d^1 coshr þ d^y2 sinhr  SðrÞ^ d^2 S^y ðrÞ:^ B ¼ d^y1 sinhr þ d^2 coshr  SðrÞ(4)^  exp½rd^1 d^2  H:c: is a two-mode squeezing
Here, SðrÞoperator. It thus follows that the joint vacuum of ^ A , ^ B is^
the two-mode squeezed state jri ¼ SðrÞj0;0i, where j0; 0i^^
is the vacuum of d1 , d2 . The entanglement of this state issimply EN ¼ 2r.
In terms of these new operators, H^ 0 ¼ !M ðb^y b^ þy^ A ^ A  ^ yB ^ B Þ and the optomechanical interactions in
Eqs. (2) and (3) take the simple form~ ^ y b^ þ H:c:;~ ^ y b^y þ H:c:; (5)H^ CR ¼ G
H^ int ¼ G
AA
qﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃﬃ~  G2  G2 . The mode ^ B completely decouwhere G
12
ples from the mechanics (it is a mechanically dark mode
[28,29]), while in the good-cavity limit of interest H^ CR canbe neglected, implying that the mode ^ A has a simple
beam-splitter interaction with the mechanics. H^ 0 þ H^ int
is trivially diagonalized,
resulting in hybridized modes
pﬃﬃﬃ^~ The existence^  ¼ ð^ A  bÞ=
2 with energies !M  G.
of three distinct eigenmodes (two hybrid, one dark) can be
useful to understand entanglement (in particular spectral
entanglement [13,30]) in the case where the mechanicalmode is driven by excessive thermal noise; we will discuss
this in a future work. We focus here on generating intracavity entanglement, which has the benefit of being insensitive to whether internal losses contribute to the damping
rate  of the cavities.
We now exploit the fact that Eq. (5) has exactly the form
used for standard cavity cooling [23,24]. Thus, if we cancouple the mechanical mode b^ to a cold reservoir, then the
beam-splitter coupling H^ int can be used to cool ^ A towards
vacuum, resulting in a stationary entangled state. A highfrequency, low-Q mechanical resonator would thus beideal. Alternatively, we will take the mechanical mode to
be coupled to a third cavity mode which is used to laser
cool its thermal occupancy towards the ground state byproviding a source of cold damping (see Fig. 1). In what

253601-2

PRL 110, 253601 (2013)

PHYSICAL REVIEW LETTERS

follows, we include the cooling cavity coupled to themechanical resonator in the definition of its effective
thermal bath; hence, the damping rate  includes the large
contribution of the optical damping. Amusingly, ourscheme is one of the few examples in optomechanics
where the enhanced mechanical damping rate resulting
from cavity cooling is actually highly beneficial.
Langevin equations and cavity cooling.—To describe thecooling potential of H^ int , we next use input-output theory
to derive the Heisenberg-Langevin equations for our linearized system. These take the standard form


d ^

pﬃﬃﬃﬃ
b ¼ i!M  b^  iðG1 d^1 þ G2 d^y2 Þ  b^in ;
dt
2


d ^
1 ^
pﬃﬃﬃﬃﬃﬃ
d1 ¼ i!M 
d1  iG1 b^  1 d^1;in ;
(6)
dt
2


d ^y
2 ^y
ﬃﬃﬃﬃﬃ
ﬃ
p
d ¼ i!M d þ iG2 b^  2 d^y2;in ;dt 2
2 2
where d^i;in , b^in describes operator-valued white noise driving the cavity and mechanical modes, and we have taken
the good-cavity limit !M   (allowing us to drop terms
due to H^ CR ). Equations (6) are readily solved to find thesteady-state occupancy and correlation of the Bogoliubov
modes. In the following analytic expressions, we take
1 ¼ 2 for simplicity and focus on the good-cavity limit(though Fig. 2 includes corrections due to H^ CR ).
Imagine first that the optomechanical interactions vanished, i.e., H^ int ¼ 0, and consider the behavior of ^ A , ^ B
(defined for a fixed r > 0). Even at zero temperature, ^ A
and ^ B will have a nonzero occupancy: the Bogoliubovtransformation of Eq. (4) implies that vacuum noise driving
the cavities acts as effective thermal noise for ^ A , ^ B .
Writing these intrinsic (H^ int ¼ 0) occupancies ash^ yj ^ j i0 ¼ n th;j we haven th;A=B ¼ n th;1=2 cosh2 r þ ðn th;2=1 þ 1Þsinh2 r:(7)Here, n th;1 (n th;2 ) represents the temperature of the thermalbath coupled to cavity 1 (2). As one increases the squeeze
parameter r, the effective heating of the ^ j modes becomes
exponentially large, implying the state of the system is far
from being an ideal two-mode squeezed vacuum state; thestate is not entangled.
Including now the effects of Hint [and taking r ¼
arctanhðG2 =G1 Þ], the dark-mode ^ B is unaffected, whereas
the occupancy of ^ A is modified to


opt

y ^
^
1þ
n
hA A i ¼
opt þ  þ  th;Aþ

opt n ;ðopt þ Þð þ Þ th;M(8)

where n th;M represents the temperature of the mechanicalbath (which includes the cooling cavity), and the effective ‘‘cold damping rate’’ of ^ A by the mechanics isweek ending
21 JUNE 2013

FIG. 2 (color online). Stationary intracavity entanglement(as quantified by log negativity EN , left scale) as a function of
the entangling interaction G2 . We take n th;1 ¼ n th;2 ¼ 0, a
mechanical frequency of !M ¼ 2 10 MHz, damping  ¼
2 0:8 MHz, and do not make the rotating wave approximation. The solid red thin curve corresponds to a fixed value of
G1 ¼ 2 2 MHz, and 1 ¼ 2 ¼ 2 50 kHz. One clearlysees a nonmonotonic dependence on G2 . The solid blue thickcurve and short-dashed blue curve instead correspond to tuningfor each G2 , such that the dissipative
G1 to the value Gopt1
entangling mechanism can be optimized. The value of Gopt
1 
G2 [cf. Eq. (10)] is indicated by the long-dashed brown curve(right scale). The solid blue (dashed blue) curve corresponds to1 ¼ 2 ¼ 50 kHz (1=2 =2 ¼ 45 kHz, 55 kHz) and n th;M ¼ 0(n th;M ¼ 0:3). The green dash-dotted line represents the maximum stationary entanglement achievable with a coherent twomode squeezing interaction EN ¼ ln2.

~ 2 =. This is the familiar equation for cavity
opt  4G
cooling in the good-cavity limit, where now the mechanics
plays the role of a cold reservoir. For n th;M ¼ 0 and weakcoupling (opt  ), the ^ A mode is cooled by a factor
=ð þ opt Þ. In the strong coupling limit, the coolingfactor saturates to a value =ð þ Þ.
Thus, while even vacuum noise tends to heat ^ A , ^ B toan exponentially large effective temperature, the optomechanical interaction of Eq. (5) can be used to cool ^ A .
Using the inequality of Duan et al. [31], one can show that
if one cools ^ A so thath^ yA ^ A i  sinh2 r;(9)then the two cavities must necessarily be entangled(see the Supplemental Material [26]). As the orthogonal
Bogoluibov mode ^ B is decoupled from the mechanics, it
is not cooled, making it impossible to achieve an ideal
two-mode squeezed vacuum state. Nonetheless, we findthat simply cooling ^ A is sufficient to generate a steadystate with significant entanglement (EN  2r  ln2 in thelarge r limit); this is despite the fact that the resulting state
has negligible overlap with a two-mode squeezed vacuum(see the Supplemental Material [26]).
To rigorously quantify the cavity-cavity entanglement,
we compute and discuss in what follows the log negativity

253601-3

EN , which is a function of the covariance matrix; details
are provided in the Supplemental Material [26].
Maximizing entanglement.—We now see that entanglement generation is more subtle than one might expect
given the simple form of Eq. (2). In particular, if G1 isfixed, the amount of stationary entanglement is a nonmonotonic function of the entangling-interaction strength
G2 (see Fig. 2). The dissipative entanglement mechanism
discussed here directly explains this behavior, as increasing G2 has two opposing effects: it not only increases
r and the delocalization of the Bogoliubov modes (enhancing entanglement), but also increases the effective temperature of these modes. This latter effect is due both to anincrease in the effective temperature of the cavity vacuum
noise [cf. Eq. (7)], and to a suppression of the cavity~ decreases withcooling effect (as the effective coupling G
increasing G2 ).
The maximum entanglement is achieved by carefullybalancing the opposing tendencies described above; without this optimization, the entanglement will remain small.
For fixed couplings G1 , G2 , one can optimize the entanglement as a function of the mechanical damping . The
maximum occurs at a nonzero dissipation strength, whichat zero mechanical temperature and in the good-cavity~ This value simply minilimit is simply given by  ¼ 2G.
^
mizes the occupancy of A , and corresponds to a simple
impedance matching condition (i.e., the rate with whichthe ^ A mode and mechanics exchange energy matches therate at which the mechanics and its bath exchange energy).
More relevant to experiment is to consider  and  fixed,
and optimize the entanglement over coupling strength.
Focusing on the most interesting regime where thecooperativity C2  G22 =ðÞ  1, and considering thegood-cavity limit and zero temperature (the mechanicalresonator is also cooled to vacuum by the third cavitymode), we find that for fixed G2 , the optimal G1 is given by
Gopt
1

rﬃﬃﬃﬃﬃﬃﬃ


2
1þ
;
G2 þ
8


~ opt
G
i:e:;


1
ln½2C2  ¼ 2r  2 ln2:2

For large C2 , the entanglement is almost that of a twomode squeezed vacuum (i.e., EN ¼ 2r). Alternatively, we
could hold the ratio = fixed and let C2 ! 1; we have(at zero temperature)ð=Þ2p
ﬃﬃﬃﬃﬃﬃ
þ
O
Eopt
:
(12)

ln
2
þ
N

C2
This is the strong-interaction limit, where the ^ A mode
hybridizes with the mechanical resonator. The maximal
cooling of ^ A is consequently set by the ratio =[cf. Eq. (8)]. The amount of entanglement here increases
monotonically from ln2 (the maximum possible with acoherent coupling) as this cooling factor is increased.
The behavior of the stationary entanglement versus
coupling strength is shown in Fig. 2, where we haveused parameters similar to those achieved in recent stateof-the-art experiments on microwave-circuit optomechanical systems [1,5]. We assume that a !M ¼ 10 MHz
mechanical resonator is first cavity cooled to near its
ground state, with a final damping rate of  ¼ 0:8 MHz(which is predominantly due to the cold optical damping
used for the cooling). By then optimally tuning the couplings to the target modes G1 , G2 to optimize the dissipative entanglement mechanism (while keeping them
& 2:2 MHz), one can obtain a relatively large EN  2:1.
This exceeds by an order-of-magnitude the intracavity
entanglement obtained in previous studies of the same
system [12], as well as the maximum of ln2 possible
with a coherent two-mode squeezing interaction. If thisentanglement was used for a teleportation experiment,
the maximum possible fidelity would be 0.89 [32,33];
this reduces the error by a factor of 3 compared to what
would be possible with EN ¼ ln2. Figure 2 also shows that10
C2 2 1=4
:
22
(10)

Note that for large C2 , this optimal value can easily corre~ > , . Thus, in thisspond to a strong interaction G
optimal regime, the effects of our ‘‘engineered reservoir’’(cold mechanical resonator) on the target cavity modes
cannot be described by a Markovian dissipator in a masterequation; this is in stark contrast to standard dissipationby-entanglement schemes.
For the optimal value of G1 above, the entanglementtakes simple forms in two relevant limits. If we hold C2
fixed while taking the limit = ! 1, we have (at zerotemperature)
Eopt
N

week ending
21 JUNE 2013

PHYSICAL REVIEW LETTERS

Entanglement EN

PRL 110, 253601 (2013)

(11)

8
6

γ /κ = 10004γ /κ = 100

2γ /κ = 10

0
1

100

104

106

108

Cooperativity C2

FIG. 3 (color online). Intracavity stationary entanglement(quantified by EN ) versus cooperativity C2 , where we use an
optimized choice for G1 [as given by Eq. (10)], and have taken
1 ¼ 2 and zero temperature. The solid lines correspond todifferent choices of the damping ratio = as indicated; increasing = increases the amount that one can cool the delocalized^ A mode, and hence enhances entanglement. For large C2 , thesecurves asymptote to the value in Eq. (12). The dashed blue line isthe asymptotic expression of Eq. (11). The green dashed-dottedline indicates EN ¼ ln2 (the maximum EN possible with a twomode squeezing interaction).

253601-4

PRL 110, 253601 (2013)

PHYSICAL REVIEW LETTERS

large values of EN are possible even when n th;M  0 and1  2 .
In Fig. 3 we show how the stationary entanglement
grows to dramatically large values with C2 for an optimized choice of G1 . While the parameters needed for such
EN may be out of reach in current-generation optomechanics experiments, they may be more feasible by implementing a superconducting circuit realization of ourscheme [14,15].
Conclusions.—We have presented a general method forthe dissipative generation of entanglement in a three-mode
optomechanical system. The entanglement generated here
could be verified by measuring the covariance matrix ofthe two target cavity modes using homodyne techniques(see, e.g., Ref. [11]). Alternatively, one could directly
use the cavity output spectra at resonance to measure theoccupancy of the ^ A mode; verifying that it violates the
Duan inequality of Eq. (9) would also confirm the generation of entanglement (see the Supplemental Material [26]).
We thank S. Chesi and L. Tian for useful conversations.
This work was supported by the DARPA ORCHID program under a grant from the AFOSR.

[1] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman,
K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, andR. W. Simmonds, Nature (London) 475, 359 (2011).
[2] J. Chan, T. P. M. Alegre, A. H. Safavi-Naeini, J. T. Hill, A.
Krause, S. Gröblacher, M. Aspelmeyer, and O. Painter,
Nature (London) 478, 89 (2011).
[3] S. Weis, R. Riviere, S. Deleglise, E. Gavartin, O. Arcizet,
A. Schliesser, and T. Kippenberg, Science 330, 1520
(2010).
[4] A. H. Safavi-Naeini, T. P. M. Alegre, J. Chan, M.
Eichenfield, M. Winger, Q. Lin, J. T. Hill, D. E. Chang,and O. Painter, Nature (London) 472, 69 (2011).
[5] J. D. Teufel, D. Li, M. S. Allman, K. Cicak, A. J. Sirois,
J. D. Whittaker, and R. W. Simmonds, Nature (London)
471, 204 (2011).
[6] E. Verhagen, S. Deleglise, S. Weis, A. Schliesser, and T.
Kippenberg, Nature (London) 482, 63 (2012).
[7] J. T. Hill, A. H. Safavi-Naeini, J. Chan, and O. Painter,Nat. Commun. 3, 1196 (2012).
[8] C. Dong, V. Fiore, M. C. Kuzyk, and H. Wang, Science
338, 1609 (2012).
[9] F. Massel, S. U. Cho, J.-M. Pirkkalainen, P. Hakonen, T. T.
Heikkilä, and M. A. Sillanpaa, Nat. Commun. 3, 987(2012).
[10] S. Mancini, V. Giovannetti, D. Vitali, and P. Tombesi,
Phys. Rev. Lett. 88, 120401 (2002).

week ending
21 JUNE 2013[11] M. Paternostro, D. Vitali, S. Gigan, M. S. Kim, C.
Brukner, J. Eisert, and M. Aspelmeyer, Phys. Rev. Lett.
99, 250401 (2007).
[12] S. Barzanjeh, D. Vitali, P. Tombesi, and G. J. Milburn,Phys. Rev. A 84, 042342 (2011).
[13] S. Barzanjeh, M. Abdi, G. J. Milburn, P. Tombesi, and D.
Vitali, Phys. Rev. Lett. 109, 130503 (2012).
[14] N. Bergeal, R. Vijay, V. E. Manucharyan, I. Siddiqi, R. J.
Schoelkopf, S. M. Girvin, and M. Devoret, Nat. Phys. 6,296 (2010).
[15] B. Peropadre, D. Zueco, F. Wulschner, F. Deppe, A. Marx,
R. Gross, and J. J. Garcia-Ripoll, Phys. Rev. B 87, 134504(2013).
[16] J. F. Poyatos, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 77,
4728 (1996).
[17] M. B. Plenio and S. F. Huelga, Phys. Rev. Lett. 88, 197901(2002).
[18] B. Kraus and J. I. Cirac, Phys. Rev. Lett. 92, 013602(2004).
[19] A. S. Parkins, E. Solano, and J. I. Cirac, Phys. Rev. Lett.
96, 053602 (2006).
[20] B. Kraus, H. P. Büchler, S. Diehl, A. Kantian,A. Micheli, and P. Zoller, Phys. Rev. A 78, 042307(2008).
[21] 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).
[22] C. A. Muschik, E. S. Polzik, and J. I. Cirac, Phys. Rev. A
83, 052312 (2011).
[23] F. Marquardt, J. P. Chen, A. A. Clerk, and S. M. Girvin,
Phys. Rev. Lett. 99, 093902 (2007).
[24] I. Wilson-Rae, N. Nooshi, W. Zwerger, and T. J.
Kippenberg, Phys. Rev. Lett. 99, 093901 (2007).
[25] A. Tomadin, S. Diehl, M. D. Lukin, P. Rabl, and P. Zoller,
Phys. Rev. A 86, 033821 (2012).
[26] See Supplemental Material at http://link.aps.org/
supplemental/10.1103/PhysRevLett.110.253601
fordetails on the calculation, a discussion of the two
mechanical resonator system, and details on the cavity
output spectrum.
[27] M. Schmidt, M. Ludwig, and F. Marquardt, New J. Phys.
14, 125005 (2012).
[28] Y.-D. Wang and A. A. Clerk, Phys. Rev. Lett. 108, 153603(2012).
[29] L. Tian, Phys. Rev. Lett. 108, 153604 (2012).
[30] C. Wipf, T. Corbitt, Y. Chen, and N. Mavalvala, New J.
Phys. 10, 095017 (2008).
[31] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Phys.
Rev. Lett. 84, 2722 (2000).
[32] G. Adesso and F. Illuminati, Phys. Rev. Lett. 95, 150503(2005).
[33] A. Mari and D. Vitali, Phys. Rev. A 78, 062340(2008).

253601-5

Supplemental material for “Optomechanical entanglement via reservoir engineering”
Ying-Dan Wang and Aashish A. Clerk
Department of Physics, McGill University, 3600 rue University, Montreal, QC Canada H3A 2T8

SYSTEM STABILITY CONDITION

The linearized Heisenberg-Langevin equations for our three mode system are given in Eqs. (6) of the main text. If
we drop the noise terms, these equations take the form ~v = M · ~v where ~v = (b̂, dˆ1 , dˆ†2 ) and M is a 3 × 3 matrix. For
the system to be stable, we require that the eigenvalues of M all have a negative real part. This requirement leads
to the well-known Routh-Hurwitz stability conditions. In the simple case κ1 = κ2 , the stability conditions reduce tothe following necessary and sufficient condition:
G̃2 > −

κγ
4

(S1)

Thus, if κ1 = κ2 and G2 ≤ G1 , the system is always stable (regardless of the magnitude of the cavity and mechanicaldamping rates).
In contrast, if κ1 6= κ2 , then the stability conditions are modified such that the system can be unstable at G̃ = 0.
While the general form of the conditions is somewhat unwieldy, in the interesting large cooperativity limit (i.e. Ci ≡G2i /(κi γ)  0), they reduce to:κ22 − κ21
G̃2 > C̄γ max κ1 − κ2 ,(S2)2γ + κ1 + κ2
with C̄ ≡ (G21 + G22 )/(γ(κ1 + κ2 )). Note that in the γ  κ1 , κ2 limit of interest, this condition is very sensitive to thesign of κ1 − κ2 . On a heuristic level, damping asymmetry leads to negative damping terms in the equations of motion
for the modes β̂A , β̂B . For κ1 < κ2 (and G̃ > 0), it is β̂A which experiences negative damping; this is overwhelmedby the positive cold damping Γopt provided by the interaction with the mechanics (c.f. Eq. (8) in the main text), andhence there is no possibility of instability. In contrast, if κ1 > κ2 , it is the dark mode β̂B which experiences negative
damping; as there is nothing to offset this, the system can become unstable.

DEFINITION OF THE LOGARITHMIC NEGATIVITY

In this paper, unless specified, we use the logarithmic negativity EN to quantify the degree of entanglement; thisquantity is a rigorous entanglement monotone, and is zero for separable states. For two-mode Gaussian states of sortrealized by the two cavity modes dˆ1 , dˆ2 in our system, it can be calculated using the expression [1]EN = max[0, − ln 2η − ]

(S3)

with
1
η =√
2
−

q
p
Σ − Σ2 − 4 det V

(S4)

and
Σ = det B + det B0 − 2 det C.

(S5)
E
Here V is the 4 × 4 covariance matrix of the two modes of interest, defined via Vjj 0 = 12 ∆ξˆj ∆ξˆj 0 + ∆ξˆj 0 ∆ξˆj , with

 √
h
i

 √
ˆ
∆ξˆj = ξˆj − hξˆj i and ξ~ = {x̂1 , p̂1 , x̂2 , p̂2 }. Here x̂i = dˆi + dˆ†i / 2 and p̂i = −i dˆi − dˆ†i / 2 and dˆi , dˆ†j = δij . The
D

matrix B, B0 and C are 2 × 2 matrices related to the covariance matrix V as


B C
V=
.
CT B0

(S6)

2
The evaluation of entanglement needs the full information of the covariance matrix. In our main text, the relevant
system consists of two cavity modes dˆ1 and dˆ2 , or equivalently β̂A and β̂B . Thus we also the correlation of β̂A andβ̂B , beside their occupancy (Eqs. (7) and (8) in the main text). The only non-zero correlator (again working in the
good-cavity limit ωM  κ1 , κ2 where the counter-rotating terms ĤCR have negligible effect) isDEβ̂A β̂B =κ (κ + γ)(n̄th,1 + n̄th,2 + 1) sinh 2r.
Γopt γ + 2κ (κ + γ)

(S7)

ENTANGLEMENT OF MECHANICAL RESONATORS MEDIATED BY A CAVITY MODE

While the main text focuses on a 3-mode optomechanical system having two cavity modes coupled to a single
mechanical mode, our entanglement-by-dissipation scheme is very general, and can be applied in principle to any
set of three parametrically-coupled bosonic modes. In particular, it could be realized in a 3-mode optomechanicalsystem where two mechanical modes are coupled to a single cavity mode. We analyze this case below. Note thatafter this work was submitted, a paper by Tan et al. appeared which also analyzes dissipative entanglement of two
mechanical resonators coupled to a cavity [2]; unlike our analysis before, this work does not explicitly consider theeffect of non-RWA terms, terms we find to be especially problematic.
The starting Hamiltonian is now



X  †
Ĥ = ωcav â† â +ωi b̂i b̂i + gi b̂†i + b̂i â† â + Ĥdiss .
(S8)
i=1,2

Here b̂i is the annihilation operator of the i-th mechanical resonator (frequency ωi , damping rate γi ), â is the annihilation operator of the auxiliary cavity mode (frequency ωcav , damping rate κ), and gi are the optomechanical couplingstrengths. Ĥdiss describes the dissipation of each mode, as well as the driving of the cavity mode.
To achieve an entangling interaction, we will drive the cavity mode at two frequencies ωd1 = ωcav − ω1 andωd2 = ωcav + ω2 (i.e. the red(blue) sideband associated with mechanical resonator 1 (2)). We can then write the cavityˆ −iωcav t , where the classical cavity amplitude ā(t) = ā1 eiω1 t + ā2 e−iω2 t e−iωcav tannihilation operator as â = ā(t) + deis obtained from the classical equations of motion, with the sideband amplitudes ā1,2 determined independently by
the two cavity drives. We take |ā1,2 |  1 which allows us to linearize the optomechanical interaction in the usual way(i.e. we drop interaction terms not enhanced by the classical cavity amplitude). Working in an interaction picture
ˆ the linearized Hamiltonian in the displaced frame can be written as
with respect to Ĥ0 = ω1 b̂†1 b̂1 + ω2 b̂†2 b̂2 + ωcav dˆ† d,
Ĥ = Ĥint + ĤCR + Ĥdiss withĤint = G1 dˆ† b̂1 + dˆb̂†1 + G2 dˆb̂2 + dˆ† b̂†2 .
(S9)Here Gi = gi āi (we take gi , āi to be positive without loss of generality) and ĤCR describes counter-rotating termswith an explicit time-dependence. Eq. (S9) takes the same form as Eq. (2) in the main text and can be used forentangling b̂1 and b̂2 by reservoir engineering.
As discussed in the main text, the dissipative entanglement mechanism is optimal when the mode acting as theengineered reservoir is both cold and has a large damping rate compared to the target modes, as this allows optimalcooling of the Bogoliubov mode β̂A . This occurs naturally here, as the reservoir mode is a cavity mode, which given
its higher frequency (compared to mechanical modes) will be naturally closer to the ground state, and which naturally
has a large damping rate (i.e. its damping κ will be much larger than the damping rates γ of the target mechanicalmodes, by a factor ∼ 104 to 105 in typical optomechanical experiments). An additional benefit here is that since themechanical damping rate γ is small, the cooperativity of the system can reach a much larger value (in the strong
coupling limit, C can reach 106 ). Thus the entanglement can easily reach EN ∼ 7 (as shown in Fig. 3), corresponding
to teleportation fidelity of 99.9%. In contrast, the two-cavities-plus-mechanics scheme in the main text requires one
to optically damp the mechanical damping γ to a large value.
While the above features are extremely attractive, the downside of this scheme comes when one considers the nonresonant, counter-rotating terms. For the setup discussed in the main text (two cavities, one mechanical resonator), thecounter-rotating terms only lead to a very small heating of the coupled β̂A mode, much like standard optomechanicalcavity cooling (c.f. Eq. (5) in the main text). In contrast, the situation here (two mechanical resonators, one cavity) is
more complex. In general there are two sets of non-resonant, counter-rotating terms. The first are “diagonal terms”,

3
where, e.g., the cavity drive at ωcav − ω1 induces Stokes scattering with mechanical resonator 1. In our rotating frame,
these take the form:

 

A
HCR
= ā∗1 g1 dˆb̂1 e−2iω1 t + h.c. + ā∗2 g2 dˆb̂†2 e2iω2 t + h.c. .
(S10)We also have non-resonant “off-diagonal” terms where now, e.g., the cavity drive at ωcav − ω1 induces both Stokes
and anti-Stokes scattering involving mechanical resonator 2. These take the form:


B
HCR
= ā∗1 g2 ei(ω2 −ω1 )t dˆb̂†2 + e−i(ω1 +ω2 )t dˆb̂2+ā∗2 g1 ei(ω2 −ω1 )t dˆb̂1 + ei(ω1 +ω2 )t dˆb̂†1 + h.c.
(S11)
To understand the effects of these terms, we re-express them in terms of the Bogoluibov modes β̂A , β̂B using Eq.
(4) in the main text (with the proviso that in these equations, dˆi should be replaced with b̂i ). We thus have:A
HCR= G̃dˆβ̂A cosh2 re−2iω1 t − sinh2 re2iω2 t+G̃dˆβ̂ † sinh r cosh r e2iω2 t − e−2iω1 t + h.c.
(S12)
B

and





†
†
B
HCR
= ā∗1 g2 ei(ω2 −ω1 )t dˆ −β̂A sinh r + β̂Bcosh r + e−i(ω1 +ω2 )t dˆ −β̂Asinh r + β̂B cosh r





†
†
+ā∗2 g1 ei(ω2 −ω1 )t d β̂A cosh r − β̂Bsinh r + ei(ω1 +ω2 )t dˆ β̂Acosh r − β̂B sinh r
+h.c.

(S13)

The crucial difference here is that unlike the the case discussed in the main text (2 cavities, 1 mechanical resonator),
the non-resonant terms when written in the β̂A , β̂B basis enter with factors of Gi cosh r, Gi sinh r (i = 1, 2) (recall
that G1 = G̃ cosh r, G2 = G̃ sinh r). Their effects will thus be exponentially larger than in the case considered in themain text, as in that case, the counter-rotating terms only enter with a small prefactor ∝ G̃ (compare against Eq. (5)in the main text). As a result, if one wants a large squeeze parameter r, one will need to be extremely far into the
good-cavity limit to suppress the deleterious effects of the non-resonant terms.
To estimate how far one needs to be in the good cavity limit in this 2-mechanical resonator setup, note that the
B
A
that contain dˆβ̂A , dˆβ̂B (and their Hermitian conjugates) will cause negative damping of the β̂A , β̂B, ĤCR
terms in ĤCR
modes. If this negative damping becomes too large (i.e. larger than the intrinsic damping γ of these modes), this cancause instability. Insisting that the induced negative damping is less than γ in the large cooperativity, large-r limitof interest (i.e. where the entanglement can be large) leads to the requirement:
G21 e2r  κ 2
.γ
κ
ω̃

(S14)

where ω̃ = min [2ω1 , 2ω2 , |ω1 − ω2 |]. Assuming the optomechanical couplings G1 has been chosen to optimize EN
in the absence of counter-rotating terms and with G2 fixed (as per Eq. (10) in the main text), we can re-write thiscondition in terms of the cooperativity C2 = G22 /(κγ)
κ
.
ω̃



1
C2

3/4

Thus, obtaining significant entanglement (which necessarily requires a large C2 ) also requires one to be extremely
deep into the good-cavity limit, i.e. ωM  κ. For example, in order to achieve EN ∼ 5, the ratio κ/ω̃ has to besmaller than 10−2 , a value that would be extremely challenging for current optomechanics experiments. Microwave
cavity experiments have achieved a ratio ωM /κ ∼ 50 [3]; with such numbers, an EN & 1.5 could be possible. This isalready more than a factor of two greater than the maximal entanglement EN = ln 2 possible with a direct two-mode
squeezing interaction.

4
THE BOUND OF THE BOGOLIUBOV MODE OCCUPANCY BY THE GENERALIZED DUAN
INEQUALITY

Using the definition of Bogoliubov mode Eq. (4) in the main text, one gets
Dn
oE
†
β̂A
, β̂A
≡ 1 + 2ñA∗
= V11 cosh2 r + V22 sinh2 r + (C12 + C12) sinh r cosh r,
Dn
oE
Dn
oE
Dn
oE
with V11 ≡dˆ†1 , dˆ1 , V22 ≡dˆ†2 , dˆ2
and C12 ≡dˆ1 , dˆ2where {· · · } denotes anti-commutator.

According to generalized Duan inequality [4], any separable (non-entangled state) will necessarily satisfy*2 + *
2 +
1
1
1
D≡
ax̂1 + x̂2+
ap̂1 − p̂2≥ 2 + a2 ,
a
a
a

(S15)

(S16)


 √

 √
where x̂i = dˆi + dˆ†i / 2 and p̂i = −i dˆi − dˆ†i / 2 and a is an arbitrary nonzero real number.
If we set a2 = coth r, we can easily verify that1 + 2ñA =Dsinh 2r.
2(S17)

For this choice of a, using the Duan bound on D thus yields that for any separable state, the occupancy of the β̂Amode must satisfy:ñA ≥ sinh2 r.

(S18)

Note from Eq. (7) in the main text that at zero temperature at zero optomechanical coupling, ñA is exactly sinh2 r,
indicating no violation of the Duan bound (consistent with no entanglement in this case). However, any amount ofcooling (as per Eq. (8)) in the main text will then lower the value of ñA , leading to a violation of the Duan boundand hence entanglement of the two target cavity modes.
COOLING ONE BOGOLIUBOV MODE VS. COOLING TWO BOGOLIUBOV MODES

As explained after Eq. (4) of the main text, if both Bogoliubov modes (β̂A and β̂B ) are cooled to vacuum, a twomode squeezed vacuum of dˆ1 and dˆ2 can be achieved. Two-mode squeezed vacuum is a highly-entangled state withentanglement (evaluated by log negativity) EN = 2r. This is the central idea lying in previous works of dissipation
generated entanglement (see. e.g. Ref [21-22] in the main text). However, in our setup, one of the Bogoliubov
modes β̂B is decoupled from the intermediate mode (cf. Eq. (5)) which makes it impossible to be cooled by lasercooling mechanism. Thus the ideal steady-state to be achieved with our protocol, ρ̂0 , is characterized by the following
correlations
D
E
D
E
D
E
†
†
β̂A
β̂A = 0, β̂Bβ̂B = sinh2 r, β̂B β̂A = 0(S19)
where the occupany of the β̂B mode increases exponentially with the squeezing parameterr. The corresponding√√
covariance matrix of the state ρ̂0 in the basis of {x̂1 , p̂1, x̂2 , p̂2 } (where x̂i = (dˆi + dˆ†i )/ 2 and p̂i = −i(dˆi − dˆ†i )/ 2) is
3 + cosh 4r
0−8 cosh3 r sinh r
0
3

1
0
3 + cosh 4r
0
8 cosh r sinh r.
V(ρ̂0 ) = 
(S20)
3


0
4 cosh 2r + cosh 4r − 1
0
2 −8 cosh r sinh r
0
8 cosh3 r sinh r
0
4 cosh 2r + cosh 4r − 1hi
The overlap between this state and an arbitrary two mode squeezed vacuum |r̃, θi = exp r̃eiθ dˆ1 dˆ2 − h.c. |0, 0i (where|0, 0i is the vacuum state of the modes dˆ1 and dˆ2 ) can be evaluated as
F ≡ Tr (ρ̂0 ρ̂TMSV ) =41
2 2 + cosh (r − r̃) + cosh (r + r̃) − 2 cos θ sinh 2r sinh 2r̃ ,cosh r(S21)

5
where ρ̂TMSV = |r̃, θihr̃, θ| is the density matrix of the 2-mode squeezed vacuum. For large squeezing r  1, theoverlap between the two states is negligible F ∼ exp (−6r) → 0, as mentioned after Eq. (9) in the main text.
Of course the cooling of β̂B mode can also be achieved by introducing a second intermediate mode and two-tone
drives on the cavity modes. However, we find this extra complication is not necessary for our purpose. In fact,
although our scheme is insufficient to generate a 2-mode squeezed vacuum, surprisingly it is sufficient to generate a
steady state with significant entanglement EN = 2r − ln 2, which is almost as good as a 2-mode squeezed vacuum inthe large r limit. To understand the physical meaning of the state ρ̂0 , we can compare Eq. (S20) with the covariancematrix of a 2-mode squeezed thermal state [5]ρTMST = S̃ (r) ρ1,th ⊗ ρ2,th S̃ † (r)

(S22)

with
ρi,th =∞
X(n̄ei )ni

|ni i hni |e ni +1
ni =1 (1 + n̄i )(S23)

the thermal state of mode i (i = 1, 2) with average thermal occupancy n̄ei . Interestingly, we find that, Eq. (S20) canbe reproduced by setting n̄e1 = 0 and n̄e2 = sinh2 r. This shows the state obtained in our work is a two-mode squeezed
thermal state with the two modes at different effective temperatures.

MEASURING ENTANGLEMENT VIA THE β̂A OCCUPANCY

As shown in the previous section, if one can verify that the occupancy of the β̂A mode has been reduced below thevalue sinh2 r (where r = arctanh G2 /G1 ), then one has violated the Duan inequality, indicating that cavities 1 and 2
are necessarily entangled. We now show that this mode occupancy can be directly obtained from the output spectraof the two cavities. Consider the simplest case where all dissipative baths are at zero temperature, and where we are
deep in the good-cavity limit; for concreteness, also focus on the output spectrum of cavity 1. The solution to the
Heisenberg-Langevin equation of motion for dˆ1 may be written:


 √

√
dˆ1 [ω] = χ1 [ω] −iG1 χM [ω] −iG̃β̂A [ω] − γ b̂in [ω] − κdˆ1,in [ω]where the susceptibilities are:1−i(ω − ωM ) + κ/21χM [ω] =−i(ω − ωM ) + γ/2
χ1 [ω] =

(S24)
(S25)

Further, the output field from cavity 1 dˆ1,out is given by the standard input-output relation√dˆ1,out = dˆ1,in + κdˆ1 ,and its power spectral density defined asZ
S1 [ω] =dteiωt hdˆ†1,out (t)dˆ1,out (0)i.

Note that the above equations (like the equations in the main text) are in the rotating frame set by the cavity drive
frequencies; hence, ω = ωM corresponds to the cavity resonance.
We see immediately that the only contribution to S1 [ω] will be from the non-zero occupancy of β̂A . We focus for
simplicity on the regime where the effective coupling G̃ is less than γ, a regime where large entanglement is possible
but the β̂A mode does not hybridize with the mechanical resonator; the entanglement in this regime is described by
Eq. (11) of the main text. In this regime, the total damping rate of the β̂A mode will be κ + Γopt , where the effective
cold damping rate due to the mechanics is Γopt = 4G̃2 /γ (c.f. Eq. (8) in the main text). A straightforward calculation
then shows that
S1 [ωM ] = 16C1

4C̃
†
hβ̂A
β̂A i
1 + 4C̃(S26)

6
where the co-operativities C1 and C̃ are defined as C1 = G21 /(κγ), C̃ = G̃2 /(κγ). We thus see that the output
spectrum of cavity 1 at resonance (i.e. the number of photons leaving the cavity at its resonance frequency) gives a
direct measure of the occupancy of the coupled Bogoluibov mode β̂A .

[1] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
[2] H. Tan, G. Li, and P. Meystre, Phys. Rev. A 87, 033829 (2013).
[3] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman, K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, andR. W. Simmonds, Nature 475, 359 (2011).
[4] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 84, 2722 (2000).
[5] P. Marian, T. A. Marian, and H. Scutaru, Phys. Rev. A 68, 062309 (2003).