PHYSICAL REVIEW A 91, 013807 (2015)

Bipartite and tripartite output entanglement in three-mode optomechanical systems
Ying-Dan Wang,1,2 Stefano Chesi,3 and Aashish A. Clerk2
1

State Key Laboratory of Theoretical Physics, Institute of Theoretical Physics, Chinese Academy of Sciences, Beijing 100190, China
2
Department of Physics, McGill University, 3600 rue University, Montreal, Quebec, Canada H3A 2T8
3
Beijing Computational Science Research Center, Beijing 100084, China(Received 5 July 2014; revised manuscript received 8 December 2014; published 5 January 2015)
We provide analytic insight into the generation of stationary itinerant photon entanglement in a three-mode
optomechanical system. We identify the parameter regime of maximal entanglement and show that strong
entanglement is possible even for weak many-photon optomechanical couplings. We also show that strong
tripartite entanglement is generated between the photonic and phononic output fields; unlike the bipartite photonphoton entanglement, this tripartite entanglement diverges as one approaches the boundary of system stability.
DOI: 10.1103/PhysRevA.91.013807

PACS number(s): 42.50.Wk, 07.10.Cm, 42.50.Ex

I. INTRODUCTION

Entanglement is one of the most fascinating and potentially
useful aspects of quantum systems. Of particular interestis the generation of entangled itinerant quanta (which canbe easily spatially separated) and of true multipartite entanglement (involving irreducible correlations between threeor more subsystems). These goals have been the subject ofconsiderable theoretical and experimental work in a variety ofsystems spanning quantum optics setups [1,2], cold atoms [3],superconducting circuits [4–6], and spin qubits [7]. Optomechanical systems [8], where mechanical motion interacts withelectromagnetic fields, could be another powerful platform
to realize these goals. A key advantage here is the potential
to use mechanical motion to entangle disparate subsystems(e.g., microwave and optical photons). A number of schemes
to generate entangled photons in optomechanics have beenstudied theoretically [9–15]. Recent experiments have also
demonstrated mechanically mediated entanglement between
two microwave pulses [16].
Here, we analyze theoretically both itinerant and multipartite entanglement in a three-mode optomechanical system
in which two cavities are coupled to a single mode of amechanical resonator. Such a setup has been realized in several
recent experiments [17–19]. Previous theory work, largely
numerical, examined bipartite output entanglement in thissystem [11,13–15] but mostly focused on transient [15] orexperimentally -challenging strong-coupling regimes [13,14].
In contrast, we focus here on identifying the optimal conditions
to generate stationary output entanglement with weak manyphoton optomechanical couplings. We find that maximum
entanglement is achieved through a simple matching conditionof optomechanical cooperativities. Surprisingly, this operation
point coincides with the least favorable regime for thegeneration of intracavity entanglement. We provide a complete
yet simple analytic understanding of the physics, which also
allows us to illustrate the trade-off between large entanglement
and thermal resilience and to uncover a sensitive dependence ofentanglement on the time delay between emitted wave packets.
We also address the generation of tripartite entanglementin such a hybrid three-mode system, considering correlations
between both output photons and phonons. While usuallyignored, the mechanical output field could be accessed
experimentally, using, for example, optomechanical crystal
1050-2947/2015/91(1)/013807(10)geometries with phonon waveguides [20] or by having themechanical dissipation be dominated by a third auxiliary
cavity used for cooling. We find that true hybrid tripartite
entanglement is indeed created: the output state corresponds to
a “twice-squeezed vacuum,” involving two two-mode squeezing operations. We also quantify this entanglement using theGaussian Rényi-2 measure [21]. Besides being of fundamental
interest, such tripartite entangled states have application toa variety of quantum information processing tasks such as
teleportation, telecloning, and dense coding [22,23]. The
present setup is especially suited to the continuous generationof nonlocal multipartite entanglement, as the phonons andphotons from the two cavities are all emitted into spatiallyseparated outputs. Note that Genes et al. [10] also studied
tripartite entanglement in an optomechanical system, although
in a setting where the entangled subsystems were not all
spatially separated or itinerant.
It also worth noting that, while our emphasis here is on
optomechanics, the physical features uncovered here apply
generally to other bosonic three-mode systems (e.g., realizedwith superconducting circuits [24,25]).
Our paper is organized as follows: in Sec. II, we describe
the physical system and its stability conditions. In Sec. III we
obtain the optimal condition for bipartite output entanglement
based on an analytical solution. In Sec. IV, we analyze
in detail various imperfections such as thermal noise, nonrotating-wave-approximation corrections, internal losses, andfinite bandwidth. In Sec. V, we describe how to achieve
nonlocal tripartite output, which we quantify with the “residual
Gaussian-Rényi-2” (GR2) entanglement measure. Section VI
summarizes our results.
II. PHYSICAL SYSTEM

We consider a three-mode optomechanical system in which
two cavities are coupled to a single mode of a mechanicalresonator (see Fig. 1).
The standard optomechanical interacting Hamiltonian††Ĥ = ωm b̂† b̂ +[ωi âi âi + gi (b̂† + b̂)âi âi ]
(1)
i=1,2

governs the system’s dynamics, where âi is the annihilation
operator for cavity i (frequency ωi ), b̂ is the annihilation
operator of the mechanical mode (frequency ωm ), and gi is013807-1©2015 American Physical SocietyYING-DAN WANG, STEFANO CHESI, AND AASHISH A. CLERK

FIG. 1. (Color online) System schematics. Two driven cavities(cavities 1 and 2) interact with a common mechanical resonator (modeb̂). This can generate entanglement in the optical outputs. An auxiliary
third cavity (cavity a) can be used to cavity cool the mechanics andto make the mechanical output mode accessible.

the optomechanical coupling strengths. In order to generate
steady-state entanglement, cavity 1 (2) is driven at the red(blue) sideband associated with the mechanical resonator:
ωd1 = ω1 − ωm and ωd2 = ω2 + ωm . By working in an interaction picture with respect to the free Hamiltonian andfollowing the usual linearization procedure [8], we have
Ĥ R = Ĥint + ĤCR (t), withĤint = (G1 b̂† d̂1 + G2 b̂d̂2 ) + H.c. = G̃(b̂† β̂A + H.c.)(2)†

and ĤCR (t) = G̃(b̂† β̂A e2iωm t + H.c.). Here d̂i = âi − āi , withāi being the classical cavity amplitude. Gi = gi āi is the
dressed coupling(we take gi ,āi > 0 without loss of gener√†
2
ality), G̃ = G1 − G22 , and β̂A = d̂1 cosh r + d̂2 sinh r [r =tanh−1 (G2 /G1 )] is a Bogoliubov mode.
We first focus on the resolved-sideband regime ωm 
κ1 ,κ2 and make the rotating-wave approximation (RWA) byneglecting ĤCR . The non-RWA corrections are treated inSec. IV C. The combined swapping and entangling interactionsin Ĥint lead to a net entangling interaction between the two⎛

χ1

PHYSICAL REVIEW A 91, 013807 (2015)

cavity modes; as discussed in [26], this interaction has afundamentally dissipative nature. It is useful to diagonalizeĤint in terms of three normal modes [14,26]: one “me†chanically dark” Bogoliubov mode β̂B = d̂1 sinh r + d̂2 cosh r(which is robust to mechanical thermal noise) and two coupledeigenmodes involving both the mechanics and cavities.
The dissipation of each mode is treated via standard inputoutput theory [27] and is given explicitly by the following
Langevin equations:dγ√†â = − â − i(G1 d̂1 + G2 d2 ) − γ â in ,
dt
2
√
κ1
d
d̂1 = − d̂1 − iG1 â − κ1 d̂1in ,(3)dt
2√ in,†
κ2 †d †d̂2 = − d̂2 + iG2 â − κ2 d̂2 ,dt2where γ is the mechanical damping rate and κ1,2 are theexternal damping rates of the two cavities. We will first neglect
the effect of internal damping rates κ1,2
, which are discussed
in
in detail in Sec. IV. In Eq. (3), d̂j are the input operators(d̂min = â in ), whose correlation functions are d̂jin (t)d̂k (t  ) =Nj δ(t − t  )δi,j . Here Nj is the average thermal population ofeach mode, and in the following discussion, we assume thatthe cavity thermal populations are always zero (due to the highfrequency) while keeping the mechanical thermal population
Nm finite.
The √frequency-resolved output modes d̂iout [ω] ≡ dωeiωtd̂iout (t)/ 2π are related to the input by d̂iout [ω] =3inj =1 Sij [ω]d̂j [ω] (j = 3 denotes the mechanical fields). The
scattering matrix S[ω] is obtained straightforwardly from the√
input-output relations d̂iout = κi d̂i + d̂iin [27] and the system
Langevin equations (with RWA) as follows:
in,†

 C2

√
C1 C2
χ1 χ2 χm4
 C1
−χ2 4 χ1 χm + 1
iC2
χ2 χm
2

χ χ −1
4 2 m

1⎜ √C1 C2
S[ω] = 1 +⎝ − 4 χ1 χ2 χm
C 1 χ1 χm − C 2 χ2 χm + 1
iC1
χ1 χm
2
with χi = 2κi /(κi − 2iω) (i ∈ {1,2,m}, κm ≡ γ ). Notice thatat ω = 0 the matrix only depends on the cooperativities.
Given the blue-detuned laser drive, the first question
involves the stability of our linearized system. The RouthHurwitz conditions [28] yield two requirements [14,26] toguarantee stability. The first is that γtot > 0, where γtot =γ + 4G21 /κ1 − 4G22 /κ2 is the effective damping rate of the mechanical resonator interacting with the two cavities. Focusing
on the interesting and relevant regime of strong cooperativities
Ci ≡ 4G2i /(γ κi )  1 and κi  γ , the two requirements canbe combined into
G21 /G22 > max(κ2 /κ1 ,κ1 /κ2 ).

(5)

III. BIPARTITE CAVITY OUTPUT ENTANGLEMENT

We start by discussing the entanglement oflight leaving the two cavities. The canonical mode

⎞

iC1
χ1 χm
2
⎟
iC2
− 2 χ2 χm ⎠ ,

(4)

−χm

operators
1D̂iout [ω,σ,τi ] = √
σ

ω+ σ2
ω− σ2



dω e−iω τi d̂iout [ω ](6)describe output temporal modes in a bandwidth σ centeredabout the frequency ω. τi sets the absolute time at which thewave packet of interest is emitted from cavity i; without lossof generality, we set τ2 = 0 . The covariance matrix of theoutput operators can be computed via the scattering matrix.
Then, the two-mode entanglement can be quantified using thelogarithmic negativity ENout [ω,σ,τ1 ] [29,30].
For clarity, we restrict ourselves in this section to thezero-bandwidth case σ → 0; the result is then independent ofτ1 (D̂iout [ω] ≡ D̂iout [ω,0,τi ]eiωτi , ENout [ω] ≡ ENout [ω,0,τ1 ]). We
first discuss a simple analytical solution of the bipartite cavity
output entanglement and then derive the optimal condition
for output entanglement. In the following section, we will
relax the assumption of zero bandwidth σ . We will also

013807-2

BIPARTITE AND TRIPARTITE OUTPUT ENTANGLEMENT . . .

PHYSICAL REVIEW A 91, 013807 (2015)

discuss other possible imperfections and their influence on
output entanglement: mechanical thermal noise, cavity internal
losses, and non-RWA corrections.

Entanglement EN

10.00
5.00

A. Mapping to a two-mode squeezed thermal state

Using the standard approach (i.e., plugging the system
covariance matrix into the definition of logarithmic negativity),
the entanglement of the three-mode system can be easily
computed numerically. However, we find that mapping theoutput state to a two-mode squeezed thermal state can yielda simple analytical expression which reveals a number ofinteresting properties of the output entanglement.
A two-mode squeezed thermal state can be written as

 †
(7)
ρ̂12 = Ŝ12 (R12 ) ρ̂1th (n̄1 ) ⊗ ρ̂2th (n̄2 ) Ŝ12 (R12 ).
Here



Ŝ12 (R12 ) = exp R12 D̂1out [0]D̂2out [0] − H.c.

(8)

is the two-mode squeeze operator, with R12 being the squeezing parameter, and ρ̂jth (n̄j ) describes a single-mode thermal
state with average population n̄j . The output state is thus
completely characterized by just three parameters: n̄1 , n̄2 , andR12 , and the entanglement of this state (if ENout  0) is simply
given by
ENout = − ln nR − n2R − (1 + 2n̄1 )(1 + 2n̄2 ) ,(9)with nR = (n̄1 + n̄2 + 1) cosh 2R12 . If any of the populationsn̄i is zero (which is the relevant case as discussed below),
we find that any finite squeezing R12 > 0 yields nonzero
entanglement.
From the Langevin equation (3), it is easy to verify
that D̂iout [ω] are all zero and only three correlators of
the output cavity modes are nonzero: (D̂1out [ω])† D̂1out [ω],(D̂2out [ω])† D̂2out [ω], and D̂1out [ω]D̂2out [ω]. Thus the covariance matrix of the system is the same as a two-mode squeezed
thermal state. Since two Gaussian states with the same
covariance matrix represent the same state, this output cavity
state can be mapped to a two-mode squeezed thermal state
whose nonzero correlation functions are simply out†D̂i [ω] D̂iout [ω] = n̄i + (n̄1 + n̄2 + 1) sinh2 |R12 [ω]|,eiθ[ω](n̄1 + n̄2 + 1)D̂1out [ω]D̂2out [ω] = −2× sinh 2|R12 [ω]|,Intracavity

1.00
0.50
0.10
0.05
0.0

where
D = (C1 + C2 )2 + 2 (C1 + C2 ) (1 + 2Nm ) + 1 + 4C1 C2 ,
E = 4 C1 C2 (C1 + C2 + 1 + 2Nm ) .
(13)

0.4

0.6

0.8

1.0

FIG. 2. (Color online) Entanglement of the two cavity output
fields on resonance (i.e., ω = 0; thick blue solid/dashed line) andcorresponding intracavity entanglement (thin red solid line) as
functions of C2 , with C1 = 4000 and ωm  κ (allowing use of theRWA). The orange dot denotes the value of ln 2C1 as in Eq. (18) with
κ1  κ2 . Except for the blue dashed line (mechanical bath occupancy
Nm = 100), results are for zero temperature. The output EN curves
only depend on Ci . For the intracavity curve (thin red solid line), we
have also assumed strong coupling, taking G1 : κ : γ = 100 : 10 : 1.

Notice that such a mapping is unique; that is, no other choiceof n̄1 , n̄2 , and R12 yields the same covariance matrix. From
Eqs. (9), (11), and (12), we see that the output entanglement
depends only on cavity cooperativities Ci = 4G2i /(γ κi ) andthe mechanical bath temperature (we always assume zerothermal occupancy for the two cavities).
Plugging Eqs. (11)–(13) into Eq. (9), one can obtain thefull expression of the output entanglement including thermal
effects. By focusing first on zero temperature we obtainn̄1 = 0,n̄2 = 4C2

γ2
4C2
=
,
(1 + C1 − C2 )2
γ 2tot(14)√
4 C1 C2 (C2 + C1 + 1)
tanh 2R12 = 2.
C1 + C22 + 6C1 C2 + 2(C1 + C2 ) + 1(15)

We stress that the effective thermal occupancies n̄j are notsimply equal to bath temperatures; in particular, n̄2 = 0 evenwhen all input noises are vacuum. The zero-temperature output
entanglement is


(12)

0.2

C2 /C1

ENout [0] = ln

(10)where n̄i ≡ n̄i [ω], R12 [ω] is complex, in general, and θ [ω] ≡arg(R12 [ω]). Using the correlators with σ = 0 and ω = 0 (i.e.,
output light near the cavity resonances), we find that R12 isreal and the output state can be characterized by√
E 2 − D 2 ∓ 4 (C2 + Nm ) 1±2Nm+n̄1,2 =
− ,C 1 − C2 + 1
2
2 (C1 − C2 + 1)2(11)
where the upper (lower) sign applies to n1 (n2 ), andtanh 2R12 = E/D,

Output field(1 + C1 − C2 )2√A + B + 2C2 (1 + 2C1 ) − 4 AB, (16)

with A = C2 (C1 + C2 ), B = (1 + C1 )2 + C1 C2 . Note that an
analytic expression for entanglement can also be obtained at
finite temperature, but it is too cumbersome to reveal muchinsight. In the following section, we will discuss in detail theeffect of finite temperature in two different regimes. The results
at both zero and finite temperature are plotted in Fig. 2.
Equations (9), (14), and (15) give us a simple way tounderstand entanglement generation. The ideal situation iswhen the two-mode squeeze parameter R12  1 and thethermal occupancy n̄i → 0. From Eq. (15), we see that alarge R12 generically requires large cooperativity. However,Eq. (14) indicates that this limit also yields a large effective
temperature for cavity 2 (heuristically, Ĥint turns vacuum noise
into thermal noise). This heating degrades the purity of the state013807-3

YING-DAN WANG, STEFANO CHESI, AND AASHISH A. CLERK[31]; however, as n̄1 remains zero, it only slightly degrades theentanglement [26].
B. Optimal output entanglement

Equation (16) shows that the zero-temperature entanglement is only a function of C1 , C2 ; thus weak coupling does
not prevent strong entanglement. For fixed C1 , the expression
is optimized when
C2 = C1 + 1 ≡ C2,opt .

(17)

Heuristically, this condition corresponds to having a totalmechanical damping γtot = 0. While ENout [0] appears to be only
a function of Ci , the ratio κ1 /κ2 also plays an independent role
via the stability condition of Eq. (5). If κ1  κ2 , C2,opt is also
the maximum value of C2 for which the system is stable. In
contrast, if κ2  κ1 , the system becomes unstable before C2
reaches C2,opt ; hence one cannot achieve the optimal amount ofentanglement. We thus see that, in addition to achieving large
Ci , it is advantageous to have κ1  κ2 . Figure 3 illustrates the
behavior of ENout [0] in these two cases (C2,opt is indicated asa gray vertical line). Note that a similar dependence on κ1 /κ2
was observed numerically in [14], but in a regime far from
optimal entanglement (i.e., for C2
C2,opt ).
In the large C1 limit, the maximum achievable entanglement
for the two cases reduces to  2C1
ln 1+2Nm ,κ1  κ2 ,
out
EN [0] max  κ2 −κ1 2
4κ22 κ1 Nm − ln κ2 +κ1 + C1 (κ1 +κ2 )3 , κ1 < κ2 ,
(18)

10
8 (a)
6
4
2
0
0

κ 1> κ 2 2.5 10
8 (b)
6
4
2
0
0
12000

C2

κ1< κ2 2.5

0 6748 12000

0Re(Ev) (MHz)

ENout[ω=0]with Nm = Nm (1 + κ2 /κ1 ) + 1. In both cases, the output
entanglement ENout is maximal at the boundary of systemstability, similar to the behavior of a nondegenerate parametricamplifier (NDPA) [32]. However, unlike a NDPA, ENout does
not diverge at this boundary. Further, while ENout is maximal

PHYSICAL REVIEW A 91, 013807 (2015)

at this boundary, the intracavity entanglement is extremely
suboptimal at this point. For κ1  κ2  γ and C2  C2,opt 
1, one has extremely large entanglement of the output fields
while simultaneously having almost zero entanglement of theintracavity fields (see Fig. 2).
IV. DISCUSSION OF IMPERFECTIONS

So far we have given the analytical solution and the optimal
condition for cavity output entanglement in the ideal case,where we have neglected the cavity internal loss, non-rotatingwave terms, and finite bandwidth. However, a real system
has various imperfections. We now focus on the case κ1 =
κ2 = κ, Cj  1 and carefully examine the influence of various
imperfections in two generic regimes.
The first is that of equal couplings G1 = G2 = G (i.e., C1 =C2 ), which, as discussed, is essentially optimal for maximizing
ENout [0] (it is also an ideal point to generate quantum-limited
amplification [33]). In this regime, the effective coupling G̃
in Eq. (2) vanishes, meaning that the three normal modes ofĤint are degenerate. Consequently, the σ = 0 entanglement
spectrum ENout [ω] has a single peak at ω = 0 [see thick curves
3/4
in Fig. 4(a)] with a width
in the weak-coupling case√of ∼γ 5C 1/12
(red top curve) and ∼ G(2κ γ )in the strong-couplingcase (green bottom curve). We stress that achieving large
optimal EN in this regime requires only strong cooperativities
and not the more stringent strong-coupling condition Gj > κj[see Eq. (18)].
Keeping κ1 = κ2 and Cj  1, another generic regime iswhere G2 /G1 is sufficiently smaller than 1 such that theeffective coupling G̃ is larger than κ,γ ; this necessarily
requires G1 > κ. In this regime, the three normal modes ofĤint are spectrally resolved, and ENout [ω] correspondingly has
three peaks [14] [see the thin blue curve of Fig. 4(a)]. The
entanglement at ω = 0 is necessarily much smaller than the
optimal value in Eq. (18). One finds that C1  1 is not by itself
enough to ensure large EN in this regime; one also needs tobe deep in the strong-coupling regime, G1  κ. As discussed
in [14], this “resolved-modes” regime does, however, offer
enhanced protection against mechanical thermal noise, as thecentral peak is due to the mechanically dark normal mode βB .
A. Thermal resilience of entanglement

C2

FIG. 3. (Color online) Output entanglement of the two cavity
fields at resonance (i.e., ω = 0) as a function of C2 for C1 = 12 000
in two different cases: (a) κ1 > κ2 (κ1 = 1.5κ2 = 2π × 0.3 MHz)and (b) κ1 < κ2 (κ1 = 0.75κ2 = 2π × 0.3 MHz). If the instability
condition were not considered, the entanglement would reach amaximum finite value at C2 = C1 + 1 (gray vertical lines). The
green lines show the real part of the three eigenvalues (Ev ) of
the susceptibility matrix; the system becomes unstable when oneof the real parts becomes positive, which is indicated by the blue
dashed lines. The red dots indicates the maximum entanglementin each case by also considering the stability constraints. We usedγ = 2π × 10 KHz. The dashed portion of the entanglement curve(after it intersects with the blue vertical line) indicates that the system
has entered the unstable regime.

Thermal fluctuation severely limits the coherence of thelow-frequency mechanical resonator, and it is important
to improve the thermal resilience for protocols where themechanics acts as a coherent mediator, such as for statetransfer [34,35]. The bipartite entanglement mediated by
the mechanics also suffers from mechanical thermal noise.
While the analytical solution including thermal noise has been
derived in the previous section, we discuss here the thermal
resilience of the output entanglement in the two representative
regimes mentioned above in order to elucidate the function of
the mechanically dark mode [14,34–36].
For the first regime G1 = G2 = G and κ1 = κ2 = κ (i.e.,
C1 = C2 ), the entanglement value at ω = 0 is
2C
+ O(C −1 ),ENout [0] ≈ ln
(19)
1 + 2Nm

013807-4

BIPARTITE AND TRIPARTITE OUTPUT ENTANGLEMENT . . .

10

15

0

ω/κ

normal modes requires G1 > κ cosh r  κ, i.e., deep in thestrong-coupling regime for a large r (r is defined for the
Bogoliubov modes after Eq. (2)). For typical parameters, themaximum entanglement in the resolved-peak regime is far
less than the optimal value that would be achieved if C2 wereincreased to C1 [see Fig. 4(b)]. On the other hand, the output
entanglement at the central peak can be written as (assuming
Nm e2r
C̃)

(b)

10

0.1

5

-15 -10 -5 0

Entanglement ENout[0]

8
6
4
2
0
-0.1 0(a)outEntanglement EN [ω]

15

5
0
0

5 10 15

PHYSICAL REVIEW A 91, 013807 (2015)

ENout [0] ≈ 4r − 2e2r (2Nm + 1)/C̃,

200 400 600 800 1000

Nm

(20)

with the effective cooperativity C̃ = 4G̃ /γ κ. This showsthat, as the central peak corresponds to the dark mode, themechanical noise is suppressed by the effective cooperativity
and the entanglement is robust to thermal fluctuations [14], asshown in Fig. 4(b). The entanglement at the side peaks (brightmodes involving the mechanics) is still sensitive to thermal
noise and in the large r limit:



4γ
1
out
EN [ω = ±G̃] ≈ 2r − ln
Nm +
,
(21)
κ
2

(c)

12 7
10 6
0
8
6

40

ωm /κ

4
2
0 -2
10 0.1

1

80

11 0 102 103

ωm /κ

(d)

12

out

Entanglement ENout[0]14 8Entanglement EN [0,σ,τ1]

2

8
4
0
10-7 10-5

10-3

σ/ κ

0.1

1

which is valid if γ Nm /κ

FIG. 4. (Color online) Different aspects of the output entanglement in 3 different cases: strong equal-coupling (red thick upper
line, G1 = G2 = 13.3κ), weak equal-coupling (green thick lower
line, G1 = G2 = 0.1κ) and resolved normal-modes (blue thin line,
G1 /κ = 13.3 and G2 /κ = 6.7. For the blue and red curves, γ /κ =1.67 × 10−3 while for the green curve γ /κ = 3.3 × 10−5 . (a) Output
spectrum in the limit of small mode bandwidth (σ → 0) as a functionof mode center frequency ω. The inset shows a zoom-in of thegreen curve (weak coupling case). (b) Thermal resilience. The output
entanglement at ω = 0 is plotted vs Nm for the three cases. (c)The influence of counterrotating terms. For the three cases, thecurves show how the output entanglement ĒNout [ω = 0] including thecounterrotating terms approaches the RWA results neglecting thoseterms (marked by the dots) as ωm /κ increases. The inset shows a good
fit of the analytical expression (29) (dashed lines) with the numericalresult (solid line) in the weak-equal-coupling case, and the gray line
indicates the RWA result. (d) Effects of nonzero mode bandwidth σon the output entanglement of modes with center frequency ω = 0.
The solid lines are for zero time delay between the two cavity output
modes [τ1 = 0; see Eq. (6)], while the dashed lines are the result
including an optimal time delay τ1 = κ/(4G2 ).

We stress that a robust peak in ENout [ω] remains in the weakcoupling case G < κ as long as C  1.
The entanglement at ω = 0 in the equal-coupling regime
shows a prominent decrease with temperature due to themechanical thermal noise. At low temperature, the equalcoupling regime (G1 = G2 ) yields a decent entanglement evenfor weak coupling (κ > G), as shown by the green curve in
Fig. 4(b).
√We turn now to another regime: G̃ = G21 − G22 > κ,
where the three normal modes of the interaction Hamiltonianin Eq. (1) are spectrally resolved. In this regime (as discussed
in [14]), ENout [ω] has three resolved peaks, with each peakcorresponding to a normal mode. The “dark mode” βB
corresponds to the peak at ω = 0, while the peaks at ±G̃ arethe coupled “hybrid” modes [see thin blue curve in Fig. 4(a)];all have width ∼κ. In the simple case κ1 = κ2 , resolving the1.

B. Influence of cavity internal losses

Except the loss through output channels, there is also
extra zero-point noise incident on the cavity which can be
described by the internal damping rate κi . This can be treatedas an additional independent input channel in the Langevinequations for the cavities. In contrast to the generation of
intracavity entanglement [26], the output entanglement issensitive to the internal losses of the cavities [see Eq. (3)].
a. Equal-coupling case. With finite internal loss κ1 = κ2 =κ  and also assuming G  γ , κ, κ  , the entanglement reads


1
κ
out
+
,
(22)
EN [0] ≈ − lnκ + κ
2C
which recovers the two limits
ENout [0] ≈ ln 2C(23)

at small internal loss (κ/κ   C  1) andENout [0] ≈ ln(1 + κ/κ  )(24)

at large internal loss (κ/κ 
2C). In the case of a cavity
with a tunable
external
dampingrate, the optimal κ satisfies√κ = κ  ( 2C  − 1), with C  = 4G2 /κ  γ . The corresponding
entanglement is ENout [0] ≈ 12 ln(C  /2). This result is the sameas that for the optimal intracavity entanglement, as shownin Eq. (11) of [26] (notice that the definition of C in thispaper differs by a factor of 4 compared to that in [26]).
Thus, the maximum entanglement is the same for both spectral
entanglement and intracavity entanglement.
b. Resolved-peak case. Including internal loss and assuming κi = κ  , G̃  κ, γ , κ  , we have
κ + e4r κ 
er sinh r.
(25)
ENout [0] ≈ 4r − ln+
κ + κ
C̃
When κ/κ 
e4r C̃/(e2r + C̃), this reduces to ENout [0] ≈ln(1 + κ/κ ). The maximum of entanglement in Eq. (25) is013807-5

YING-DAN WANG, STEFANO CHESI, AND AASHISH A. CLERK

out[0] = 4r − ln
EN,opt1+



2 e4r − 1
1 − e2r+√C̃ e2r + 1
C̃

5


,
(26)

with

C. Influence of counterrotating wave terms

So far, we have only discussed the dynamics with the RWA,
considering the good cavity limit κωm . Here we will givethe full result including the non-RWA terms and show that theRWA condition can be more stringent in certain regimes.
First, we notice that the counterrotating terms are time
independent in the rotating frame with respect to the cavity
drives. In this frame, the full Hamiltonian is written as††
Ĥ = ωm (b̂† b̂ + d̂1 d̂1 − d̂2 d̂2 ) + Ĥint + ĤCR , with
Ĥint = G1 b̂† d̂1 + G2 b̂d̂2 + H.c.,(27)††
ĤCR = (G1 b̂† d̂1 + G2 b̂d̂2 ) + H.c.
Thus, a closed set of equations in the frequency domain canbe obtained:γ†− i(ω − ωm )b̂[ω] = − b̂[ω] − iG1 (d̂1 [ω] + d̂1 [ω])2√†
− iG2 (d̂2 [ω] + d̂2 [ω]) − γ b̂in [ω],κ1−i(ω − ωm )d̂1 [ω] = − d̂1 [ω] − iG1 (b̂[ω] + b̂† [ω])2√
− κ1 d̂in,1 [ω],
(28)
κ2 †
†
−i(ω − ωm )d̂2 [ω] = − d̂2 [ω] + iG2 (b̂[ω] + b̂† [ω])2√ †− κ2 d̂in,2 [ω].
These equations can be solved analytically. Figure 4(c) shows
a comparison of the results for ENout [ω = 0] with and without
RWA. The lower curve corresponds to the case where G̃ > κand the normal modes are resolved: counterrotating terms
suppress the entanglement but become insignificant onceωm > κ. The two upper curves correspond to the case withequal coupling. In the good cavity limit, the maximum entanglement in the equal-coupling regime is much larger than the
resolved-peak regime. Consequently, non-RWA corrections
play a larger role, and one can only achieve the RWA result
deep into the good cavity limit ωm  κ.
For both equal-coupling and resolved-peak cases, theleading non-RWA correction to ENout [0] is
κ2
.
2
16ωm

3
2
0

However, notice that in order to have resolved peaks, G̃  κ +κ  . This means the optimal condition is normally not satisfied
unless κ 
γ.

out

(b)
0.5

1


κopt = κ ( C̃(1 + e2r )/2 − 1).
ĒNout [0] − ENout [0] ≈ −eEN [0]

(a)

4

θ/π



|R12|

obtained as follows:

PHYSICAL REVIEW A 91, 013807 (2015)(29)

Here ĒNout [0] is the entanglement including the counterrotating
terms. This approximate expression is shown by the orange
dashed line in the inset of Fig. 4(c). Thus the condition to

0

-0.5
-0.02

0 0.02
ω/κ

-0.02

0
ω/κ

0.02

FIG. 5. (Color online) The dependence of the squeezing parameter R12 [ω] = |R12 |eiθ on frequency. The left panel shows themagnitude change of R12 , while the right panel shows the phasechange of R12 . The dashed line in the right panel is the linear
dependence given by Eq. (31). The parameters used are G/κ = 0.1,
γ /κ = 3.3 × 10−5 .

justify the use of RWA is

 
out
 ωm 
1
e EN
 
,
 κ 
4 ENout(30)

which is looser in the case of resolved peaks (due to the smaller
ENout ).
D. Time delay to improve the entanglement bandwidth

So far, we have considered the output entanglement withzero bandwidth, i.e., ENout [ω,σ = 0,τ1 ]. However, in practice,
the filter function has a finite bandwidth. Assuming the center
frequency is set at ω = 0 (the cavity resonance frequency
in the laboratory frame), the entanglement has a nontrivial
dependence on the bandwidth σ and the relative time delay τ1(taking τ2 = 0).
The solid lines of Fig. 4(d) show the entanglement of the two
output cavity modes for zero time delay ENout [ω = 0,σ,τ1 = 0]
as a function of bandwidth σ . While the equal-coupling caseyields large entanglement at σ = 0, it is much more sensitive
to the increase in the bandwidth. ENout [0,σ,0] decays on a scaleσ ∼ C −1/4 γ . In contrast, in the resolved normal-mode case,
ENout [0,σ,0] is less sensitive to increasing the mode bandwidthand is only suppressed significantly when σ ∼ κ.
The strong sensitivity to nonzero σ in the case of equal
coupling is related to the change in the squeezing phaseat different ω. As discussed before, the cavity output state
can be characterized by a two-mode squeezed thermal state[see Eq. (7)] with a complex squeezing parameter R12 [ω]whose phase is θ [ω] = arg−D̂1out [ω]D̂2out [ω] [see Eq. (10)].
At ω = 0 and σ = 0, R12 [0] is real [see Eq. (15)], i.e.,θ = 0. For ω = 0, a frequency-dependent phase arises [seeFig. 5(b)]. In the large-squeezing limit, i.e., |R12 |  1, sucha phase variation leads to a rapid decrease in entanglement
as σ increases. Assuming κ1 = κ2 and G1 = G2 , where themagnitude of squeezing is maximized, the peak width of|R12 |[ω] is comparable to that of ENout [ω]. The corresponding
phase variation is approximately linear in the vicinity of ω = 0[as shown by the dashed line in Fig. 5(b)],δθ ≈ (κ/4G2 )δω.

(31)

This suggests [as per Eq. (6)] that the optimal entanglement
is between a cavity-2 mode emitted at τ2 = 0 and a cavity-1

013807-6

BIPARTITE AND TRIPARTITE OUTPUT ENTANGLEMENT . . .

mode emitted at τ1 = κ/4G2 . The entanglement with timedelay is shown by the red dashed curve in Fig. 4(d). On aheuristic level, one could think that our system first generates
entangled phonon-cavity-2 photon pairs via the G2 interaction
in Eq. (2); next, the G1 interaction swaps the phonon state
into the photon state in cavity 1. This physical picture is inagreement with finding that the optimal entanglement involves
a positive delay for the cavity-1 output mode, as in Fig. 4(d).
Finally, we note that for ω away from zero, the squeezing
phase θ has a nonlinear frequency dependence, and hence theoptimal filter function for larger σ will not correspond to asimple delay as in Eq. (31).
V. GENUINE TRIPARTITE ENTANGLEMENT OF THE
ITINERANT MODES

As mentioned in the Introduction, alternative methods exist
to access the mechanical output field. We discuss them in detailbelow, and motivated by these possibilities, we analyze themultipartite entanglement properties of all three output fields inour system. Note that previous work studied the nonstationary
tripartite entanglement of intracavity fields generated by theclosed-system Hamiltonian Ĥint in Eq. (2) [37]; in contrast, ourfocus here is on the steady-state output entanglement in our
dissipative three-mode system.
A. Phonon output through an auxiliary cavity

Experimental setups using optomechanical crystals [20,38]
can access the mechanical output field via engineered phonon
waveguides. For a general optomechanical system where aphonon waveguide is absent, we find it is still possible toaccess the phonon output by having its damping dominated by
the optical damping of an auxiliary cooling cavity (see Fig. 1).
This auxiliary cavity will have a large damping rate and will
be coupled to the mechanics via a linearized optomechanicalcoupling (the cavity is strongly driven by a red-detuned laser):Ĥa = Ga (b̂d̂a† + b̂† d̂a ),(32)

where d̂a is the annihilation operator of the auxiliary cavity and
Ga is the corresponding coupling. The auxiliary cavity, due tothe large damping rate κa , can be described by the following
steady-state relation:d̂a = −i22Gab̂ − √ d̂ain .
κa
κa

(33)

Plugging this into the equation of motion for the mechanicalresonator, one obtains
2G2
d
2iGa
†
b̂ = − a b̂ − i(G1 d̂1 + G2 d̂2 ) + √ d̂ain ,
dt
κa
κa

(34)

where we have neglected other damping channels of theresonator, supposing that G2a /κa is sufficiently large. Comparing Eq. (34) with the Langevin equation without an√†
auxiliary cavity, dtd b̂ = − γ2 b̂ − i(G1 d̂1 + G2 d̂2 ) − γ d̂ain , we
can identifyb̂in = −i d̂ain , γ = 4G2a /κa .

PHYSICAL REVIEW A 91, 013807 (2015)√d̂aout = −i γ b̂ − i b̂in ; thusd̂aout = −i b̂out .

(36)This shows that a strongly damped auxiliary cavity can serveas the output of the mechanical mode.
B. Three-mode output state

We focus on zero frequency and zero bandwidth andconsider the ideal case where all dissipative baths are at zerotemperature. In this case, the three-mode output state is apure state | s . Since n̄1 = 0 [see Eq. (14)], cavity 1 must
be in the ground state if we consider the “unsqueezed” state| s  ≡ S̃12 (−R12 )| . However, | s  has residual entanglement between cavity 2 and the mechanical mode. This is clear
from the finite value of n̄2 , which diverges at the instability
condition γtot = 0 (the finite population n̄2 is a consequence of
tracing out the mechanical mode). One can compute explicitly
the correlations of the state | s  and identify it as a squeezed
vacuum of cavity 2 and the mechanical mode, with squeezing
parameter R2m = sinh−1 n̄2 . Hence the output state can be
written as a twice-squeezed vacuum|  = Ŝ12 (R12 )Ŝ2m (R2m )|000,

(37)

out
0]D̂2out [0] − H.c.) is the twowhere Ŝ2m (R2m ) ≡ exp(iR2m D̂mmode squeeze operator entangling the output of the mechanicsoutat ω = 0 (denoted by D̂m[0], which is defined in a way
similar to that of the cavity output) with that of cavity 2.
Since the effective temperature represented by n̄2 diverges at
the instability point [see Eq. (14)], the squeezing parameter
R2m = sinh−1 n̄2 diverges as well. Equation (37) demonstrates
that the effective temperature n̄2 which degraded the cavitycavity entanglement in Eq. (14) is a direct consequence ofentanglement between cavity 2 and the mechanics. It also
demonstrates the asymmetry between the three modes (i.e.,there is no direct squeezing between the mechanical andcavity-1 outputs).
Furthermore, Eq. (37) shows that we have genuine tripartite
entanglement (GTE): none of the parties can be separated fromany other party in a mixture of product states, implying a fully
inseparable state [39,40]. To see this, note that the total stateis pure, and if one traces over one subsystem, the remaining
two are in a mixed state. GTE is also evident from writing thestate in the Fock-state basis |n1 ,n2 ,nm ,qpp Cp+q Nm  2 N1  2| =|p,p + q,q,(38)(1 + N2 )(p+q+1)/2
pq
pwhere Cp+q are binomial coefficients and Ni  =(D̂iout [0])† D̂iout [0] (i = 1,2,m) is the average photon orphonon number of each mode

(35)

The input-output relation of the auxiliary cavity is d̂aout =√ind̂a + κa d̂a , which, together with Eqs. (33) and (35), gives
013807-7

4C1 C2,(1 + C1 − C2 )24C2 (C1 + 1)
N̂2  =,(1 + C1 − C2 )2
4C2
N̂m  =.
(1 + C1 − C2 )2
N̂1  =(39)

YING-DAN WANG, STEFANO CHESI, AND AASHISH A. CLERKPHYSICAL REVIEW A 91, 013807 (2015)

nonzero and diverges as one increases C2 to the instability
point C1 + 1 (due to the divergence of n̄2 and R2m ). These
calculations also reveal the absence of any direct entanglement
between the mechanical and cavity-1 outputs.
The Rényi-2 entropy is given by S2 (ρ) = − ln Tr(ρ 2 ) andallows us to define an entanglement measure ε2 (ρA:B ) for
bipartite states ρAB . For pure states one simply has ε2 (ρA:B ) =S2 (ρA ), with ρA being the reduced density matrix of subsystem
A. The tripartite entanglement is then characterized throughthe residual entanglement ε2 (ρi:j :k ), given byε2 (ρi:j :k ) = ε2 (ρi:j k ) − ε2 (ρi:k ) − ε2 (ρi:j )  0,
FIG. 6. (Color online) Three-mode entanglement, evaluated by
the residual GR2 entanglement with cavity 1 as the focus mode
2 (ρ1:2:3 ), vs the cooperativity of cavity 2. Here C1 is fixed to10, 50, 100, 200 (from top to bottom). The gray dashed lines indicate
the onset of system instability. The three-mode entanglement diverges
at γtot = 0 or C2 = C1 + 1. The C1 = 100 (red) curve is replotted in
the inset, together with the two residual GR2 entanglements where
the focus mode is cavity 2 [2 (ρ2:1:3 )] or mechanics [2 (ρ3:2:1 )].

|  only involves Fock states |n1 ,n2 ,nm , where n2 = n1 +
nm ; there is thus a perfect correlation between the three systems
that is only evident by looking at all three modes. If any
two modes are traced out, the remaining third mode is in an
impure thermal state. Another interesting observation is that,
although the tripartite steady-state equation (37) is generated
by interaction between intracavity modes and the input-output
relations, the output state described by Eq. (37) could also be
achieved by a unitary evolution generated by the followinginteraction between the output modes:
 out
†
out
H̃ = R12 sin R2m D̂m[0]D̂2out [0] − R2m D̂m[0] D̂1out [0]+ iR12 cos R2m D̂1out [0]D̂2out [0] + H.c.

(40)

This result can be obtained by Wei-Norman decompoout
sition in terms of a closed algebra, {D̂m[0]† D̂1out [0] +out[0]D̂2out [0] + H.c.,D̂1out [0]D̂2out [0] − H.c.}.
H.c.,D̂mwhere i = j = k = 1,2,3 represent the three different modes(the third mode denotes the mechanics in our case). In Eq. (41),ε2 (ρi:j k ) is the bipartite entanglement partitioning the global
system into A = i and B = j k, while ε2 (ρi:j ) and ε2 (ρi:k )consider the reduced density matrices of subsystems ij andik, respectively. In general, three different values of ε2 (ρi:j :k )are obtained, depending on the choice of the “focus mode” i.
There are special cases when ε2 (ρi:j :k ) is invariant under mode
permutation [21]; as we explain below, that is not the case forour system.
For pure tripartite Gaussian states, ε2 (ρi:j :k ) can be evaluated analytically [21]. The first term in Eq. (41) isε2 (ρi:j k ) = S2 (ρi ) = ln ai ,

ai =

β
,
8ak2
⎪
⎪
2
2

⎪
⎩ ai 2−aj 2 ,
ak −1



det Vi .

(43)

Using the covariance matrix for our system, calculated
from the Langevin equation, it is possible to obtain explicit
formulas for ε2 (ρ1:23 ), ε2 (ρ2:13 ), and ε2 (ρ3:12 ) in terms of thecooperativities C1,2 . Although we omit them here, it is worth
mentioning that these three quantities are all nonzero, whichshows that none of the three systems are separable.
The last two terms in Eq. (41) are given by [21]ε2 (ρj :k ) = 12 ln gi (i = j = k),

The tripartite entanglement can be measured by the residue
Gaussian Rényi-2 entanglement [21]. We find that it is indeed

gi =

(42)

where ai is related to the covariance matrix Vi of subsystem i:

C. Three-mode entanglement measure basedon Rényi-2 entropy

⎧
⎪
1,
⎪
⎪
⎨

(41)

where


ai2 + aj2 − 1,αk < ak < ai2 + aj2 − 1,ak 

ak  αk ,

with
"
# 

 


#
2
2
# 2 ai2 + aj2 + ai2 − aj2 + ai2 − aj2  ai2 − aj2 + 8 ai2 + aj2

αk = $,
2 ai2 + aj2√
β = 2a12 + 2a22 + 2a32 + 2a12 a32 + 2a22 a32 + 2a12 a22 − a14 − a24 − a34 − δ − 1,
δ = [(a1 − a2 − a3 )2 − 1][(a1 + a2 − a3 )2 − 1][(a1 − a2 + a3 )2 − 1][(a1 + a2 + a3 )2 − 1].
013807-8(44)

BIPARTITE AND TRIPARTITE OUTPUT ENTANGLEMENT . . .

PHYSICAL REVIEW A 91, 013807 (2015)

Since ε2 (ρj :k ) is the GR2 two-mode entanglement measure in the subspace where mode i is eliminated, it serves as analternative way to quantify the two-mode entanglement in our system, besides the logarithmic negativity used in Sec. III. We
obtain all the two-mode entanglements as follows: (1 + C2 )2 + C12 + 2C1 + 6C1 C2
C22 + 6C1 C2 + C121,
ε2 (ρ1:2 ) = ln
≈ ln
+O
2
2
2
C
(1 + C2 ) + C1 + 2C1 − 2C1 C2(C1 − C2 )2 
(1 + C1 )2 + C22 + 6C2 + 2C1 C2
1,(45)
ε2 (ρ3:2 ) = ln
≈OC2
(1 − C2 )2 + C12 + 2C1 + 2C1 C2ε2 (ρ1:3 ) = 0,
where the last approximation is in the large-C limit. These results show that the entanglement between cavity 1 and the mechanics
is always zero, although the mechanics is entangled with the composite system of the two cavities [ε2 (ρ3:12 ) = 0]. In particular,
there is entanglement between the mechanics and cavity 2, although it is much smaller than the entanglement between cavities 1
and 2. These results are all in agreement with those based on the logarithmic negativity.
We now turn to the evaluation of the tripartite entanglement ε2 (ρi:j :k ), which is permutationally invariant only when all ε2 (ρi:j k )and ε2 (ρi:j ) are nonzero [21]. The fact that ε2 (ρ1:3 ) = 0 implies that ε2 (ρi:j :k ) depends on the focus mode; thus we consider all
three choices, i = 1,2,3. Using the results for ε2 (ρi:j k ), ε2 (ρi:j ) discussed above, Eq. (41) gives



 2
2C2
C1 + (1 + C2 )2 + 2C1 − 2C1 C2
−1
−
2
tanh
,
ε2 (ρ1:2:3 ) = ln
(1 + C1 − C2 )2(1 + C1 )2 + 6C1 C2 + C22
 2
C2 + (1 + C1 )2 + 6C2 (1 + C1 )ε2 (ρ2:1:3 ) = ln
C22 + (1 + C1 )2 + 2C2 (3 + C1 )
(46)
 2

[C1 + (1 + C2 )2 − 2C1 (C2 − 1)][C12 + (1 − C2 )2 + 2C1 (1 + C2 )]+ ln,(1 + C1 − C2 )2 [C12 + (1 + C2 )2 + 2C1 (1 + 3C2 )]



 2
4C2
C2 + (1 + C1 )2 + 6C2 − 2C1 C2
−1
−
2
tanh
.
ε2 (ρ3:1:2 ) = ln
(1 + C1 − C2 )2(1 + C1 − C2 )2As announced, the three residual GR2 entanglements areunequal. Nevertheless, they are all larger than zero, which
confirms the presence of genuine tripartite entanglement.
Furthermore, as shown in Fig. 6, all ε2 (ρi:j :k ) diverge at the
instability point γtot = 0. This divergence is similar to the caseof a parametric amplifier and is related to the divergences ofn̄2 and R2m discussed in subsection V B.
VI. CONCLUSION

We have studied in detail the bipartite and tripartite entanglement of the output fields in a three-mode optomechanicalsystem. For bipartite photonic entanglement, we have found an
explicit analytical expression based on identifying the output
with a squeezed thermal state. This has allowed us to findan optimal “impedance matching” condition for entanglement
generation: tuning the system parameters to achieve thiscondition yields a dramatic increase of the maximum output
entanglement. Strikingly, our optimal condition corresponds toa regime where entanglement between the two optical cavities
is nearly zero [26] but at the same time strongly entangled
photons are emitted into the output channels.
We also show that large values of the output entanglement
can be achieved even in the standard “weak-coupling” regime
with large cooperativity (the same regime used for cavity[1] T. Aoki, N. Takei, H. Yonezawa, K. Wakui, T. Hiraoka, A.
Furusawa, and P. van Loock, Phys. Rev. Lett. 91, 080404 (2003).

cooling in numerous experiments): one simply needs totune parameters to achieve the aforementioned impedance
matching condition. Furthermore, we have considered realistic
conditions including different types of imperfections: finite
temperature, internal losses, sizable counterrotating terms, andfinite mode bandwidth. By restricting ourselves to the optimal
equal-coupling regime, we find that it is crucial to consider
time delays of emitted wave packets; we find the optimalvalue of this delay and show that it substantially improves
the attainable entanglement. Finally, we show how to generate
genuine tripartite entanglement of spatially separated, itinerant
quanta: we provide a simple analytic characterization of theoutput state and use a rigorous entanglement monotone toquantify the three-mode entanglement. Our results can also be
applied to other parametrically coupled three-mode bosonicsystems, in addition to optomechanical systems.
ACKNOWLEDGMENTS

We thank M. Woolley for useful discussions. This work
was supported by the DARPA ORCHID program under agrant from the AFOSR. S.C. and Y.D.W. acknowledge support
from the Chinese Youth 1000 Talents Program. Y.D.W. also
acknowledges the grant from NSFC (No. 11434011).

[2] J.-W. Pan, Z.-B. Chen, C.-Y. Lu, H. Weinfurter, A. Zeilinger,
and M. Żukowski, Rev. Mod. Phys. 84, 777 (2012).

013807-9

YING-DAN WANG, STEFANO CHESI, AND AASHISH A. CLERK[3] J. Raimond, M. Brune, and S. Haroche, Rev. Mod. Phys. 73, 565(2001).
[4] M. Neeley, R. C. Bialczak, M. Lenander, E. Lucero,
M. Mariantoni, D. Sank, H. Wang, M. Weides, J. Wenner, Y. Yin,
T. Yamamoto, A. N. Cleland, and J. M. Martinis, Nature(London) 467, 570 (2010).
[5] L. DiCarlo, M. Reed, L. Sun, B. L. Johnson, J. M. Chow, J. M.
Gambetta, L. Frunzio, S. M. Girvin, M. H. Devoret, and R. J.
Schoelkopf, Nature (London) 467, 574 (2010).
[6] E. Flurin, N. Roch, F. Mallet, M. H. Devoret, and B. Huard,
Phys. Rev. Lett. 109, 183901 (2012).
[7] H. Bernien, B. Hensen, W. Pfaff, G. Koolstra, M. S. Blok,
L. Robledo, T. H. Taminiau, M. Markham, D. J. Twitchen,
L. Childress, and R. Hanson, Nature (London) 497, 86(2013).
[8] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt,
arXiv:1303.0733.
[9] M. Paternostro, D. Vitali, S. Gigan, M. S. Kim, C. Brukner,
J. Eisert, and M. Aspelmeyer, Phys. Rev. Lett. 99, 250401(2007).
[10] C. Genes, A. Mari, P. Tombesi, and D. Vitali, Phys. Rev. A 78,
032316 (2008).
[11] C. Wipf, T. Corbitt, Y. Chen, and N. Mavalvala, New J. Phys.
10, 095017 (2008).
[12] S. G. Hofer, W. Wieczorek, M. Aspelmeyer, and K. Hammerer,
Phys. Rev. A 84, 052327 (2011).
[13] Sh. Barzanjeh, M. Abdi, G. J. Milburn, P. Tombesi, and D. Vitali,
Phys. Rev. Lett. 109, 130503 (2012).
[14] L. Tian, Phys. Rev. Lett. 110, 233602 (2013).
[15] M. C. Kuzyk, S. J. van Enk, and H. Wang, Phys. Rev. A 88,
062341 (2013).
[16] T. A. Palomaki, J. D. Teufel, R. W. Simmonds, and K. W.
Lehnert, Science 342, 710 (2013).
[17] C. Dong, V. Fiore, M. C. Kuzyk, and H. Wang, Science 338,
1609 (2012).
[18] J. T. Hill, A. H. Safavi-Naeini, J. Chan, and O. Painter, Nat.
Commun. 3, 1196 (2012).

PHYSICAL REVIEW A 91, 013807 (2015)[19] R. Andrews, R. W. Peterson, T. P. Purdy, K. Cicak, R. W.
Simmonds, C. A. Regal, and K. W. Lehnert, Nat. Phys. 10,
321 (2014).
[20] A. H. Safavi-Naeini and O. Painter, New J. Phys. 13, 013017(2011).
[21] G. Adesso, D. Girolami, and A. Serafini, Phys. Rev. Lett. 109,
190502 (2012).
[22] S. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513 (2005).
[23] S. Koike, H. Takahashi, H. Yonezawa, N. Takei, S. L. Braunstein,
T. Aoki, and A. Furusawa, Phys. Rev. Lett. 96, 060504 (2006).
[24] N. Bergeal, R. Vijay, V. E. Manucharyan, I. Siddiqi, R. J.
Schoelkopf, and M. Devoret, Nat. Phys. 6, 296 (2010).
[25] A. Baust et al., arXiv:1405.1969.
[26] Y.-D. Wang and A. A. Clerk, Phys. Rev. Lett. 110, 253601(2013).
[27] C. Gardiner and P. Zoller, Quantum Noise, 3rd ed. (Springer,
New York, 2004).
[28] E. X. DeJesus and C. Kaufman, Phys. Rev. A 35, 5288 (1987).
[29] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
[30] M. B. Plenio, Phys. Rev. Lett. 95, 090503 (2005).
[31] M. J. Woolley and A. A. Clerk, Phys. Rev. A 89, 063805 (2014).
[32] D. F. Walls and G. J. Milburn, Quantum Optics (Springer, Berlin,1994).
[33] A. Metelmann and A. A. Clerk, Phys. Rev. Lett. 112, 133904(2014).
[34] Y.-D. Wang and A. A. Clerk, Phys. Rev. Lett. 108, 153603(2012).
[35] Y.-D. Wang and A. A. Clerk, New J. Phys. 14, 105010 (2012).
[36] L. Tian, Phys. Rev. Lett. 108, 153604 (2012).
[37] A. Ferraro, M. G. A. Paris, M. Bondani, A. Allevi, E. Puddu,
and A. Andreoni, J. Opt. Soc. Am. B 21, 1241 (2004).
[38] M. Eichenfield, J. Chan, R. M. Camacho, K. J. Vahala, and O.
Painter, Nature (London) 462, 78 (2009).
[39] G. Giedke, B. Kraus, M. Lewenstein, and J. I. Cirac, Phys. Rev.
A 64, 052303 (2001).
[40] P. van Loock and S. L. Braunstein, Phys. Rev. Lett. 84, 3482
(2000).

013807-10