PHYSICAL REVIEW A 87, 063846 (2013)

Two-mode back-action-evading measurements in cavity optomechanics
M. J. Woolley1,2 and A. A. Clerk1
1

2

Department of Physics, McGill University, Montréal, Québec, Canada H3A 2T8
School of Engineering and Information Technology, UNSW Canberra, Canberra, Australian Capital Territory 2600, Australia(Received 22 April 2013; published 28 June 2013)
We study theoretically a three-mode optomechanical system where two mechanical oscillators are coupled
to a single cavity mode. By using two-tone (i.e., amplitude-modulated) driving of the cavity, it is possible
to couple the cavity to a single collective quadrature of the mechanical oscillators. In such a way, a backaction-evading measurement of the collective mechanical quadrature is possible. We discuss how this can allow
one to measure both quadratures of a mechanical force beyond the full quantum limit, paying close attention
to the role of dissipation and experimental imperfections. We also describe how this scheme allows one togenerate steady-state mechanical entanglement; namely, one can conditionally prepare an entangled, two-mode
squeezed mechanical state. This entanglement can be verified directly from the measurement record by applying
a generalized version of Duan’s inequality; we also discuss how feedback can be used to produce unconditional
entanglement.
DOI: 10.1103/PhysRevA.87.063846

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

I. INTRODUCTION

The idea of a back-action-evading quantum measurement
is by now a well-known concept [1–4]. The simplest and moststudied realization is the continuous position measurement ofa simple harmonic oscillator. If one tries to simultaneously
follow both quadratures of the oscillator’s motion (i.e., boththe amplitudes of the sine and cosine components of the
motion), then the unavoidable effects of quantum back-action
imply that one cannot improve the measurement indefinitely
by increasing the measurement strength. This leads to thequantum limit on continuous position detection [3,5]. If instead
one measures just a single quadrature of the motion (say theX̂ quadrature), there is no such limit as one has “evaded”
the back-action: The measurement back-action only heats
the conjugate P̂ quadrature, which has no effect on themeasurement, as the oscillator dynamics does not couple it
to X̂.
Such back-action-evading (BAE) measurement schemes
allow one to measure one quadrature of a narrow-bandwidth
force acting on a mechanical oscillator with arbitrary precision.
They also naturally lead to squeezing of the measured
quadrature, albeit a conditional squeezing, where one musthave access to the full measurement record to produce
unconditional squeezing. Back-action-evading schemes arealso known as “quantum-nondemolition (QND) in time”
measurements, as one is making a QND measurement ofan observable which is explicitly time dependent. These
ideas have recently been discussed [6] within the specificcontext of quantum optomechanics [7], where a mechanicaloscillator is dispersively coupled to a driven cavity. A
BAE single-quadrature measurement was even implemented
experimentally by Herzberg et al., using a nanomechanicaloscillator coupled to a driven microwave-frequency cavity [8].
One could now ask about exploiting the idea of back-action
evasion in more general settings. The simplest generalization
would involve two harmonic oscillators, implying four quadratures of motion. Such a two-mode BAE scheme would involvemeasuring two quadratures of motion, with the back-action
only driving the unmeasured conjugate quadratures. If one
1050-2947/2013/87(6)/063846(24)could suitably couple the system, one could imagine measuring
both quadratures of an applied force without any quantum
limit.
That such a scheme is possible was recently discussed indetail by Tsang and Caves [9,10], though similar ideas werediscussed earlier (see, for example, Ref. [11] or Appendix D ofRef. [12]). From a practical perspective, it appears challenging,as it requires one of the two oscillators involved to have

(a)

(b)

(c)

FIG. 1. (Color online) (a) The system studied consists of two
mechanical oscillators, each coupled to a cavity (or circuit) mode.
Measurement is performed via phase-sensitive detection near thecavity resonance of the cavity output field. (b) Frequencies in thissystem, defined with respect to the cavity resonance frequency
ωc . The blue lines indicate the standard mechanical sidebands, at±ωa and ±ωb . The cavity is driven at ±ωm = ±(ωa + ωb )/2, as
indicated by vertical red lines on either side of the cavity resonance.
(c) Representation of the two-mode BAE scheme in the ideal, fully
symmetric case. The observable X̂+ is directly measured, and theobservable P̂+ is heated by the corresponding back-action. P̂−is dynamically coupled to X̂+ and therefore effectively measured
provided that   γ . In this limit, X̂− (which is dynamically coupled
to P̂− ) is also heated by the back-action of the measurement.

063846-1©2013 American Physical Society

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

effectively a negative mass (and hence a negative frequency).
Similar to standard single-oscillator back-action evasion, one
might expect that there is also a squeezing aspect to such
two-mode BAE. This was discussed by Hammerer et al. [13]:The squeezing here of the delocalized measured quadraturesnaturally gives rise to an EPR-style entanglement between thetwo modes.
In this paper, we now extend the above idea of two-modeback-action evasion to quantum optomechanics, specificallyconsidering a system where two mechanical oscillators arecoupled to a single electromagnetic mode, as depicted inFig. 1. Such two-mode optomechanical systems have recently
been realized experimentally [14]. We show that the two-mode
BAE scheme can be realized in a relatively simple fashion insuch a system, using a generalization of the two-tone driving
scheme used to achieve single-mode back-action evasion [1,6].
This is in contrast to previous proposals, where to achieve
the requisite negative mass, one of the two modes was notmechanical, but rather an atomic ensemble [13] or anothercavity mode [9]. Unlike the single-mode back-action-evasion
case, we expect that this scheme will not be susceptible toparametric instabilities [15], since here the number of photons
in the cavity will not oscillate at double the mechanicalresonance frequencies.
Note that applications of BAE measurement ideas to asystem of two mechanical oscillators have previously beenstudied [16–22], motivated by attempts at gravitational wave
detection. The approaches presented in those works differ
significantly from our scheme; in particular, they involvean explicit coupling between the mechanical oscillators,
something that is not required with our scheme.
We provide a thorough discussion of the optomechanicaltwo-mode BAE scheme, including a discussion of expected
experimental imperfections and methods for countering these.
In addition, we provide a full analysis of the force sensitivityof this scheme and the possibility of beating conventional
quantum limits. Unlike previous work [9,10], we include
dissipation in this analysis.
Finally, we also consider the generation of stationary,conditional entanglement using this scheme and the possibilityof turning this into unconditional entanglement using feedback. This is in contrast to more conventional measurementbased entanglement-generation schemes, which rely on strong
measurements and postselection; such schemes are wellestablished in optics [23–25] and have recently been discussed
in the context of superconducting circuits [26,27]. Our scheme
also differs from that in Ref. [13] as we generate stationary
entanglement and do not require a single-shot strong feedback
operation.
The remainder of this paper is structured as follows. In
Sec. II we derive the Hamiltonian and Heisenberg-Langevinequations describing this system and introduce a method
for compensating for system asymmetries. In Sec. III we
calculate the back-action heating effect and the measured noise
spectrum, with a view to assessing the 5 limits of our system.
In Sec. IV we rigorously compare the force sensitivityof our scheme, including imperfections, against conventional
quantum limits. The situation is significantly more complicated when one includes the presence of dissipation (unlike

Refs. [9] and [10]). We find that for a signal force appliedat the mechanical resonance frequency, the number of added
noise quanta can be made to go to zero (below the quantum
limit of 1/2),
n̄add [ωa ] → 0,

(1)even in the presence of dissipation. However, this is only possible in a narrow bandwidth about the mechanical resonanceand requires a careful tuning of optomechanical couplings.
For a signal force off-resonant from the mechanical resonancefrequency (detuned by an amount ), we find that the added
noise can be made to go to1n̄add [ωa + ] = √ ,
2 2

(2)

well below the standard quantum limit in this case (n̄add [ωa +] = /γ , where γ is the average mechanical damping rate)and even below the full quantum limit (i.e., allowing for
optimal detector noise correlations).
The conditional dynamics of the mechanical oscillators,
under a continuous measurement via the coupled cavity, aredescribed in Sec. V. The steady-state conditional variances are
determined, and these are used to determine in which regimes
the mechanical oscillators are conditionally entangled. The
presence of all-mechanical entanglement is determined from
Duan’s inequality [28], which here (in the case of symmetric
optomechanical couplings and mechanical damping rates)
becomes the requirement2(n̄th + 1/2)
< 1,
(3)
VX+ + VP− =ηCwhere VX+ (VP− ) is the variance in the sum (difference) of
the mechanical X̂ (P̂ ) quadratures, n̄th is the effective thermal
occupation of the mechanical environment [see Eq. (83a)], η isthe quantum efficiency of the measurement [see Eq. (73)], andC is the so-called optomechanical cooperativity, essentiallyquantifying the measurement strength [see Eq. (29)]. Bothquantities on the left-hand-side of Eq. (3) may be obtainedin a straightforward manner from the measurement record.
Therefore, the mechanical oscillators should be entangled
provided that the cooperativity obeys
C>2(n̄th + 1/2)
.
η

(4)

From an experimental perspective, this would appear to
be a modest requirement, well within the reach of current
optomechanics experiments based on microwave circuits[29]. In Sec. VI, the total unconditional variances of thecollective quadratures are calculated, and a feedback scheme
is introduced that allows the conditional entanglement to beconverted to unconditional entanglement.
II. SYSTEM
A. Hamiltonian

The system (see Fig. 1) is composed of two mechanicaloscillators, with resonance frequencies ωa and ωb , eachindependently and dispersively coupled (with strengths ga and063846-2

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

gb , respectively) to a common cavity mode having resonancefrequency ωc . The Hamiltonian isĤ = ωa â † â + ωb b̂† b̂ + ωc ĉ† ĉ + ga (â + â † )ĉ† ĉ
+ gb (b̂ + b̂† )ĉ† ĉ + Ĥdrive + Ĥdiss ,(5)

where â and b̂ denote mechanical lowering operators, ĉ denotes
the electromagnetic lowering operator, and Ĥdrive accounts fordriving of the electromagnetic mode. The term Ĥdiss accounts
for dissipation, with the modes subject to damping at rates γa ,γb , and κ, respectively.
To realize the two-mode BAE dynamics suggested by Tsang
and Caves [9], we need the two mechanical oscillators tohave equal and opposite frequencies and, moreover, havethe cavity (the detector in this scheme) only couple to asingle collective quadrature of the two mechanical oscillators(an operator which is explicitly time dependent). Both theserequirements can be realized by simply adapting the two-tone(i.e., ampltiude-modulated) driving scheme used for ordinary
single-quadrature BAE in Refs. [6] and [8]. One applies adrive to the cavity at both ωc + ωm and ωc − ωm , whereωm = (ωa + ωb )/2 is the average of the two mechanicalfrequencies; i.e.,Ĥdrive = (E+∗ e−iωm t + E−∗ eiωm t )eiωc t ĉ + H.c.

(6)

Working in an interaction picture with respect to Ĥ0 =ωm (â † â + b̂† b̂) + ωc ĉ† ĉ, one finds that the frequencies ofthe mechanical oscillators are now ± as desired, where2 = ωa − ωb . We further specialize to an optomechanicalsystem in the good cavity limit (ωm  κ) and assume large
driving amplitudes |E+ |,|E− |, such that one can linearize theoptomechanical interaction. In this standard regime, one finds(as desired) that in the interaction picture, the cavity only
couples to a collective mechanical quadrature with a timeindependent coupling. Specifically, defining the the quadrature
operators for each mode in the interaction picture in thestandard manner (α = a,b,c),√X̂α = (α̂ + α̂ † )/ 2,(7a)√†P̂α = −i(α̂ − α̂ )/ 2,(7b)
and collective (canonically conjugate) mechanical quadratureobservables via√
X̂± ≡ (X̂a ± X̂b )/ 2,(8a)√
P̂± ≡ (P̂a ± P̂b )/ 2,(8b)
the optomechanical Hamiltonian takes the form
Ĥ = (X̂+ X̂− + P̂+ P̂− ) + GX̂+ X̂c − Gd X̂− X̂c + Ĥdiss .
(9)Here, the many-photon optomechanical couplings are givenby√G ≡ 2 (ga + gb ) c̄,(10a)√
Gd ≡ 2 (gb − ga ) c̄,
(10b)

2
2
where c̄ = |E± |/ ωm + κ /4 is the steady-state amplitude of
the cavity field at the driven sidebands. Note that we want to
be in a regime where κ  , such that the cavity responds
fast enough to monitor the collective mechanical oscillation.

PHYSICAL REVIEW A 87, 063846 (2013)

The details of this derivation are given in Appendix A.
We stress here that the interaction terms in Eq. (9), evenaccounting for contributions from the off-resonant sidebandsof the driving fields, will not include terms oscillating at
double the mechanical resonance frequencies. Accordingly,
the cavity photon number will not oscillate at this frequency
and the parametric instabilities associated with single-modeback-action evasion [15] should be avoided in the two-modecase.
For equal optomechanical couplings (i.e., Gd = 0) andequal mechanical damping rates (γa = γb ≡ γ ), this systemperfectly realizes the BAE scheme of Tsang and Caves. To
see this, consider the Heisenberg-Langevin equations for oursystem in this ideal symmetric limit. Defining the vector ofcollective mechanical quadrature operators asV = (X̂+ ,P̂− ,X̂− ,P̂+ )T ,(11)

a standard calculation yieldsd 
V = M · V + FBA + N · ξ .
dt

(12)

Here, M describes the oscillator dynamics in the rotating
frame,
⎡
⎤
−γ /2

0
0
0
0 ⎥
⎢ − −γ /2
M=⎣
,
(13)
0
0
−γ /2
 ⎦
0
0
− −γ /2
while FBA is the back-action force from the cavity:
FBA = (0,0,0,−GX̂c )T .

(14)

Further,ξ = (X̂+,in ,P̂−,in ,X̂−,in ,P̂+,in )T ,N=√

γ I4 ,(15)

describe the Langevin noise arising from the dissipative bathsof each mechanical resonator, with In denoting the n × n
identity matrix. Finally, the equations of motion for the cavity
quadratures are˙ = − κ X̂ + √κ X̂ ,
X̂
(16a)
c
c
c,in
2
√
κ
P̂˙ c = −GX̂+ − P̂c + κ P̂c,in .
(16b)2
We thus have the desired behavior.
(i) The collective X̂+ and P̂− quadratures are linked dynamically the same way as X̂ and P̂ for a harmonic oscillator,
but are dynamically independent from their conjugate pair (P̂+and X̂− ). In the language of Tsang and Caves, the commuting
observables X̂+ and P̂− form a quantum-mechanics-freesubspace.
(ii) From Eq. (16b), we see that the cavity measures X̂+(and hence the dynamically linked Pˆ− quadrature) and its backaction only drive the unmeasured mechanical quadratures P̂+and X̂− [cf. Eq. (14)].
One thus has a setup which allows the BAE measurementof both of the collective quadratures, X̂+ and P̂− . Information
about the motion of these collective quadratures appears at thecavity output at sidebands ± from the cavity resonance.

063846-3

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

with the input noise vector given by

B. Compensating for unequal optomechanical couplings

While the scheme is easy to understand in the perfectly
symmetric case, any real expeirment will have to contend
with both an asymmetry in the optomechanical couplings(ga = gb ) and an asymmetry in the mechanical damping(γa = γb ). Asymmetric couplings cause Gd in Eq. (9) to benonzero, implying that the cavity measures both the X̂+ andX̂− collective quadratures. As a result, the measurement is nolonger BAE (i.e., the unwanted measurement of X̂− causes
back-action to drive P̂− , which then corrupts the measurementof X̂+ ). Even if there is no coupling asymmetry, if γa = γb ,the damping terms in the Heisenberg-Langevin equations
will link X̂+ and X̂− , also ruining the back-action-evasion.
Both these effects can easily be seen by writing out the
Heisenberg-Langevin equations in the general asymmetriccase; see Appendix B, especially Eqs. (B2) and (B5).
We show in later sections that while small deviations fromthe ideal symmetric system cause a departure from perfect
back-action evasion, one can still beat quantum limits on
detection as well as generate squeezing and entanglement.
Here, we point out that it is possible to exactly null thedeleterious effects of asymmetric optomechanical couplings
by introducing additional parametric drives on each mechanical resonator. Note that attempting to compensate forthe asymmetric optomechanical couplings by modifying theweights of the two cavity drives only leads to more unwanted
interaction terms.
To begin our analysis of the case ga = gb , it is useful
to introduce a new set of canonically conjugate collective
mechanical quadratures, such that in this basis, the cavity
only couples to a single (collective) quadrature. Defining therotation angle θ as
θ = tan−1(Gd /G) ,(17)

ξ = (X̃+,in ,P̃−,in ,X̃−,in ,P̃+,in )T .

The appropriate matrices are specified in Eqs. (B9) and (B10).
The cavity is now only coupled to a single collective
quadrature X̃+ , with a modified coupling constant G̃. The freeevolution of the oscillators in this quadrature basis is, however,
not of the desired form; the back-action evasion is ruined by
the terms in the last line of Eq. (19), which dynamically link themeasured observables (X̃+ ,P̃− ) to their (perturbed) conjugate
pair (P̃+ ,X̃− )It is now natural to ask whether it is possible to modify
the mechanical system in some simple way so as to eliminate the last line in Eq. (19). Such a modification would
require introducing both position and momentum couplings
between the two mechanical oscillators, something that is
not experimentally feasible. A more modest approach would
involve simply adding terms to the Hamiltonian of eachmechanical oscillator (but not introduce any coupling). We
consider parametrically driving each mechanical oscillator,
and thus adding a term Ĥpa to the system Hamiltonian of
the formĤpa = − (â 2 + â †2 + b̂2 + b̂†2 ).

(18a)

P̃± ≡ cos θ P̂± ∓ sin θ P̂∓ .

(18b)

The Hamiltonian of Eq. (9), in terms of the observables
introduced in Eqs. (18a) and (18b), takes the form
˜ X̃+ X̃− + P̃+ P̃− ) + G̃X̃+ X̂c + ĤdissĤ = (22
˜ X̃−+ 12 p(+ P̃−2 − X̃+− P̃+2 ),(19)

where we have definedp ≡ tan 2θ =ga2 − gb2
2Gd /G=,
22ga gb
1 − Gd /G2

(20a)

2
2
2
˜ ≡  cos 2θ =  1 − Gd /G = 2 ga gb , (20b)
1 + G2d /G2
ga2 + gb2

G̃ ≡ G sec θ = G

1 + G2d /G2 .

(20c)

One can recast the Heisenberg-Langevin equations describing
the system in terms of these rotated observables. The system
retains the form of Eq. (12), now in terms of the vector ofrotated observables,V = (X̃+ ,P̃− ,X̃− ,P̃+ )T ,

(21)

(23)

This is written in the interaction picture; in the nonrotating
laboratory frame, this corresponds to an equal-strength parametric modulation of each oscillator’s spring constant at thefrequency 2ωm . Writing this in the tilde quadrature basis ofinterest, we haveĤpa =

2
2(P̃−2 + P̃+2 − X̃+
− X̃−
).

(24)

˜
We thus see that by taking = p/2,
these terms combinewith the last term in Eq. (19) to give a Hamiltonian which has
the general BAE form

we introduceX̃± ≡ cos θ X̂± ∓ sin θ X̂∓ ,(22)

˜ X̃+ X̃− + P̃+ P̃− ) + G̃X̃+ X̂c
ĤC = Ĥ + Ĥpa = (2
˜ P̃−2 − X̃++ p() + Ĥdiss .
(25)
Note that the required value of the parametric modulation
corresponds to requiring the amplitude of the spring constant˜ imodulation of each oscillator to satisfy ki /ki = 4p/ω
(i, = a,b). In this Hamiltonian, the measured observables
X̃+ and P̃− are dynamically independent of their conjugate
quadratures, as required. The extra terms in the last lineof Eq. (25) only modify the dynamics of the unmeasured
subsystem. We thus see that two-mode BAE is possible
even without having perfectly symmetric optomechanicalcouplings.
III. NOISE SPECTRA

To gain a better appreciation of the two-mode BAE scheme
introduced above, we calculate the noise spectra of the measured and perturbed collective mechanical quadratures, and of
the homodyne current from measurement of the cavity output.
These may be calculated from the Heisenberg-Langevinequations in the usual manner [4]; the details are given in
Appendix C . Note that all quantum noise spectra quoted
throughout this paper are symmetrized quantum noise spectra.
Of particular interest here are the effects of asymmetries (either

063846-4

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

in the optomechanical couplings or in the mechanical damping
rates), which will cause deviations from the perfectly BAE
result. Readers not interested in the effects of imperfections(i.e., asymmetries) may wish to read Sec. III A and then jump
to Sec. IV, where we discuss beating the quantum limit on
force sensing.
A. Ideal symmetric caseWe start with the ideal case of symmetric optomechanicalcouplings (ga = gb ) and symmetric mechanical damping(γa = γb = γ ). The symmetrized noise spectral density of anobservable Ẑ ≡ [V ]i is defined by
SZ [ω] ≡

1
2

+∞
−∞

dt eiωt {Ẑ(t),Ẑ(0)} .

(26)

For the measured observable X̂+ it is
SX+ [ω] =1γ /2(n̄a + n̄b + 1)
22
(γ /2) + (ω − )2
γ /2,+(γ /2)2 + (ω + )21[n̄a + n̄b + 1 + C(2n̄c + 1)]
2
γ /2
γ /2
,
+
×
2
2
2
(γ /2) + (ω − )
(γ /2) + (ω + )2(28)

where n̄c is the number of thermal quanta associated withcavity dissipation, and we have quantified the optomechanicalcoupling G in terms of the cooperativity parameter [30],defined as
2G2
.
C=
γκ

noise spectrum SI [ω] as[ω] = SI [ω]/K2 ,(30)
SXmeas+√where K ≡ 2BG/ κ/2 is a measurement gain with B
denoting the amplitude of the local oscillator in the homodyne
detection scheme. Further details of the calculation of thisquantity are given in Appendix C. Taking the most interesting
limit κ  ,γ (i.e., the cavity is fast enough to measure thecollective quadrature), we find1/γ(2n̄c + 1).
(31)
8C
The second term here is the imprecision noise of the measurement (i.e., output shot noise). In the perfect BAE limit
considered here, it can be made arbitrarily small without any
resulting back-action penalty by increasing the cooperativity
via the driving strength E± . Note that the standard quantum
limit result at resonance (applicable when one is limited byback-action) is that the added noise must be at least as large
as the zero-point motion of the oscillator. We thus see that one
beats this standard quantum limit constraint when
SXmeas[ω] = SX+ [ω] ++C > (1 + 2n̄c )/8.

(27)

where n̄a and n̄b are the temperatures of the two mechanicalbaths (expressed as a number of thermal quanta). We are working here in the interaction picture as always (i.e., ω is measured
with respect to the average mechanical frequency ωm ) andhave assumed the good cavity (or “resolved-sideband”) limitκωm . The spectrum of Eq. (27) consists of two Lorentzians,
centered at ±, and weighted by the thermal and quantum
fluctuations of the mechanical oscillators. Since we are in
the perfect back-action-evasion limit, there is no dependence
on the optomechanical coupling G. Turning to the perturbed
observable P̂+ , the spectrum is given, in the limit κ  , by
SP+ [ω] =PHYSICAL REVIEW A 87, 063846 (2013)(32)

Typically, the thermal occupation of the cavity, n̄c , is small,such that the constraint of Eq. (32) is not experimentally
demanding. Indeed, henceforth, we set n̄c = 0. The results
accounting for nonzero n̄c are readily obtained by multiplying
cavity-dependent contributions to the noise spectra by(2n̄c + 1).
B. Asymmetric optomechanical couplings andmechanical damping

In the presence of coupling asymmetry, the measured
collective quadrature is now X̃+ (rather than X̂+ ) as definedin Eq. (18a), and so it is the noise spectrum of this quantitythat we seek. Our goal is to understand whether coupling anddamping asymmetries can be tolerated, or even exploited, forthe purpose of force sensing beyond the usual quantum limits.
In the presence of asymmetry, the two sets of collective
quadratures (X̃+ ,P̃− ) and (X̃− ,P̃+ ) are coupled to one another,increasing the complexity of the spectra. In general, we
write the noise spectral density as the sum of a contribution
from thermal and quantum fluctuations of the mechanicaloscillators, a contribution from the back-action of the cavity
detector, and an imprecision noise contribution, as
imp

SX̃meas[ω] = SX̃th+ [ω] + SX̃ba+ [ω] + SX̃ [ω].
+
+

(29)

The cooperativity here is essentially the measurement strength,
and thus describes the back-action heating rate of the perturbed
collective quadrature. It is equivalent to the quantity nBA used
in Ref. [6] to denote back-action heating of the perturbed
quadrature (as a number of quanta) in optomechanical singlemode back-action evasion.
We now turn to the spectrum of the cavity output. We
consider the case of a single-sided cavity, with homodyne
detection of the cavity output field. The measured spectrumof X̂+ is then defined in terms of the total homodyne current(33)

Now in the fully symmetric case considered previously, we
have SX̃th [ω] = SX+ [ω] as given in Eq. (27), SX̃ba [ω] = 0 (since++
the measurement is BAE), and SX̃meas [ω] is given by Eq. (31).
+
Now the simplest contribution to the spectrum in theasymmetric case is the imprecision noise contribution. Thisis given by [cf. Eq. (31)]
1/γ
,
8C̃
where γ is now the average mechanical damping rate,

063846-5

imp

SX̃ [ω] =(34)

γ = (γa + γb )/2,

(35)

+

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

and we have introduced the rotated cooperativity parameter by


2G̃2
2
2 arctan pC̃ == C[1 + (Gd /G) ] = C sec, (36)γκ
2which follows from Eqs. (20a), (20c), and (29).
The thermal and back-action contributions to the noise
spectrum exhibit complicated dependencies on system asymmetries. They are given, in a general form, by Eqs. (C6)and (C10), respectively. For weak asymmetries, however,
one finds that the measured noise spectrum of X̃+ still has
˜ as in theapproximately Lorentzian peaks at ω = ± (not ),
symmetric case. The system asymmetries will be quantified
through the dimensionless parameters p and d, which quantify
the degree of coupling and damping asymmetry, respectively.
The parameter p is defined in Eq. (20a), while d is defined byd = (γa − γb )/(γa + γb ).

(37)

Since the general expressions for these noise contributions
are complicated, we focus on two special cases: the noise at
the effective mechanical resonance (, where   γ ) and thenoise far detuned from the effective mechanical resonance( + , where     γ ). Furthermore, we focus on
the limit γ / → 0. First we consider the system in theabsence of the compensating parametric driving, such thatthe Hamiltonian is given by Eq. (9). In Sec. III B3 we consider
the compensated system, with Hamiltonian given by Eq. (25).
1. Noise at mechanical resonanceFirst we consider the noise contributions at the effective
mechanical resonance frequency. The noise contribution fromthermal and quantum fluctuations is a complicated functionof the coupling asymmetry, as this asymmetry modifies theeffective mechanical susceptibility. However, in the case where
the coupling asymmetry is zero (p = 0), we simply have
SXth+ [] =1
1(1/2 + n̄a ) +(1/2 + n̄b ) .
γa
γb

(38)

That is, the noise contribution from each mechanical oscillator
is simply given by the thermal and quantum fluctuations of thebath to which it is coupled, scaled by its damping rate. Anothercase of interest, for reasons that shall become clear later, is thatof matched asymmetries (p = −d), and in this case we find
SX̃th+ [] =1[(1/2 + n̄a )(1 + d 3 ) + (1/2 + n̄b )(1 − d 3 )].
γ(39)

Next we consider the noise contribution due to the backaction of the cavity on the mechanical oscillators. At theresonant peaks, this is
SX̃ba+ [] =1 (p + d)2 (1 + p2 )C̃.
γ (1 − d 2 + p2 )2(40)

That is, the asymmetries lead to back-action heating which
is proportional to the measurement strength C. As shown in
Fig. 2, in general, the back-action noise contribution increases
in the presence of either coupling or damping asymmetries (asindicated by the solid lines). Surprisingly, however, by having
appropriately tuned damping and coupling asymmetries suchthat p + d = 0, one can cancel this extra heating effect at

FIG. 2. (Color online) Back-action noise, SX̃ba+ [ω], for a range
of optomechanical coupling and mechanical damping asymmetries,
centered on ω = +. The plotted spectra are offset and normalized
by the noise spectral density at  in the ideal, fully symmetriccase; cf. Eq. (27). In the fully symmetric case (black line), there
is no back-action contribution to the spectrum. Asymmetries leadto back-action contributions to the spectra (solid lines). However, by
matching asymmetries the back-action contribution on resonance may
be nullified (dashed lines); see Eq. (40). The asymmetries for eachline are [cf. Eq. (37)]: d = 0,Gd /G = 0 (black); d = 0.05,Gd /G =0.025 (blue, solid); d = 0.05,Gd /G = −0.025 (blue, dashed); d =0.2,Gd /G = 0.1 (purple, solid); d = 0.2,Gd /G = −0.1 (purple,
dashed); d = 1,Gd /G = 0.414 (red, solid); d = 1,Gd /G = −0.414(red, dashed). The other parameters used here (and in subsequent
figures) are based on the experiments performed by Teufel andco-workers [29], as described in Sec. V G.

resonance (as indicated by the dashed lines). The measured
observable is driven by back-action from both the coupling
asymmetry and the damping asymmetry; at resonance theseare precisely out of phase and may coherently cancel. Unfortunately, though, this cancellation occurs only in a small
bandwidth about the mechanical resonance, as is described
quantitatively in Sec. IV. Recall that this back-action noise
cancellation occurs in the limit γ /  → 0. To higher-orderin γ / , and in the case of matched coupling and damping
asymmetries (p = −d), we have
1
SX̃ba+ [] =
γ



γ


2


1
(1 + d 2 )3 (1 + d 2 − 1 + d 2 )C. (41)8This result is useful in determining the ultimate limit to force
sensing in this system.
2. Noise away from mechanical resonance

Next we consider the noise spectrum far from the effective
resonance frequency (ω = ±) at a detuning , where    γ . The thermal contribution in this limit is1p(1 + d)1−+ n̄a
2
1 + p2



p
1
1+
. (42)+ n̄b+ (1 − d)2
1 + p2γ
SX̃th+ [ + ] =42

In the case of no asymmetries (p,d = 0), Eq. (42) reduces toa simple result consistent with Eq. (27).

063846-6

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

Next, the back-action contribution to the noise spectrum atω = ±( + ) is
SX̃ba+ [ + ] =

p2
γ
C̃.
42 1 + p2

(43)

Clearly, far from the effective mechanical resonance, theback-action heating due to the presence of damping asymmetryis strongly suppressed compared to the heating due to coupling
asymmetry. As we see in Sec. IV, this will allow for excellent
force sensing of a force applied to one of the mechanicaloscillators which is strongly detuned from resonance. In theabsence of coupling asymmetry, the dominant back-action
noise contribution is at higher order in γ /,


γ 2 2
γ
SXba+ [ + ] =
d C.
(44)
42 2
That is, the back-action heating due to damping asymmetryis suppressed by the small factor (γ /2)2 .
3. Noise with compensation

Now we turn to the noise spectra in the presence of thecompensation described in Eq. (25). The question is whether
such compensation is useful for the purpose of force sensing,and to answer this question, we calculate the measured noise
spectra. For weak asymmetries, the noise spectrum remains
well-approximated by Eq. (27) if we make the replacement
˜ That is, the effective collective mechanical oscillator → .
˜ (rather than ), but again we shall quote
frequency is now 
˜ assuming ˜  γ ) andresults in two cases: at resonance (,˜ + , assuming ˜    γ ).
far-detuned from resonance (At resonance, the thermal noise contribution is1
˜ =n̄a + n̄b + 1
SX̃th+ []γ [1 − d 2 (1 − p2 )]+ (n̄b − n̄a )d

1 − d 2 (1 + p2 ).
1 − d 2 (1 − p2 )(45)This reduces to the result of Eq. (38) in the case p = 0 (where˜ → ). The back-action contribution is˜ =SX̃ba+ []PHYSICAL REVIEW A 87, 063846 (2013)

contribution is a complicated function of the asymmetries. At
resonance it is possible to cancel this contribution by matching
asymmetries [see Eq. (40)], while the contribution far fromresonance is relatively insensitive to damping asymmetry [seeEq. (44)].
IV. FORCE SENSITIVITY

A key motivation for the two-mode BAE scheme is thepossibility of continuously monitoring both quadratures of anarrow-band force without any back-action-related quantum
limit. We now analyze this possibility in detail, paying close
attention to the role of mechanical dissipation, something thatwas not discussed in Ref. [10].
A. Conventional quantum limits

We consider the standard situation where a narrow-bandforce F (t) is applied to one of the two mechanical oscillators
in our setup (say the a mechanical oscillator):
ĤF = F (t)x̂a = f (t)(â + â † ),(47)

where f (t) = xa F (t). In the standard way, we express thetotal noise in the measured homodyne signal as an equivalent
added thermal noise on the driven oscillator. The measured
force noise spectral density has contributions both from theadded noise of the measurement as well as the inherent
quantum (zero-point) and thermal fluctuations of the measured
system, and takes the general form
SFmeas [ω] ≡ h̄mγa ω(1 + 2n̄a [ω] + 2n̄add [ω]),(48)

where n̄a [ω] describes the thermal occupation of the driven
oscillator and n̄add [ω] is the noise added by the measurement.
The most straightforward way to monitor F (t) would be tocontinuously monitor the position x̂a of the driven mechanicaloscillator. In this case, the full quantum limit on continuous
position detection [5] directly leads to a quantum limit on force
detection:n̄add [ω]  1/2.

(49)

2

d
1
C̃.
γ [1 − d 2 (1 − p2 )]2(46)

Clearly, the back-action heating due to coupling asymmetry
is now attenuated compared with the uncompensated case; cf.
Eq. (40).
Far from the effective mechanical resonance, the thermal
˜ + ], is given by Eq. (42) withnoise contribution, SX̃th [+˜ + ],
p = 0, and the back-action noise contribution, S ba [X̃+is given by Eq. (44) provided that we make the replacement
C → C̃. In both cases, the noise contribution due to coupling
asymmetry is suppressed by the use of compensation.

We refer to this as the full quantum limit on continuous force
detection: The added noise at each frequency is at best equal
to the zero-point noise.
Reaching the above-defined conventional quantum limit at
frequencies far from the mechanical resonance requires strong
correlations between the position detector’s back-action andimprecision noises [5], something that can be difficult toachieve. In the absence of such correlations, one is subject
to a more severe constraint, the so-called “standard quantum
limit.” For frequencies ω = ωa +  far from the mechanicalresonance (i.e., ||  γa ), this standard quantum limit takes
the form

C. Summaryn̄add [ωa + ]  /γa .

The noise spectrum of the measured observable may be
expressed as a sum of contributions from thermal and quantum
fluctuations, back-action noise, and imprecision noise, as per
Eq. (33). In the fully symmetric case there is no back-action
contribution, and the noise spectrum is given by Eqs. (27)and (31). In the asymmetric case, the back-action noise

By now calculating the added noise n̄add [ω] of the two-mode
BAE scheme which is the focus of this paper, we can determine whether one can surpass these conventional quantum
limits. Naively, one might think that as one is evading themeasurement back-action, there should be no quantum limiton the added noise. The situation is, however, more complex.

063846-7(50)

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

The auxiliary oscillator in our scheme (mechanical oscillatorb) will have its own quantum and thermal fluctuations which
contribute to the output noise and thus the added noise of themeasurement. One cannot simply set this to zero by making
γb = 0, because if γa = γb , the complete BAE nature of thescheme is lost (as shown above). Despite these complications,
we find that there are ways to perform force sensing using
our scheme that yield added noise numbers far below theconventional quantum limits.

This second contribution may be further decomposed into
imp
back-action noise, n̄ba
add [ω], and imprecision noise, n̄add [ω],contributions. These quantities follow from Eqs. (48) and (52).
While full expressions for the frequency-dependent added
noise n̄add [ω] may be readily derived, for simplicity we focus
on the two main cases of interest: that of a mechanically resonant signal force, and that of a far-detuned mechanical signal
force. Henceforth, in this section we consider results in thezero-temperature limit (n̄a ,n̄b = 0), allowing the possibilityof force sensing near and beyond quantum limits.

B. Results for two-mode BAE

The two-mode BAE scheme is capable of measuring both
quadratures of the signal force F (t) as long as it is contained ina bandwidth  about the mechanical resonance frequency
ωa ; we thus restrict attention to this case. The linearized
equations of motion readily yield the linear-response relation
between the measured observable X̃+ and the frequency
components of the applied force (see Appendix D). For
frequencies δ satisfying |δ|, we find

1. Detection of a mechanically resonant force

The resonant signal force corresponds to having F (t)centered in a narrow bandwidth γa about the mechanicalresonance frequency ω1 . We can consider resonant force
detection both without and with the compensation scheme
described in Sec. III B. From Eq. (51), we see that we needto know the transfer function χF at the effective resonancefrequency ω2 (δ = 0). One finds, in the regime ω2  γ , thatX̃+ [±(ω2 + δ)] = xa χF [±(ω2 + δ)]F [±(ω1 + δ)], (51)where ω2 =  and ω1 = ωa in the uncompensated case,and χF [ω] is a transfer function (which satisfies χF [−ω] =χF [ω]∗ ) whose form follows from the equations of motion;
we discuss its properties further in what follows. Note thatin the compensated case the relevant frequencies are shifted
˜ and ω1 = ωa = ωa + (1 − cos 2θ ), withsuch that ω2 = θ as defined in Eq. (17), and for this reason results in thissection are expressed in terms of the general frequencies ω1and ω2 . We consider both the cases of a mechanically resonant
signal force (δ = 0) and a mechanically nonresonant signal
force (δ =   γ ). We stress here that the left-hand side ofEq. (51) involves the frequency components of X̃+ [ω] in therotating frame used to write the system Hamiltonian in Eqs. (9)and (25), whereas the frequency components of F [ω] are inthe (nonrotating) laboratory frame.
Equation (51) tells us that the full spectral information of
the signal force is reproduced in the dynamics of the collective
mechanical quadrature X̃+ , the only proviso being that themechanical resonance frequency ω1 is effectively shifted toω2 . It follows that the measured force noise spectral density isgiven by
SFmeas [ω1 + δ] =h̄2 /(xa )2 measS[ω2 + δ],|χF [ω2 + δ]|2 X̃+(52)

where the calculation of the output noise SX̃meas [ω2 + δ] is+
discussed in Sec. III; recall that these noise spectra were
explicitly symmetric in frequency.
Writing the measured force noise spectral density in theform of Eq. (48), we can express the added noise of our schemeasiχF [±ω2 ] = ∓ gr (p,d),γ

where gr describes the modification of the susceptibility due
to asymmetries; its full form as a function of the asymmetry
parameters p and d is given by Eq. (D5) or (D6) in theuncompensated or compensated cases, respectively.
In the perfectly symmetric case (i.e., both mechanical oscillators have identical optomechanical couplings and damping
rates), gr = 1, implying that X̂+ responds to the applied force
analogously to an oscillator driven on resonance. In this case,from Eq. (31) we find that the added noise of the force-sensing
scheme on resonance isn̄add [ωa ] =

imp

(53)

where n̄aux
add [ω] is the added noise associated with the inherent
thermal and quantum fluctuations of the auxiliary mechanicaloscillator, and n̄cav
add [ω] is the added noise associated withthe coupling of the mechanical oscillators to the cavity.

1
1
1
+
→ .
2 8C
2

(55)

As expected, there is no back-action contribution to the added
noise, and the added noise is a monotonically decreasing
function of the measurement strength C. Nonetheless, one
cannot beat the standard quantum limit in this case. Even
though the measurement back-action goes to zero, the added
noise is still limited by the zero-point fluctuations of theauxiliary mechanical oscillator.
Now we consider force sensing in the presence of asymmetries. First we consider the case of damping asymmetry
without coupling asymmetry (d = 0,p = 0). In this case theadded noise due to the auxiliary mechanical oscillator is simplyn̄aux
add [ωa ] =cavn̄add [ω] = n̄aux
add [ω] + n̄add [ω]ba= n̄aux
add [ω] + n̄add [ω] + n̄add [ω],

(54)

11+d
.
21−d

(56)

This contribution may be understood in the following manner.
As the damping rate of the auxiliary mechanical oscillator isincreased beyond that of the driven oscillator (i.e., d from0 → −1), its noise spectrum is broadened. Therefore, most of
its noise is outside the force detection bandwidth of the driven
oscillator, and the added noise of the force detection scheme

063846-8

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

at resonance is reduced. If we consider the opposite limit, inwhich the damping rate of the auxiliary mechanical oscillator is
reduced below that of the driven oscillator (i.e., d from 0 → 1),
the added noise contribution is increased. The force noise
spectrum due to the auxiliary oscillator is narrowed, becoming
sharply peaked about resonance. In either case, however, thedamping asymmetry alone will lead to a back-action noise
contribution, as given by Eq. (40), that will typically dominate
the noise due to the auxiliary oscillator.
Next we consider the case of coupling asymmetry without
damping asymmetry (p = 0,d = 0). Again, coupling asymmetry leads to back-action heating and therefore reduces force
sensitivity. However, compensation, as per Eq. (25), can beused to almost restore the system’s force-sensing capability.
The added noise at the shifted mechanical resonance frequency, ωa , is now

arctan p111n̄add [ωa ] = sec2
+
→ ,
(57)
2
2
2
8C̃
with the limit (equal to the full quantum limit) being approached in the high-cooperativity, low-coupling asymmetry
regime. Again, the full quantum limit cannot be surpassed,
however.
The situation becomes more complex with both damping
and coupling asymmetries. In general, these asymmetries
cause back-action heating and reduce the system’s force
sensitivity. Indeed, using Eq. (40), the added noise due toback-action, to second-order in p and d, is
2n̄ba
add [ωa ] = (p + d) C.

(58)

The added noise now has terms which increase with themeasurement strength C.
However, from Eq. (40), we know that it is possible tocancel back-action on resonance by matching asymmetries(i.e., setting p = −d). This is promising for force sensing,but we must also check the added noise due to the auxiliary
mechanical oscillator. This contribution is


11−d
d
2
n̄aux
.
(59)
(1
+
d
+
d
[ω
]
=
)
1
+
√
a
add
21+d
1 + d2
Clearly, the added noise contribution of the auxiliary mechanical oscillator goes to zero when d → 1, and combined with
Eq. (40), this suggests the added noise of the scheme can go tozero with matched asymmetries. This suggestion is based on
calculations performed in the limit γ /  → 0. We can checkthis possibility by calculating the added noise to higher orderin γ /. The added noise contribution due to back-action withmatched asymmetries (p = −d) isn̄ba
1 + d 2 − d(1 − d − 1 − d 2 )]add [ωa ] = [1 − 
(1 + d 2 )3 1 γ 2
C.
(60)
×
1+d 8 From Eqs. (41) and (60) we find, in the extreme damping
asymmetry case (d = 1 and p = −1), that the added noisedue to the cavity isn̄cav
add [ωa ] =

γ2
3
.
C+
√
2
2
128( 2 − 1)C(61)

PHYSICAL REVIEW A 87, 063846 (2013)

While there is a back-action contribution here (first term), it is
proportional to (γ / )2 and hence small. One readily finds an
optimal measurement strength,
√

3
,
(62)
C0 = √8
2−1 γ
at which point the total added noise is√3γ√n̄add [ωa ]|C0 =
.
8
2−1

(63)

The added noise remains close to zero provided that   γ ,and so the full quantum limit can be beaten for resonant force
sensing via matched asymmetries.
Unfortunately, this result applies over only a very small
bandwidth about the mechanical resonance frequency. The
added noise may be evaluated as a function of frequency, andwe are interested in the bandwidth B over which the added
noise remains below the full quantum limit. It may be shown
that, in the case of extreme damping asymmetry (d → 1), thisbandwidth is

γ
1 γ
,
(64)
B= √ √ ∼γ
2 2 Cwith the second scaling attained for an optimal measurement
strength (cooperativity); see Eq. (62). That is, the usefulbandwidth becomes very small.
Plots of the added noise, as a function of the damping rateof the auxiliary oscillator (with the damping rate of the driven
oscillator held constant), are shown in Fig. 3. The added noise
is shown for three coupling asymmetries, and in each casethe added noise at resonance may be reduced to the level ofthe auxiliary oscillator fluctuations by matching coupling anddamping asymmetries.

FIG. 3. (Color online) Added noise for resonant (ω = ωa ) andnonresonant (ω = ωa + ) force sensing, n̄add [ω], as a functionof the damping rate of the auxiliary mechanical oscillator, γb /γa ,while γa is held fixed. The dashed black line is a lower boundon the contribution due to auxiliary oscillator thermal fluctuations;
see Eqs. (59) and (68a), for example. The blue, purple, and red
curves correspond to the added noise for resonant force sensing
at Gd /G = 0,−0.162 and −0.414, respectively. The back-action
contribution to the added noise can be made to go to zero on resonanceby matching asymmetries, irrespective of the size of the damping
asymmetry. The cyan line is the added noise for off-resonant force
sensing, with no coupling asymmetry and the measurement strength
optimized.

063846-9

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

2. Detection of a mechanically nonresonant force

We now turn to the case where the signal force F (t) is
detuned from the mechanical resonance, and again we consider
both the original and the compensated systems. Specifically,
we consider the case where F (t) is contained in a bandwidth
B centered on a frequency ω1 + , where ||  γa and B ||. From Eq. (51), we see that we need to know the transfer
function at ω2 + ,1
gn (p,d),χF [±(ω2 + )] =
2

+

(65)

where gn describes the modification of the transfer function
due to asymmetries; its full form as a function of the asymmetry
parameters p and d is given by Eq. (D7) in the original caseand by Eq. (D8) in the compensated case. As compared withthe case of a mechanically resonant signal force [see Eq. (54)],
the symmetric “gain” is now 1/(2) rather than 1/γ . Also notethat force detection, in this case, is not limited by a resonant
detection bandwidth of order γ .
In the perfectly symmetric case (i.e., both mechanical oscillators have identical optomechanical couplings and damping
rates), gn = 1. Then we find that the added noise is2 2 1
1
1→ ,(66)
n̄add [ωa + ] = +2γ
8C2with the limit being approached in the high-cooperativity
regime. Clearly, the scheme surpasses the standard quantum
limit on nonresonant force detection [cf. Eq. (50)]. It, however,is only equal to the full quantum limit [cf. Eq. (49)] thatcomes from allowing detector noise correlations. Again, theresidual added noise of 1/2 a quantum in Eq. (66) is due to thezero-point noise of the auxiliary mechanical oscillator.
We can easily determine the effects of asymmetries on ourforce-sensing scheme. From Eq. (43), to second-order in theasymmetry parameters and in the limit γ / → 0, we find2n̄ba
add [ωa + ] = p C.

These arguments can be made precise. The following results
describe the case with compensation included and so reduceto the case of no compensation when there is no coupling
asymmetry to begin with (p = 0 such that ωa → ωa ). The
added noise associated with the cavity is
γ 2 d2
4 arctan p
n̄cav
sec
C
[ω
+
]
=
add a
42 1 + d2

(67)

As already noted, in this highly detuned regime, it is only
the optomechanical coupling asymmetry that yields an appreciable deviation from being quantum-limited. This suggeststhat one could improve on the symmetric added noise resultof Eq. (66) by exploiting damping asymmetry to reduce theadded noise contribution from the auxiliary oscillator.
Indeed, these added noise contributions due to the auxiliary
mechanical oscillator are given by
11−d[1 + p(p + 1 + p2 )], (68a)21+darctan p
11−dauxn̄add [ωa + ] =
sec2
,
(68b)
21+d
2

n̄aux
add [ωa + ] =in the original and compensated cases, respectively. In bothcases, they go to zero in the limit d → 1. Combined with theinsensitivity of the back-action noise to damping asymmetry
as expressed in Eq. (67), this would appear to enable force
sensing beyond the full quantum limit, provided that there
is no coupling asymmetry. In fact, we might expect that thisforce sensitivity can be achieved in the presence of coupling
asymmetry if compensation is used.

1
42 1
.
γ 2 1 + d 8C

This added noise is minimized by a cooperativity,
arctan p
1
42
,
cos2
C0 = 2 √γ 2 2|d|2
at which point the added noise due to the cavity is


|d|
1

2 arctan p
n̄cav
sec
.
[ω
+
]|
=
√
C
0
add a
2
2 |1 + d|

(69)

(70)

(71)

Taking the limit d → 1, the added noise due to the auxiliary
oscillator vanishes as per Eq. (68b), and the total added noise
becomes
112 arctan pn̄add [ωa + ]|C0 = √ sec→ √ , (72)
2
2 2
2 2
with the result tending to the limit as the coupling asymmetry,p → 0. Now, one cannot only beat the standard quantum limit
for nonresonant force detection [cf. Eq. (50)] by a greater
amount, but one can also surpass the full quantum limit on
force detection [cf. Eq. (49)]. Without compensation, the added
noise in Eq. (72) grows more rapidly as |p| increases.
C. Summary

In summary, for sensing of a mechanically resonant force,
force sensing at the full quantum limit (being equal to thestandard quantum limit, on resonance) is possible in thesymmetric case. Even with coupling asymmetry, the possibilityof force sensing near this quantum limit is retained via
compensation. Further, it is possible to surpass this fullquantum limit via matched asymmetries, though only in avery narrow bandwidth.
For sensing of a mechanically nonresonant force we cansurpass the standard quantum limit in the symmetric case.
Furthermore, one may exploit damping asymmetry to surpassthe full quantum limit. In both cases, these results canbe retained in the presence of coupling asymmetry using
compensation. Additionally, force sensing in this case is not
limited to a resonant detection bandwidth.
V. CONDITIONAL VARIANCES AND ENTANGLEMENT

In Sec. III we calculated the noise spectra of measured
and perturbed observables in our system, with a view toanalyzing its force-sensing capabilities in Sec. IV. These are
all unconditional quantities; that is, the description of the stateof the system is not updated based on the measurements madeupon it. However, since we are continuously monitoring thesystem, we can also discuss its conditional dynamics, arisingfrom updating the best estimate of the state of the system
based on the measurement record. In particular, a continuous

063846-10

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

measurement conditionally projects the system into an eigenstate of the measured observable. This shall enable, in thissystem, conditional mechanical two-mode squeezing and theentanglement of the two mechanical oscillators.
A. All-mechanical entanglement

The preparation and verification of macroscopic, allmechanical entanglement is a fundamental goal of the studyof mechanical systems in the quantum regime [7]. Such astate is also a physical approximation of an Einstein-PodolskyRosen channel [31], a key ingredient in quantum information
processing protocols with continuous variables [32]. Electromechanical entanglement [33] and the entanglement ofphonons in bulk [34], at the single-phonon level, have bothbeen demonstrated, though this is not the case for mechanicalcontinuous variables. The entanglement of collective spin
operators of atomic ensembles [35] and that of the motional
states of trapped ions [36] have been achieved, however.
There has been a great deal of discussion of the possibilityof entangling macroscopic mechanical degrees of freedom
with other degrees of freedom (see [37] and references
therein), including other mechanical degrees of freedom.
Particular attention has been paid to the possibility of allmechanical entanglement with the interaction being mediatedby an electromagnetic field. Proposals studied include placing
mechanical oscillators in a ring cavity or interferometer, bothwithout [38–40] and with [41] measurement and feedback,placing dielectric membranes in a cavity [42], and using remote
optomechanical systems, both with [43–45] and without[46,47] protocols dependent on optical measurements.
More recently, Schmidt and co-workers have proposed an
all-mechanical entanglement generation scheme based on the
detuned, modulated (and therefore, two-tone) driving of acoupled electromagnetic cavity [48]. Compared with theirapproach, our scheme has the considerable advantages thatit does not require highly detuned driving of the cavity andis robust to the initial thermal populations of the mechanicaloscillators. Further, the possibility of entanglement verificationis built into our proposal.
B. Conditional dynamics in adiabatic limit

The conditional dynamics of the system may be describedin the standard manner using a stochastic master equation[49]. The evolution of the joint density operator, σ , of thetwo mechanical modes and one electromagnetic mode, under
homodyne detection of the cavity mode quadrature ĉeiφ +ĉ† e−iφ with a quantum efficiency η, is given bydσ = −i[Ĥ,σ ] dt + κ(n̄c + 1)D[ĉ]σ dt√+ κ n̄c D[ĉ† ]σ dt + ηκM[ĉeiφ ]σ dW

PHYSICAL REVIEW A 87, 063846 (2013)

We consider the “good measurement” limit, in which thecavity damping rate exceeds the rates in the Hamiltonians of˜ G̃ (˜ rather than  in
Eq. (19) or Eq. (25); that is, κ > ,p,
the compensated case). In this limit, the cavity responds rapidly
to the dynamics of the coupled mechanical oscillators, andwe may adiabatically eliminate the cavity mode [50]. Setting
φ = π/2 for√ convenience, the cavity lowering operator is givenby ĉ = −i 2G̃X̃+ /κ. Note that it is crucial that the adiabatic
elimination is performed in terms of the rotated observables;
otherwise, the dissipative terms in Eq. (73) will involve linear
combinations of mechanical annihilation operators, which
greatly complicates the description.
Now the evolution of the (reduced) joint density operatorof the two mechanical oscillators having traced out the cavity
mode, ρ = Trcav [σ ], is given bydρ = −i[Ĥ ,ρ] dt − [X̃+ ,[X̃+ ,ρ]] dt
2+ ηM[X̃+ ]ρ dW + γa (n̄a + 1)D[â]ρ dt+ γa n̄a D[â † ]ρ dt + γb (n̄b + 1)D[b̂]ρ dt+ γb n̄b D[b̂† ]ρ dt,

where Ĥ now denotes the Hamiltonian of either Eq. (19) or
Eq. (25), excluding terms involving the electromagnetic mode
operators. We have also introduced the collective quadrature
measurement rate as ≡ γ C̃ = 2G̃2 /κ,

where Ĥ denotes the original or compensated Hamiltonian ofEq. (19) or (25), respectively, D[Â]σ ≡ Âσ Â† − 12 Â† Âσ −1
σ Â† Â is the dissipative superoperator, M[Â]σ = Âσ +2σ Â† − Tr [Âσ + σ Â† ]σ is the measurement superoperator,and dW is the Wiener increment.

(75)

where G̃ is the modified coupling constant introduced in
Eq. (20c). The measurement record increment associated with
Eq. (74) is given by
dr = X̃+ dt + dW/ 4η.
(76)
Integrating Eq. (74) according to the measurement record
increment of Eq. (76) allows us to continually update ourbest estimate of the state of the two mechanical oscillators.
C. Best estimates of quadratures and conditional variances

Now according to Eq. (74), with the quadratic Hamiltonian
obtained from Eq. (19) or Eq. (25), and linear damping,decoherence, and measurement, the steady-state reduced
density operator ρ will be Gaussian. Therefore, the conditional
state of the system will be fully described by knowledge ofits first and second moments. One can immediately write
the equations of motion for the (scalar) best estimates of
the collective quadrature observables and for the conditionalvariances corresponding to the stochastic master equation in
Eq. (74). The vector of best estimates is formed asV̄ = (X̄+ ,P̄− ,X̄− ,P̄+ )T .

†+ γa (n̄a + 1)D[â]σ dt + γa n̄a D[â ]σ dt+ γb (n̄b + 1)D[b̂]σ dt + γb n̄b D[b̂† ]σ dt, (73)

(74)

(77)

For later convenience, we introduce the following notation for
the ten independent elements of the two-mode, symmetricallyordered covariance matrix, :⎤⎡
VX̃+ +− XX ++⎢ +− VP̃−− P P ⎥
−
⎥
(78)
=⎢
⎣ XX −− VX̃− −+ ⎦ .
++ P P −+ VP̃+063846-11

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

In general, the elements of  will describe covariances of
the rotated observables, though we drop the tilde notation in
the subscripts of the off-diagonal covariances in Eq. (78) for
convenience. In the symmetric case, the tilde notation on thesubscripts of the collective quadrature variances in Eq. (78)are dropped.
The equations for the best estimates of the collective (sumand difference) mechanical quadrature observables take theform of an Ornstein-Uhlenbeck process [cf. Eq. (12)],method for determining whether the mechanical oscillators areentangled is provided by Duan’s inseparability criterion [28].
Note that this is a sufficient, but not necessary, condition for
the inseparability of a bipartite quantum state. Specifically, we
consider a generalized version of Duan’s criterion [54], which
here takes the formd 
 · ξ (t),V̄ = M · V̄ + Q(79)
dt
where the system matrix M is given by Eq. (B9) or Eq. (B12),
for the Hamiltonian of Eq. (19) or Eq. (25), respectively, themeasurement noise weighting vector is = 4η(VX̃ ,+− ,XX ,++ )T ,
(80)
Q
+

In the case of symmetric optomechanical coupling (Gd = 0),
the criterion of Eq. (84) reduces to the standard version of
Duan’s criterion,and ξ (t) = dW/dt is a white-noise process describing thenoise in the measurement.
The equation for the best estimates of the covariances of
the collective mechanical quadrature observables is
˙ = M + MT + L − KT K,

VX̃+ + VP̃− < cos 2θ =1 − G2d /G2.
1 + G2d /G2

VX+ + VP− < 1.

E. The symmetric casewhere  is the symmetrically ordered two-mode covariance matrix in the ordered basis defined by Eq. (21),(X̃+ ,P̃− ,X̃− ,P̃+ ), and again M is given by Eq. (B9) or
Eq. (B12). The matrices K and L are defined by(82a)(K)11 = 4η, (K)ij = 0,
⎤
⎡ 

n̄d
0
n̄eq 0 ⎥⎢ 0 n̄eq 0
n̄
d ⎥
,
(82b)
L=γ⎢


⎣ n̄d
0 n̄eq 0 ⎦
0 n̄eq0n̄d

1. Conditional variancesn̄th = 12 (n̄a + n̄b ),(83a)n̄d = 12 (n̄a − n̄b ).

(83b)

Note that an equation of the form of Eq. (81) has previously
been obtained in a study of the optimal control of two-mode
entanglement (parametrically coupled optical modes) under
feedback [51].
Due to the rules of the Itō calculus, Eq. (81) is adeterministic system and we may solve for the steady-statecovariance matrix. The steady-state equation corresponding
to Eq. (81) is a continuous-time algebraic Riccati equation(CARE). This is a system of nonlinear algebraic equations,though one for which numerical methods are well-developed.
D. Steady-state conditional variances and entanglement

The two-mode, symmetrically ordered covariance matrix
 enables a full characterization of the entanglement andpurity of the system [52]. Indeed, the entanglement may bedirectly quantified using the logarithmic negativity [53]. Froman experimental perspective, however, the estimation of allelements of the covariance matrix is challenging. A simpler(85)

Crucially, both quantities on the left-hand side of Eq. (84) may
be obtained in a straightforward manner from the measurement
record. The quantity X̃+ is measured directly, while in theregime   γ , P̃− is dynamically coupled to it such that it isalso effectively measured.

(81)

where we have introduced the notation n̄eq = n̄th + 1/2 + d n̄dand n̄d = n̄d + d(n̄th + 1/2), in terms of the dimensionlessdamping asymmetry introduced in Eq. (37) and the effective
mechanical thermal occupations,(84)

In the ideal, symmetric case, an analytic form for the steadystate covariance matrix can be found (see Appendix E). Of
particular interest is the variance of the measured collective
quadrature. Taking the limit γ /  → 0 while the cooperativity
C is fixed leads to√√n̄th + 1/2 1
8ηC(n̄th + 1/2) + 1 − 1
VX+ =
→
√
√ ,
4ηC
2η
C
(86)
with the latter result following for C  1/[8η(n̄th + 1/2)].
This is a very good approximation for most experimentally
relevant scenarios, including that which we consider in moredetail below. Alternatively, taking the limit C → ∞ whileγ /  is fixed, we find the form√√n̄th + 1/2 1
4ηC(n̄th + 1/2) + 1 − 1
→
VX+ =
√ .
√
4ηC
2 η
C
(87)
The scaling with C is now the same as in the case  = 0 [cf.
Eq. (E3a)], which itself scales as in the case for single-mode
BAE measurement [6]. Clearly, for a sufficiently strong measurement (large cooperativity), the conditional variance of themeasured collective quadrature observable, X̂+ , is squeezed
below the vacuum level (1/2) and tends asymptotically to
zero.
To assess whether we are truly generating a two-mode
squeezed state and entanglement, we also need to check that theP̂− quadrature is squeezed by the measurement [cf. Eq. (85)].
It can be shown (see Appendix E) that 21γ[16η2 C 2 n̄tot
VP− = VX+ + 16ηC+ (1 + 4ηC n̄tot + 1 + 8ηC n̄tot )2 ](88)

063846-12

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

PHYSICAL REVIEW A 87, 063846 (2013)

to second order in γ / , where we have introduced the notation
for the total thermal and quantum fluctuations,n̄tot = n̄th + 1/2.

(89)

Note that Eq. (88) does not imply that the difference between
VP− and VX+ diverges as C → ∞. In this limit higher-ordercontributions are important, and the perturbative result of
Eq. (88) is not valid. However, provided the ratio γ /  issmall, from Eq. (88) we do have
VP− ∼ VX+ .

(90)

That is, in the limit that the collective oscillator frequency
greatly exceeds the mechanical damping rate, the collective
quadrature P̂− that is dynamically coupled to the directly
measured observable X̂+ is also effectively measured andtherefore conditionally squeezed. Accordingly, the measurement conditionally generates a two-mode squeezed state.
That both X̂+ and P̂− are effectively measured (andtherefore conditionally squeezed) in the regime   γ can beseen in the following way. Defining primed operators by the= e−i Ĥ0 t Ẑ± e+i Ĥ0 t with Ĥ0 = (X̂+ X̂− +transformation Ẑ±P̂+ P̂− ), we can express the observables X̂+ and P̂− as



X̂+
X̂+
cos tsin t=.
(91)− sin t cos tP̂−P̂−
While the observables X̂+ and P̂− oscillate at ±, theobservables X̂+and P̂− are constants of the motion withrespect to the closed system dynamics. From Eq. (91) it isclear that by continuously monitoring X̂+ , we are continuouslymonitoring both X̂+and P̂− . Since these observables also
determine P̂− , we are effectively measuring P̂− as well(provided that   γ ).
Recalling the “subspaces” (or subsystems) introduced inSec. II, we now know the conditional variances of theobservables in the “measured” subsystem. For the observablesin the “perturbed” subsystem, the variances in the limit
γ / → 0 follow from Eqs. (E9a) and (E9b) as
VX− ,VP+ = n̄th + 1/2 + C/2.

(92)

Effectively, we are measuring both X̂+ and P̂− in this limit,
and so the conjugate observables of both are equally perturbed.

FIG. 4. (Color online) The Duan quantity, VX+ + VP− , as a function of the cooperativity C, for effective temperatures corresponding
to n̄th = 0,1,25, and in the regime   γ . The bound on the Duancriterion, as in Eq. (85), is shown as the dashed black line. The
quantity VX+ + VP− is seen to be below this line for sufficiently
large cooperativities, indicating that the mechanical oscillators
are entangled. This remains the case even for a relatively large
thermal occupation. The required cooperativities are experimentally
accessible.

the total insensitivity to temperature that has been found in ascheme to generate an EPR channel in a cascaded atomic andnanomechanical system [13]. The primary distinction is that
here we are describing steady-state entanglement, rather thaninstantaneous entanglement arising from a strong single-shot
feedback operation.
The Duan quantity of Eq. (85) is plotted in Fig. 4, as afunction of cooperativity, for a range of thermal occupations
and for parameters corresponding to the experiment of Teufeland co-workers [29]. It is seen that the mechanical oscillators
will be entangled in an achievable parameter regime, even formechanical oscillators initially in a thermal state far from theground state.
Further, one can explicitly evaluate the entanglement
measure known as the logarithmic negativity [53]. Using theresult of Eq. (86) we find
EN = 12 {log2 [ηC/(n̄th + 1/2)] − 1}(94)

in the regime specified by the criterion of Eq. (93) andzero otherwise. The onset of entanglement, as a functionof cooperativity, according to this measure is found to beconsistent with that predicted by the application of Duan’scriterion [cf. Eq. (93)].

2. Entanglement

Next we explicitly consider the entanglement of themechanical oscillators. From Eqs. (86) and (90), in theregime   γ and the strong measurement limit, we find that
VX+ + VP− < 1, implying that our mechanical oscillators areconditionally entangled. Substituting the asymptotic form ofEq. (86) into Duan’s criterion leads to the sufficient conditionon the measurement strength for the generation of mechanicalentanglement by measurement,C>2(n̄th + 1/2)
.
η

(93)

Experimentally, this is not an overly demanding condition. It isinteresting to note, however, that this condition does not exhibit

F. The asymmetric case

If the damping rates or the optomechanical coupling rates ofthe mechanical oscillators are different, an analytical solution
for the steady-state covariance matrix is no longer possible
and we must resort to numerical methods. The one exception
to this is the case where we have coupling asymmetry
alone and compensation as per Eq. (25). In this case, theasymptotic results of Eqs. (86) and (87) remain valid, implying
entanglement can still be achieved for strong measurements;
see Appendix E for more details.
Figures 5 and 6 show the steady-state conditional variance
VP̃− and the generalized Duan quantity VX̃+ + VP̃− as afunction of the coupling asymmetry Gd /G. These are plotted

063846-13

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

FIG. 5. (Color online) The variances VP̃− (dashed lines), and thegeneralized Duan quantities VX̃+ + VP̃− (solid lines) as a functionof the optomechanical coupling asymmetry Gd /G for the original[Eq. (19)] and compensated [Eq. (25)] Hamiltonians. We assumed
γa = γb , though even for large damping asymmetries, the modification to these lines is small. The bound on the generalized Duaninequality, as in Eq. (84), is shown as the curved dashed blackline. This plot was prepared with the effective occupation set ton̄th = 0 and a cooperativity of C = 100. At such low temperature, thescheme is highly robust against coupling asymmetry. Even withoutcompensation, the mechanical oscillators remain entangled for acoupling asymmetry of up to ∼60%. With compensation, the scheme
is unaffected by coupling asymmetry.

assuming γa = γb , though the entanglement generated is only
very weakly dependent on damping asymmetry, even when thisasymmetry is large. Curves are shown for both the original
and the compensated cases in both figures, for the effective
occupation n̄th = 0 in Fig. 5 and for n̄th = 5 in Fig. 6. Notethat in the original (uncompensated) case, the variance of thedirectly measured observable VX̃+ (the difference between
the solid and dashed lines) is independent of the coupling
asymmetry, while the variance of the dynamically coupledobservable VP̃− is less effectively squeezed as the coupling
asymmetry is increased.

FIG. 6. (Color online) The variances plotted in Fig. 5, but now
with an effective thermal occupation n̄th = 5. At this higher effective
temperature, the robustness of the scheme (without compensation) is
reduced. The mechanical oscillators remain entangled for coupling
asymmetries up to ∼17%. With compensation, the scheme isunaffected by coupling asymmetry.

At n̄th = 0, the entanglement generation is seen to behighly robust against coupling asymmetry, even withoutcompensation. The mechanical oscillators remain entangled up
to a coupling asymmetry of ∼60%. At n̄th = 5, the robustness
against coupling asymmetry is reduced. In this case, theoscillators only remain entangled for a coupling asymmetry
up to ∼17%. It should, however, be possible to independently
cool the mechanical oscillators via auxiliary cavities such thatthey are initially at low thermal occupations [29]. Further, withcompensation included, the entanglement generation is totallyinsensitive to coupling asymmetry. Furthermore, the coupling
asymmetry for which entanglement generation is possible,
beyond a threshold value, is not strongly dependent on thecooperativity.
G. Experimental parameters

The experimental parameters used to prepare Figs. 2–7 correspond to a micromechanical membrane in a superconducting
microwave cavity, a design due to Teufel and co-workers[29]. One could conceivably form the microwave cavity fromtwo micromechanical membranes, forming the three-mode
optomechanical system that we have studied. However, thetwo-mode BAE scheme we have described could be performed
in a number of cavity optomechanics implementations.
For the purpose of numerical calculations, the cavity
resonance frequency is taken to be 8 GHz and the mechanicalresonance frequencies are centered around 10 MHz. The
cavity damping rate is κ/2π = 200 kHz, such that the systemis operating well into the resolved-sideband regime. The
mechanical damping rate is γ /2π = 100 Hz. The zero-point
fluctuations of the membrane are given by x = 4 fm,
such that the nominal single-photon optomechanical coupling
rate in Eq. (5) is gi /2π = 200 Hz. Further, suppose thecavity is driven such that its steady-state population at thedriven sidebands is 106 photons. This corresponds to a(nominal) effective optomechanical coupling rate in Eq. (9)that is G/2π = 200 kHz. These parameters correspond to acooperativity of C = 4 × 103 . This is a much larger cooperativity than is actually required for the scheme presented.
In fact, such an optomechanical coupling rate would invalidate the adiabatic elimination performed in order to obtain
Eq. (74).
For the purpose of calculations here, we take a nominal
effective coupling G/2π = 70.7 kHz, corresponding to acooperativity of C = 500. Throughout this work, a quantum
measurement efficiency of η = 1 has been assumed. Note
that, from the perspective of entanglement generation, one
can compensate for the reduced efficiency via a larger
cooperativity [cf. Eq. (93)]. These parameters are within thelimits of validity of our analysis: Recall that we assumed
resolved-sideband operation (ωa ,ωb  κ) in order to avoid
spurious back-action heating of the measured observable and˜ G̃) such that it
the “good measurement” regime (κ > ,p,
was possible to perform the adiabatic elimination.
H. Summary

The conditional dynamics of the system under homodyne
detection of the field output from the cavity are described by

063846-14

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

Eq. (73). In the adiabatic limit, in which the cavity responds
rapidly to the mechanical motion, we have the simplified
description of Eq. (74). This equation allows us to calculate thesteady-state conditional variances of the collective mechanicalquadratures [see Eq. (81)] and then assess the entanglementof the mechanical oscillators using Duan’s criterion [seeEq. (84)]. It is found that the mechanical oscillators will beentangled for reasonable experimental parameters (see Fig. 4),even in the presence of considerable asymmetries (see Figs. 5
and 6).
VI. UNCONDITIONAL VARIANCES AND FEEDBACK

The variances and covariances calculated in the previous
section as solutions of Eq. (81) are conditional quantities.
However, the measured output noise spectrum of Eq. (31) isan unconditional quantity. This means it consists of both thefluctuations described by the best estimate of the conditionalvariance (VX̃+ ) and the fluctuations in the best estimate of the2):
observable itself ( X̄+2.
VX̃tot+ = VX̃+ + X̄+(95)

We can define the total (unconditional) variances of the othercollective quadratures in a similar fashion. From Eq. (95),the fact that the conditional variance is squeezed below thevacuum level does not necessarily imply that the unconditionalvariance will be also. Here we calculate the total variances anddescribe how one can use feedback to reduce the unconditionalvariances to the conditional variances.
A. Unconditional variances

Fluctuations in the best estimates of the collective quadrature observables are described by Eq. (79). The steady-statefluctuations are then given by the solution of the Lyapunov
equation,¯ + M¯ T = −Q
Q
T ,
M

B. Measuring conditional variances

The conditional variances may be obtained from measurements of the unconditional variances and the fluctuations in
the best estimates, as per Eq. (95). The unconditional variances
may be obtained directly from the measured noise spectral
density, or equivalently, from the measurement record. The
unconditional variance of X̃+ is obtained directly, while
the unconditional variance of P̃− follows from looking at thequadratures of the measurement signal itself, as describedaround Eq. (91).
The fluctuations in the best estimates of X̃+ and P̃− maybe obtained from the filter of Eq. (79). Written out explicitly
in terms of the measurement record of Eq. (76), the filter isd V̄ = M · V̄ + 4η[I (t) − X̄+ ]Q,dt
X̄+ (t)cos t= 4η VX+ e−γ̃ t−sintP̄− (t)This is precisely as expected; a BAE measurement should
neither heat nor cool the mechanical oscillators in the unconditional sense. The unconditional variances in the asymmetriccase, both without and with compensation, are given in
Appendix F1 . Note, however, that Eq. (97) remains valid (interms of rotated observables) for the compensated case withno damping asymmetry.

(100)

C. Feedback
1. Unconditional variances

As noted earlier, feedback may be employed to reduce thefluctuations in the best estimates of our observables and so
reduce the unconditional variances to conditional variances.
In general, this may be achieved by adding the damping
terms −(αγ /2)X̄+ and −(αγ /2)P̄− to the equations for X̄+and P̄− (respectively) to the filter of Eq. (99). This requires
the application of (asymmetric) feedback forces to the two
mechanical oscillators, given by
αγ
Fa(b) (t) = √ (cos θ ∓ sin θ ) X̄+ sin ωm t2αγ
+ √ (cos θ ± sin θ ) P̄− cos ωm t,
2

(97)

(98)

∗ I (t),where the asterisk denotes convolution and we have introduced
the notation γ̃ = γ + 4ηVX+ . The fluctuations in the best
estimates may be obtained from the filtered measurement
record of Eq. (100).

Substituting the first asymptotic form of Eq. (86) into Eq. (97)
we find that
VPtot− = VXtot+ = n̄th + 1/2.

(99)

 and V̄ are as specified following Eq. (79), andwhere M, Q,
I (t) = dr/dt is called the measurement current. Of course,
the filter of Eq. (99) is itself dependent on knowledge of thesystem’s covariance matrix. We assume the system covariances
have reached their steady-state values: That is, one should use
the calculated values that follow from the steady-state solutionof Eq. (81).
In the fully symmetric case, the measured subsystemdecouples from the perturbed subsystem. Then the filter of
Eq. (99) can be recast in the limit +− → 0 as(96)

¯ is a matrix of second moments of best estimates
where of collective quadrature observables, described in the orderedbasis defined by Eq. (77), (X̄+ ,P̄− ,X̄− ,P̄+ )T . The matrix M isgiven by Eq. (13), Eq. (B9), or Eq. (B12), in the symmetric, uncompensated asymmetric or compensated asymmetric cases, is given by Eq. (15) in the symmetricrespectively, while Q
case and by Eq. (B3) in both asymmetric cases.
In the fully symmetric case, and in the relevant regime  γ , we find the total unconditional variances to be
VPtot− ∼ VXtot+ = VX+ (1 + 2ηCVX+ ).

PHYSICAL REVIEW A 87, 063846 (2013)(101)in the laboratory frame. From the form of Eq. (99), we
expect that X̄+ and P̄− will themselves oscillate at ±, andtherefore the applied feedback force of Eq. (101) will, in
fact, have significant weight at the two mechanical resonancefrequencies. In most cases (except the case where there are
significant coupling and damping asymmetries) it suffices toapply the feedback component proportional to the estimate X̄+alone.

063846-15

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

In order to calculate the variances under feedback we
can use an equivalent classical description of the dynamics
of our quantum system [55], coupled to the filter giving
the best estimates of the collective quadrature observables(and upon which the feedback will be based). The classicalquantity representing the collective mechanical quadratureX̃+ is denoted x̃+ (and similarly for the other collective
quadratures). The total unconditional variance of the measured
observable, under feedback, is then VX̃fb = x̃+2 .
+Now we can write the measurement record incrementin terms√ of such a classical representation as dr = x̃+ dt +dW/ 4η. We can also write the measurement record
increment in terms of the best estimateof the measured√
observable as dr = X̄+ dt + d W̃ / 4η. Comparing these
expressions leads tod W̃ = 4η(x̃+ − X̄+ )dt + dW.
(102)

2. The symmetric casePutting Eq. (102) together with the system dynamics andfilter from Eqs. (12) and (99), respectively, both the conditional
and unconditional dynamics can be described by the OrnsteinUhlenbeck process,

d 
dW
,
Y = −S · Y + T ·
dt
dt

(103)

where the system matrix S and the noise weighting matrix
T are defined in Eqs. (F2) and (F5) of Appendix F. The
 in Eq. (103)state vector Y and noise increment vector d W
are formed as = d W̄ ,Y = V̄ , d W
v
dw


(104)

where V̄ is the vector of best estimates of collective quadratureobservables introduced in Eq. (77), v is the vector of classicalrepresentations of the collective quadratures,v = (x̃+ ,p̃− ,x̃− ,p̃+ )T ,(105)

 is a vector of Wiener increments corresponding to thed W̄measurement noise, = dW (1,1,1,1)T ,
d W̄

(106)

and d w
 is a vector of independent Wiener increments describing the noise associated with the coupling of the oscillators totheir mechanical environments,dw
 = (dWx̃+ ,dWp̃− ,dWx̃− ,dWp̃+ )T .

(107)

Now Eq. (103) is, in general, an eight-dimensional
Ornstein-Uhlenbeck process, with the steady-state variances
given by solutions of the Lyapunov equation,
S + ST = TTT ,(108)

where  denotes the steady-state covariance matrix in theordered basis Y . Solving Eq. (108) provides all of theconditional and unconditional variances of interest. Note,however, that the matrices S and T themselves depend on
the steady-state conditional variances, which must first be
obtained by solving the algebraic Riccati equation giving thesteady-state of Eq. (81).

In the fully symmetric case, the equations for the measured
subsystem (X̂+ ,P̂− ) in Eq. (103) decouple from those for the
perturbed subsystem (X̂− ,P̂+ ). Since θ = 0 now, the feedback
force applied to both mechanical oscillators, from Eq. (101),
should beαγ
F (t) = √ X̄+ sin ωm t.
(109)
2
The equations for the measured subsystem alone can be
written in the same form as Eq. (103), with the appropriate
matrices specified in Appendix F 3 . This now describes afour-dimensional Ornstein-Uhlenbeck process, with steadystate solutions for the covariances again given by the solutionof Eq. (108). The reduced system size in this case facilitatesa perturbative solution. Expanding in the (assumed) small
parameter (1 + α)−1 , where α is the feedback gain in Eq. (109),
we find
VXfb+ = VX+ + 4ηC
VPfb− = VP− + 4ηC

VX2+
1+α
VX2+
1+α

,

(110a),(110b)

to first-order in (1 + α)−1 and in the regime   γ . That is, inthe limit of a large feedback gain α, the unconditional variances
reduce to the conditional variances. The conditional twomode squeezing and entanglement calculated in Sec. V cantherefore be converted to unconditional two-mode squeezing
and entanglement via feedback. In the opposite limit of
no feedback (α = 0), we have VXfb+ = VPfb− = n̄th + 1/2, as
expected.
Note that in the (opposite) case  = 0, the result of
Eq. (110a) is exact, consistent with the known result for
single-mode BAE measurement [6]. In this case we also have
VPfb− = VP− , as the observable P̂− is neither measured nor
perturbed by the measurement.
3. The asymmetric case

Again, in the asymmetric case all conditional and unconditional covariances may be obtained by solving the system
of Eq. (108). However, the size of the system is such that
we cannot obtain useful analytical solutions, and we resort tonumerical solutions. Figure 7 shows unconditional variances
under feedback as a function of the feedback gain, in theuncompensated asymmetric case and for an optomechanicalcoupling asymmetry of 5%. It is seen that these steady-state
unconditional variances (VX̃fb and VP̃fb ) approach the steady+−state conditional variances (VX̃+ and VP̃− ) asymptotically as
the feedback gain is increased. However, this approach isslower than in the fully symmetric case. In the compensatedasymmetric case, the unconditional variance VX̃fb is the same+as VP̃fb , and both tend to the (common) steady-state conditional−variances under strong feedback.
D. Summary

In this section the unconditional variances of the collective
mechanical quadratures were calculated. The fact that thereis two-mode squeezing in the conditional variances does not

063846-16

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

PHYSICAL REVIEW A 87, 063846 (2013)

all-mechanical entanglement and force sensing beyond conventional quantum limits.
Note added. Recently, we became aware of a recent related
work by Zhang et al. [60]. This work described how the
Hamiltonian of Eq. (9) in our paper could also be obtained
for two condensates in an optical lattice. However, neither ananalysis of the conditional dynamics under measurement nor
a calculation of the force sensitivity of the system was given
there.
ACKNOWLEDGMENTS
FIG. 7. (Color online) The steady-state unconditional variances
under feedback, VX̃fb+ and VP̃fb− , as a function of the feedback gain α[introduced in Eq. (101)], for an optomechanical coupling asymmetryof Gd /G = 0.05, n̄th = 25, C = 500, and in the regime   γ . The
steady-state conditional variances VX̃+ and VP̃− are shown as thedashed horizontal lines. These conditional variances are unequal dueto the coupling asymmetry. The steady-state unconditional variances
approach the steady-state conditional variances asymptotically as the
feedback gain is increased, though not as rapidly as in the case ofsymmetric coupling.

imply that there is two-mode squeezing in the unconditionalvariances. However, feedback can be employed [see Eq. (101)]
to reduce the unconditional variances to their corresponding
conditional variances (see Fig. 7).
VII. CONCLUSIONSA BAE measurement of a collective quadrature of two
mechanical oscillators may be performed using a cavity
detector. This is achieved, for a system operated in the
resolved-sideband regime, by two-tone driving of the cavity
with a detuning of plus and minus the average of the two
mechanical oscillator frequencies.
For the purpose of force sensing, one can surpass the fullquantum limit on resonant force sensing provided the damping
rate of the auxiliary oscillator is much lower than that ofthe driven oscillator, and the coupling asymmetry is matchedaccordingly. This would be experimentally challenging, so
perhaps more useful is the fact that one can significantly
surpass the standard quantum limit (and also surpass thefull quantum limit) on detuned force sensing without such
constraints on system asymmetries.
Further, in the adiabatic limit and in the regime where thecollective mechanical oscillator frequency is much greaterthan the average mechanical damping rate, one can conditionally prepare an entangled two-mode squeezed state of themechanical oscillators via measurement of the cavity output
field. The presence of this entanglement may be verified
in a straightforward manner from the measurement record.
Further, simple feedback based on the measurement record
may be used to convert the conditional two-mode squeezing
and entanglement to unconditional entanglement.
Given the rapid progress in the measurement of microwave
fields [8,56,57], including demonstrations of feedback control[58,59], the proposed scheme appears to be a promisingroute towards the generation and verification of macroscopic,

This work was supported by the DARPA ORCHID program
under a grant from the AFOSR. We thank Carlton Caves,
Gerard Milburn, Warwick Bowen, Ian Petersen, and Alexandre
Blais for valuable discussions.
APPENDIX A: DERIVATION OF HAMILTONIAN

We start from the Hamiltonian of Eq. (5), with the signal
force given by Eq. (6). Moving into an interaction picture withrespect to Ĥ0 = ωm (â † â + b̂† b̂) + ωc ĉ† ĉ, we can rewrite thisHamiltonian as
Ĥ = (â † â − b̂† b̂) + ga (âe−iωm t + â † eiωm t )ĉ† ĉ
+ gb (b̂e−iωm t + b̂† eiωm t )ĉ† ĉ + Ĥdiss + Ĥdrive .

(A1)

The environments of the three oscillators are assumed to
be ensembles of noninteracting oscillators. The usual Born,
Markov, and rotating-wave approximations are made on
the system-environment interactions, and we also ignore
environment-induced level shifts, as is typical for quantum
optical master equations [4]. Therefore, the Heisenberg equations corresponding to Eq. (A1), neglecting noise terms, areγaâ˙ = −iâ − iga eiωm t ĉ† ĉ − â,
(A2a)
2
γb
b̂˙ = ib̂ − igb eiωm t ĉ† ĉ − b̂,(A2b)2ĉ˙ = −iga (âe−iωm t + â † eiωm t )ĉ − igb (b̂e−iωm t + b̂† eiωm t )ĉκ(A2c)− iE+ e−iωm t − iE− eiωm t − ĉ.
2Adopting the ansatz ĉ(t) = ĉ0 (t) + ĉ+ (t)e−iωm t + ĉ− (t)eiωm t ,assuming ωm  κ (resolved-sideband regime), explicitly separating out the Fourier components of the field at the driven
sidebands, and then equating frequency components [61], we
obtain the systemγa††â˙ = −iâ − iga (ĉ0 ĉ+ + ĉ− ĉ0 ) − â,
(A3a)
2
γb
†
†
b̂˙ = ib̂ − igb (ĉ0 ĉ+ + ĉ− ĉ0 ) − b̂,
(A3b)
2
κ
ĉ˙ 0 = −iga (â ĉ− + â † ĉ+ ) − igb (b̂ĉ− + b̂† ĉ+ ) − ĉ0 , (A3c)2
˙ĉ+ = −iE+ + iωm ĉ+ − κ ĉ+ − iga â ĉ0 − igb b̂ĉ0 ,
(A3d)
2
κ
ĉ˙ − = −iE− − iωm ĉ− − ĉ− − iga â † ĉ0 − igb b̂† ĉ0 . (A3e)2Solving Eqs. (A3d) and (A3e) for the steady-state at thedriven sidebands, assuming the optomechanical couplings are

063846-17

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

relatively small, yields−iE+≡ c̄+ ,
−iωm + κ/2−iE−≡ c̄− .
=iωm + κ/2

rotated bases and both without and with the compensating
parametric driving.

ĉ+ =

(A4a)

ĉ−

(A4b)

1. Original system, original basis

Then we can write down the effective Hamiltonian, corresponding to Eqs. (A3a)–(A3c), asĤ = (â † â − b̂† b̂) + Ĥdiss∗∗ĉâ † + c̄+ĉâ + c̄− ĉ† â)+ ga (c̄+ ĉ† â † + c̄−∗∗ĉb̂† + c̄+ĉb̂ + c̄− ĉ† b̂).
+ gb (c̄+ ĉ† b̂† + c̄−(A5)

Now we assume that the cavity sideband amplitudes will havethe same steady-state amplitudes but different phases, c̄+ =c̄e−iψ and c̄− = c̄eiψ . Here c̄ is assumed to be real, without
loss of generality, because it only changes the cavity quadrature
to which the collective mechanical quadrature is coupled, andis readily compensated by adjusting a local oscillator phase.
The Hamiltonian becomesĤ = (â † â − b̂† b̂) + 2ga c̄X̂a,ψ X̂c + 2gb c̄X̂b,ψ X̂c + Ĥdiss ,(A6)
where we have introduced quadratures of mechanical√ andelectromagnetic modes as X̂√q,ψ ≡ (q̂eiψ + q̂ † e−iψ )/ 2 andP̂q,ψ ≡ −i(q̂eiψ − q̂ † e−iψ )/ 2, with the lack of a phase
subscript implying that the phase has been set to zero. The
Hamiltonian of Eq. (A6) may be rewritten as
Ĥ = (X̂+ X̂− + P̂+ P̂− ) + G(cos ψ X̂+ − sin ψ P̂+ )X̂c
− Gd (cos ψ X̂− − sin ψ P̂− )X̂c + Ĥdiss .

where the average mechanical damping rate γ and the dimensionless damping asymmetry d are as introduced in Eqs. (35)and (37), respectively. Clearly from Eq. (B2), the asymmetricdamping directly couples the nominally “measured” (X̂+ ,P̂− )and the nominally “perturbed” (X̂− ,P̂+ ) subsystems. The
vector of noise input operators is given by Eq. (15), though thenoise input weighting matrix is⎤⎡√√
γ+ (d) √ 0γ− (d) √ 0⎢ 0
γ+ (d) √ 0γ− (d) ⎥⎥ , (B3)√N=⎢⎦⎣ γ− (d)
0
γ
(d)
0
+
√
√
γ− (d)0γ+ (d)
0where we have introduced the collective damping rate notation,√√ √γ± (d) ≡ γ ( 1 + d ± 1 − d)/2.
(B4)
The vector of back-action forces associated with the coupling
to the cavity is nowFBA = (0,Gd X̂c ,0,−GX̂c )T ,(A7)Irrespective of the value of ψ, there exists back-action in the“unperturbed” subspace due to asymmetric coupling (Gd =0), and so we may set ψ to zero. Then the Hamiltonian reduces
to that given in Eq. (9).
APPENDIX B: HEISENBERG-LANGEVIN EQUATIONSThe system of Heisenberg-Langevin equations describing
the dynamics of the two mechanical oscillators is given inSec. II by Eq. (12), repeated here for convenience,d 
V = M · V + FBA + N · ξ ,dt

In the presence of both coupling and damping asymmetries,
in the basis of the original collective observables of Eq. (11)and without compensation, the system matrix is
⎡
⎤
−γ /2

−dγ /20−γ /20−dγ /2 ⎥
⎢ −
M=⎣
,
(B2)
−dγ /2
0
−γ /2
 ⎦
0
−dγ /2
−
−γ /2

(B1)

where V is the vector of collective mechanical quadrature
operators of Eq. (11) and ξ is the vector of input noise operatorsof Eq. (15). The coupling of the mechanical oscillators to thecavity mode is included in Eq. (B1) as an inhomogeneity,
the vector of back-action forces FBA . In the fully symmetriccase (symmetric mechanical damping rates and symmetric
optomechanical coupling rates) the system matrix M is given
by Eq. (13), the back-action force vector FBA is given by
Eq. (14), and the input noise weighting matrix N is given by
Eq. (15). The Heisenberg-Langevin equations for the cavity
quadratures themselves are given by Eqs. (16a) and (16b).
Here we give the Heisenberg-Langevin equations describing
the system in the presence of both damping and coupling
asymmetries. These equations retain the form of Eq. (B1) inall cases, and they shall be specified both in the original and(B5)with the asymmetric coupling rate Gd as introduced in
Eq. (10b). The equation for the cavity quadrature X̂c is still
Eq. (16a), while the equation for P̂c is now√κP̂˙ c = −GX̂+ + Gd X̂− − P̂c + κ P̂c,in .
(B6)2
According to Eq. (B6), the cavity is no longer coupled to (i.e.,measures) X̂+ alone, but to a linear combination of X̂+ andX̂− . Accordingly, the observable P̂− (which is dynamicallycoupled to the measured observable X̂+ ) is now heated byback-action from coupling to the cavity, as is also clear from
Eq. (B5). Clearly, both types of asymmetries ruin the perfect
BAE property of the measurement scheme. The significanceof these asymmetries is addressed throughout the paper.
2. Original system, rotated basis

In the presence of optomechanical coupling asymmetry
it is convenient to move into a basis of rotated collective
quadrature observables, as defined in Eqs. (18a) and (18b). The
Heisenberg-Langevin equations retain the form of Eq. (B1),with V now denoting the vector of rotated quadrature observables as defined in Eq. (21),V = (X̃+ ,P̃− ,X̃− ,P̃+ )T ,(B7)

and the vector of input noise operators as introduced in
Eq. (22),063846-18

ξ = (X̃+,in ,P̃−,in ,X̃−,in ,P̃+,in )T ,(B8)

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

with the rotated collective input noise operators defined interms of the original collective input noise operators in theobvious manner.
The system matrix is now [cf. Eq. (B2)]
⎤
⎡
˜
˜
−γ /2

−dγ /2 −p
⎢ −
˜
˜
−γ /2
−p
−dγ /2 ⎥
⎥
(B9)
M=⎢
˜
˜ ⎦,
⎣ −dγ /2
p
−γ /2

˜
˜
p
−dγ /2−−γ /2
˜ as introduced in Eq. (20b), while the matrix N is stillwith 
given by Eq. (B3). The back-action force vector is nowFBA = (0,0,0,−G̃X̂c )T ,(B10)

while the Heisenberg-Langevin equations for the cavity
quadratures are now Eqs. (16a) and [cf. Eq. (B6)]√κP̂˙ c = −G̃X̃+ − P̂c + κ P̂c,in .
(B11)
2
From Eq. (B9), it is clear that in the rotated basis both
the damping and the coupling asymmetries directly couple
the nominally measured (X̃+ ,P̃− ) and nominally perturbed(X̃− ,P̃+ ) subsystems. However, now only the rotated observable X̃+ is coupled to the cavity (i.e., measured); indeed, thismotivated the definition of the rotated observables in the firstplace. Accordingly, and as seen in Eq. (B10), the back-actionof the cavity only directly heats P̃+ .
3. Compensated system, rotated basis

Including compensating parametric driving as per Eq. (25),the system matrix in the rotated basis of Eq. (21) is
⎤
⎡
˜
−γ /2

−dγ /2
0
⎢ −
˜
−γ /2
0
−dγ /2 ⎥⎥M=⎢˜˜ ⎦ . (B12)⎣ −dγ /2
2p
−γ /2

˜
˜
2p
−dγ /2−−γ /2The effect of the compensation is that the nominally measured
subsystem is no longer coupled to the nominally perturbed
subsystem via the coupling asymmetry, though it remains
coupled through the damping asymmetry. This partial decoupling shall prove to be useful. The other parts of Eq. (B1)remain unchanged: The vector of noise input operators isgiven by Eq. (22), the vector of back-action forces is given
by Eq. (B10), and the equations for the cavity quadratures are
given by Eqs. (16a) and (B11).
APPENDIX C: CALCULATION OF SPECTRAPHYSICAL REVIEW A 87, 063846 (2013)

For convenience, we also introduce the thermal susceptibility
matrix,χ̄ [ω] ≡ χ [ω] · N.

(C3)

The symmetrized noise spectral densities of interest, for an
observable Ẑ = [V ]i , are defined in Eq. (26). They may beobtained from frequency-domain solutions of the HeisenbergLangevin equations as
SZ [ω] = (Ẑ[ω])† Ẑ[ω].

(C4)

From Eq. (C1), there are clearly two contributions to the noise
spectra of the collective mechanical quadratures. The first term
describes the intrinsic thermal and quantum fluctuations of themechanical oscillators, while the second term describes the
back-action heating of the mechanical oscillators due to theircoupling to the cavity. We now calculate these contributions
in turn.
2. Thermal noise spectra

The thermal noise spectra ultimately depend on the correlation functions of the input mechanical noise operators; theseareẐi,in (t)Ẑi,in (t  ) = 12 (n̄a + n̄b + 1)δ(t − t  ),Ẑi,in (t)Ẑj,in (t )= 12 (n̄a − n̄b )δ(t − t  ),

(C5a)
(C5b)

for i = j in the second correlation function, and whereẐ ∈ {X̃,P̃ } and i,j ∈ {+,−}. Note that n̄a (n̄b ) is the thermal
occupation of the environment to which mechanical oscillatora(b) is coupled. Therefore, the thermal contribution to thenoise spectral density of the measured observable X̃+ is

 4

1
th
2
|χ̄1m [ω]|SX̃+ [ω] = (n̄a + n̄b + 1)
2
m=1
 2

1
∗
(χ̄1m
[ω]χ̄1(m+2) [ω]+ (n̄a − n̄b )
2
m=1

∗
+ χ̄1(m+2)[ω]χ̄1m [ω]) ,(C6)

with χ̄j k [ω] denoting elements of the thermal susceptibility
matrix introduced in Eq. (C3). This quantity is evaluatedexplicitly, both at the resonant peaks ( in the original case,˜ in the compensated case) and at a large detuning , in
Sec. III B.

1. Heisenberg-Langevin equations in the frequency domain

The back-action of the cavity on the mechanical oscillators
may be determined by calculating noise spectra of themeasured and perturbed observables. The system of Eq. (B1)
is readily solved in the frequency domain in all cases. Taking
the Fourier transform of Eq. (B1) leads toV [ω] = −χ [ω] · N · ξ [ω] − χ [ω] · FBA [ω],(C1)

where we have introduced the cavity susceptibility,

where we have introduced the susceptibility matrix,
χ [ω] = (M + iω1)−1 .

3. Back-action noise spectra

The contribution to the noise spectral density of themeasured observable due to back-action is itself dependent
on the dynamics of the cavity. Accordingly, we first take theFourier transform of the coupled cavity quadrature in Eq. (16a),√X̂c [ω] = χc [ω] κ X̂c,in [ω],
(C7)

(C2)
063846-19

χc−1 [ω] ≡ −iω + κ/2.

(C8)

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

The input noise correlation function for this cavity quadratureisX̂c,in (t)X̂c,in (t  ) = 12 (2n̄c + 1) δ(t − t  ),(C9)

where n̄c is the thermal occupation of the cavity bath.
Frequently, in experiments, n̄c will be very close to zero. From
Eqs. (B10), (C1), (C7), and (C9), the back-action contribution
to the spectrum of the measured observable, assuming thatκ  , isSX̃ba+ [ω] = γ |χ14 [ω]|2 C̃(2n̄c + 1),(C10)

where χ14 [ω] is an element of the susceptibility matrix
introduced in Eq. (C2) and C̃ is the rotated cooperativity
parameter introduced in Eq. (36). Equation (C10) is evaluatedexplicitly at resonance and far-detuned from resonance inSec. III B.
However, for the purpose of assessing the system’s usefulness for force sensing, we need to know the bandwidth
over which quantum limits can be surpassed. This requires
knowledge of the full frequency-dependent contribution to
the noise spectrum. In general, this is complicated, thoughreasonably simple forms may be obtained in the perturbative
regime (i.e., to second-order in p and d) and in the extreme
asymmetric case (d = 1 and p = −1 without compensation,or just d = 1 in the case with compensation).
4. Output noise spectrum

The total output noise spectrum of the measured collective
mechanical quadrature consists not only of the quantum andthermal fluctuations of the mechanical oscillators (calculated
in Appendix C) and the back-action heating due to thecoupling to the cavity (calculated in Appendix C), but also
includes the noise added by the detector (i.e., the cavity).
The total noise spectrum is readily calculated using the
input-output formalism of quantum optics [4], and this will
facilitate comparisons with conventional quantum limits on
measurement [5].
To calculate this noise spectrum, we need to know how
the collective mechanical quadrature couples to the cavity
field. This is given, in the frequency domain, by the Fouriertransform of Eq. (B11),√P̂c [ω] = −χc [ω]G̃X̃+ [ω] + χc [ω] κ P̂c,in [ω],(C11)
where the cavity susceptibility χc [ω] is defined in Eq. (C8).
From Eqs. (C7) and (C11), we have for the cavity modeannihilation operator√√ĉ[ω] = κχc [ω]ĉin [ω] − iχc [ω]G̃X̃+ [ω]/ 2.
(C12)
Applying the usual boundary condition for a single-sided
optical cavity,√
κ ĉ[ω] = ĉin [ω] + ĉout [ω],(C13)
we find the cavity output field to bei(ω − ωc ) + κ/2
ĉout [ω] = −c̄out [ω] −ĉin [ω]i(ω − ωc ) − κ/2−iχc (ω − ωc ) κ/2G̃X̃+ [ω − ωc ], (C14)with the frequencies now specified in the (nonrotating)laboratory frame. Here c̄out [ω] describes the output field dueto the coherent cavity driving fields, the second term describes
fluctuations of the input field filtered by the cavity, and thethird term carries the signal due to the mechanical oscillation.
It is assumed that the output cavity field is subject tohomodyne detection in the usual manner [4]. The measured
homodyne current is then†∗(t)ĉout (t) + bLO (t)ĉout (t),Iˆ(t) = bLO(C15)

where the local oscillator amplitude is bLO (t) = iBe−iωc t ,with B assumed real without loss of generality. The Fouriertransform of each component in Eq. (C15) is evaluated asa convolution integral, allowing us to calculate the spectrumof the homodyne measurement current SI [ω], as quoted in
Eqs. (30) and (31) [6].
APPENDIX D: FORCE-SENSING TRANSFER FUNCTIONS

The ability to perform a BAE measurement of an oscillating mechanical observable suggests the possibility of force
sensing beyond conventional quantum limits. Given that thenoise spectral density of the measured observable has been
calculated (see Sec. III), the next task is to calculate the transfer
function relating the signal force on one mechanical oscillator
to the measured collective mechanical quadrature. This allows
one to determine the noise added by the force-sensing schemeand facilitates a comparison with conventional quantum limits(see Sec. IV). These transfer functions are calculated here.
Suppose that the mechanical oscillator a is subject to thesignal force f (t), as per Eq. (47). In a frame rotating at the
average mechanical frequency ωm [i.e., the same frame used towrite the Hamiltonian in Eq. (9)], the signal force is describedbyĤF = f¯(t)â + f¯∗ (t)â † ,(D1)f¯(t) = f (t)e−iωa t eit .

(D2)where we have

The Hamiltonian driving of Eq. (D1) adds a vector of driving
forces F (t) to the Heisenberg-Langevin equations of Eq. (12),given by⎡⎤(cos θ − sin θ ) Im f¯(t)⎢ (cos θ + sin θ ) Re f¯(t) ⎥F (t) = − ⎣,(D3)(cos θ + sin θ ) Im f¯(t) ⎦(cos θ − sin θ ) Re f¯(t)where θ is the angle describing the coupling asymmetry
introduced in Eq. (17). In order to evaluate the force-sensing
transfer function, we can neglect the noise terms and thecoupling to the cavity in Eq. (12), leaving the system ofHeisenberg equationsd 
V = M · V + F (t),
dt

(D4)

where V is the vector of rotated quadrature observables in
Eq. (21), and the system matrix M is given by Eq. (B9)or Eq. (B12) for the uncompensated or compensated cases,
respectively. The system of Eq. (D4) with the drive of Eq. (D3)
is readily solved in the frequency domain, leading to therequired transfer functions. These results are now given, for thecases of a signal force resonant with the mechanical oscillator

063846-20

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

and for a signal force far-detuned from the mechanicalresonance frequency.

PHYSICAL REVIEW A 87, 063846 (2013)
APPENDIX E: BEST ESTIMATES OF QUADRATURES AND
CONDITIONAL VARIANCES—THE SYMMETRIC CASE
1. Quadratures

1. Detection of mechanically resonant force

First we consider a mechanically resonant signal force.
For the uncompensated system, assume the signal force f (t)is contained in a narrow bandwidth about the mechanicalresonance frequency ωa . This implies that in the rotating frame
f¯[ω] is peaked at ±. For the compensated system we assume
f (t) is contained in a narrow bandwidth about ωa , leading
˜ The transfer function between
to f¯[ω] being peaked at ±.
the signal force and the measured observable takes the form
specified by Eqs. (51) and (54).
The transfer functions are modified in the presence ofasymmetries by the function gr (p,d). For the uncompensatedsystem, this modification is given by

In the fully symmetric case (equal optomechanical coupling
rates and mechanical damping rates), the measured subsystem(X̂+ ,P̂− ) decouples from the perturbed subsystem (X̂− ,P̂+ ) in
Eqs. (99) and (B1). For the measured subsystem alone, we canstill write the filter, for instance, in the form of Eq. (79), but
now with the matrices
−γ /2

,
− −γ /2

 = 4η VX+ ,
Q
+−

M=

1 + p2 (1 + d + p) − 1 − d − dp − p22(1 − d 2 + p2 )arctan p.
(D5)
× cosec2

This factor is just 1 for vanishing asymmetries (p,d = 0),
implying that X̂+ responds resonantly to the force as expected.
In the absence of coupling asymmetry, we find gr (p = 0,d) =1/(1 + d). That is, the gain in Eq. (54) is given by the damping
rate of the driven oscillator. Accordingly, for d → −1 (thedamping rate of the driven oscillator goes to zero), the gain
diverges. For the compensated system, the transfer function
modification is given bygr (p,d) =

arctan p
1(1 − d) cos
1 − d 2 + d 2 p2
2
arctan p.
(D6)− dp sin2(E1b)where  is the measurement rate introduced in Eq. (75).

gr (p,d) =(E1a)

2. Conditional variances

Also, in the fully symmetric case the system matrix M ofEq. (B9) simplifies to that given in Eq. (13), and the system ofconditional covariance equations of Eq. (81) becomes, writingout the equations for the independent elements explicitly,
V̇X+ = 2+− − γ VX+ − 4ηVX2+ + γ n̄tot ,

(E2a)

2
V̇P− = −2+− − γ VP− − 4η+−+ γ n̄tot ,

(E2b)

2
V̇X− = 2−+ − γ VX− − 4ηXX+ γ n̄tot ,2
V̇P+ = −2−+ − γ VP+ − 4η++ + γ n̄tot + ,(E2c)

˙ +− = −VX+ + VP− − γ +− − 4ηVX+ +− ,
˙ −+ = −VX− + VP+ − γ −+ − 4ηXX ++ ,


(E2d)
(E2e)

(E2f)
˙ ++ = −XX + P P − γ ++ − 4η++ VX+ , (E2g)˙ −− = −XX + P P − γ −− − 4η+− XX , (E2h)˙ XX = ++ + −− − γ XX + γ n̄d − 4ηVX+ XX ,
(E2i)
˙ P P = −++ − −− − γ P P + γ n̄d − 4η+− ++ ,(E2j)

2. Detection of mechanically nonresonant force

The second case of interest is that in which the signal force
is far-detuned (by an amount   γ ) from the mechanicalresonance frequency. The transfer function is now specified by
Eqs. (51) and (65). The modifications to the transfer functions
due to asymmetries are given by the function gn (p,d). In theuncompensated case this is
1 + p − 1 + p2
arctan p
gn (p,d) =
,
cosec
2
2 1 + p2

(D7)

while in the compensated case we findarctan p
gn (p,d) = cos
.
2

(D8)

In both cases the modulation is independent of the damping
asymmetry and so approaches 1 for vanishing coupling
asymmetry.

with n̄tot , accounting for quantum and thermal fluctuations, as
defined in Eq. (89). From Eq. (E2d), note that the observablethat is perturbed by the measurement, P̂+ , is heated at themeasurement rate , as required by Heisenberg’s uncertainty
principle. Further note that Eqs. (E2a), (E2b), and (E2e),
including the equations for the variances of most interest (VX+
and VP− ), form a closed system and analytical steady-statesolutions are readily found.
First consider the simplest case,  = 0. Writing themeasurement rate in terms of the cooperativity parameter, = γ C, we have√4ηC(n̄th + 1/2) + 1/4 − 1/2
VX+ =
4ηC
√
n̄th + 1/2 1
(E3a)
→
√ ,
√
2 η
C
(E3b)
VP+ = n̄th + 1/2 + C,

063846-21

VX− = VP− = n̄th + 1/2.

(E3c)

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)

The variance of the measured observable (X̂+ ) falls witha stronger measurement (higher cooperativity), while theconjugate observable (P̂+ ) is heated to a greater extent. The
uncoupled observables (X̂− ,P̂− ) are unaffected by the measurement. Note that if we set the thermal occupation of eachmechanical oscillator to be the same, the results of Eqs. (E3a)and (E3b) are consistent with that for a BAE measurement ofa quadrature of a single mechanical oscillator [6].

However, for both generating entangled states and force
sensing beyond quantum limits, we require that the collective
mechanical oscillation frequency, , is nonzero. Solving forthe steady-state of Eqs. (E2a), (E2b), and (E2e) in this caseleads to a fourth-order polynomial equation in VX+ . This may
be solved analytically, and expressions for +− and VP− thenfollow in turn. We find





1
1
−1 + √
8ηγ 2 C n̄tot + γ 2 − 42 + γ 2 + 42 γ 2 + 42 + 16ηγ 2 C n̄tot ,
4ηC
2γ

γ 
+− =
VX+ + 4ηCVX2+ − n̄tot ,
2
 2

2
γ
2
VP− = 2n̄tot − VX+ − 4ηCVX+ + ηCn̄tot − VX+ − 4ηCVX2+ ,
VX+ =with n̄tot given by Eq. (89). The results of Eqs. (E4a)–(E4c)also apply to the compensated asymmetric system, withthe understanding that the observables referred to are thenthe rotated observables, and we must replace parameters
without tildes with the corresponding parameters with tildes.
More useful results based on Eqs. (E4a)–(E4c) are given inSec. V E1. Note also that the second expression in Eq. (E3a)follows as a limit of Eq. (E4a) in the regime 
γ . Further,
substituting the full result for VX+ from Eq. (E4a) into theexpression for VP− in terms of VX+ from Eq. (E4c), andre-expressing the result in terms of VX+ , we find the resultof Eq. (88).
In the compensated case (with no damping asymmetry),
the solutions for the steady-state covariances are given by
Eqs. (E4a)–(E4c), with the understanding that they are now
solutions for the rotated observables, and parameters without
tildes must be replaced with the corresponding parameterswith tildes. Consequently, the asymptotic results of Eqs. (86)and (87) remain valid.
Given the solutions in Eqs. (E4a)–(E4c), Eqs. (E2g)–(E2j)form a linear system
Ai · Ri = Bij ,(E5)

where the appropriate matrices areR1 = (++ ,−− ,XX ,P P )T ,
⎡
⎤
−γ
0
− 
−γ −  ⎥⎢ 0
A1 = ⎣

 −γ
0 ⎦
− −
0
−γ
⎡
0
VX+ 0
0 +−
⎢ 0
− 4ηγ C ⎣
0
0 VX+
+− 0
0
⎤
⎡
0
⎢ 0 ⎥

,
B11 = ⎣
−γ n̄d ⎦−γ n̄d

(E6a)

⎤
0
0⎥
,
0⎦
0

(E6b)

⎤
VX+
VP ⎥
˜⎢
B12 = B11 − 2p⎣ − ⎦,
+−
+−

(E4a)
(E4b)
(E4c)

⎡

(E6d)

where B11 or B12 is the appropriate inhomogeneity in thesymmetric or compensated case, respectively. Subsequently,
Eqs. (E2c), (E2d), and (E2f) form a linear system of the formof Eq. (E5), now withR2 = (VX− ,VP+ ,−+ )T ,
⎡
⎤
−γ
0
2
−γ −2 ⎦ ,A2 = ⎣ 0
− 
−γ
⎡
⎤
2
4ηγ CXX − γ n̄tot
2
B21 = ⎣ 4ηγ C++− γ n̄tot − γ C ⎦ ,4ηγ CXX ++⎤⎡2−−˜ ⎣ 2++⎦,B22 = B21 − 2pXX + P P

(E7a)
(E7b)

(E7c)

(E7d)

again with B21 or B22 being the appropriate inhomogeneity
in the symmetric or compensated case, respectively. Notethat in the compensated case, as compared with the fully
symmetric case, there are simply additional inhomogeneousterms, leading to excess heating in the perturbed subsystem.
The solutions to these linear systems take simple forms in
the fully symmetric case, particularly in the physically relevantcase n̄d = 0 (both oscillators at the same temperature). The
linear systems defined by Eqs. (E6a)–(E6d) and (E7a)–(E7d)
are readily solved in turn, leading to++ ,−− ,XX ,P P = 0,(E8)and subsequently to(E6c)063846-22

VX− = n̄th +

22
1
+C 2
,
2
γ + 42(E9a)

TWO-MODE BACK-ACTION-EVADING MEASUREMENTS IN . . .

γ 2 + 22
1
+C 2
,
2
γ + 422γ
−+ = C.
 γ 2 + 42
VP+ = n̄th +(E9b)

matrix form
S≡

(E9c)

In the limit γ  , VX− ∼ n̄th + 1/2 and VP+ ∼ n̄th + 1/2 +C. In this limit, the measurement is of X̂+ alone, and so the
conjugate observable Pˆ+ is heated by an amount corresponding
to the cooperativity. In the opposite (and more interesting)limit,   γ , we find VX+ ∼ VP− ∼ n̄th + 1/2 + C/2. In
this case, X̂+ is directly monitored, but it is dynamicallycoupled to P̂− . Therefore, both of these collective quadratures
are effectively measured, and the corresponding conjugate
quadratures are heated by an equal amount.
APPENDIX F: UNCONDITIONAL VARIANCESAND FEEDBACK
1. Unconditional variances

The expressions for the unconditional variances, in theabsence of feedback, still take reasonably simple forms in theasymmetric case, provided we consider the regime where  
γ . In the original (uncompensated) case we find, assuming thatthe asymmetries are small,(1 + p2 )(1 + dp) − d 2 /2,
(1 + p2 )(1 − d 2 + p2 )1 + p2 − d 2 /2.
= 2ηCVX̃2+(1 + p2 )(1 − d 2 + p2 )PHYSICAL REVIEW A 87, 063846 (2013)

2X̃+= 2ηCVX̃2+

(F1a)

P̃−2

(F1b)

M + F + Q1
F

−Q1
,
M

(F2)

where (as usual) M is given by Eq. (13), Eq. (B9), or Eq. (B12),
F is the feedback matrix defined in block matrix form as(αγ /2)I2 022
,
(F3)
F≡
022
022
with In denoting the n × n identity matrix, 0mn representing
an m × n zero matrix, and 43 ],
Q1 ≡ 4η[ Q|0(F4) as defined in Eq. (80). Further, the input noisewith Q
weighting matrix in Eq. (103) may be expressed in block
matrix form as
044
Q2
√
,
(F5)
T≡
044
γ n̄tot I4
where we have introducedQ2 ≡ 4η Diag[VX+ ,+− ,XX ,++ ].

(F6)

3. Feedback: The symmetric caseIn the fully symmetric case we need only consider theobservables in the measured subsystem, such that Eq. (103)describes a four-dimensional Ornstein-Uhlenbeck process.
The state vector is thenY = (X̄+ ,P̄− ,x+ ,p− )T ,(F7)

the input noise vector isNote that there is a discrepancy between Eqs. (F1a) and (F1b)
only when there is both a coupling asymmetry and a damping
asymmetry.
2. Feedback: The general case

The equation describing both the best estimates of thecollective mechanical quadratures and the classical representations of these observables takes the form of Eq. (103). In
the general (asymmetric) case, the measured and perturbed
subspaces are coupled, such that the state vector contains
all collective quadratures and the noise vector includes
measurement noise and the noise due to the independent
mechanical environments, as specified by Eq. (104). The
system matrix in Eq. (103) may be expressed in the block[1] V. B. Braginsky, Y. I. Vorontsov, and K. P. Thorne, Science 209,
547 (1980).
[2] C. M. Caves, K. S. Thorne, R. W. P. Drever, V. D. Sandberg, andM. Zimmermann, Rev. Mod. Phys. 52, 341 (1980).
[3] V. B. Braginsky and F. Y. Khalili, Quantum Measurement(Cambridge University Press, Cambridge, UK, 1992).
[4] D. F. Walls and G. J. Milburn, Quantum Optics (Springer, Berlin,1994).
[5] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and
R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010).

 = (dW,dW,dWx+ ,dWp− )T .
dW
The system matrix is now, explicitly,⎡
γ /2 + αγ /2 + 4ηVX+ − −4ηVX+⎢
 + 4η+−γ /2 −4η+−
S≡⎢
⎣
αγ /2
0
γ /2
0
0

and the input noise matrix is now, explicitly,T ≡ Diag[ 4ηVX+ , 4η+− , γ n̄tot , γ n̄tot ].

(F8)
⎤
0
0 ⎥
⎥,
− ⎦
γ /2
(F9)

(F10)

Of course, Eqs. (F9) and (F10) could be expressed in block
matrix form, as in Eqs. (F2) and (F5).

[6] A. A. Clerk, F. Marquardt, and K. Jacobs, New J. Phys. 10,095010 (2008).
[7] G. J. Milburn and M. J. Woolley, Acta Phys. Slovaca 61, 483(2011).
[8] J. B. Hertzberg, T. Rocheleau, T. Ndukum, M. Savva, A. A.
Clerk, and K. C. Schwab, Nat. Phys. 6, 213 (2010).
[9] M. Tsang and C. M. Caves, Phys. Rev. Lett. 105, 123601(2010).
[10] M. Tsang and C. M. Caves, Phys. Rev. X 2, 031016 (2012).
[11] B. O. Koopman, Proc. Natl. Acad. Sci. USA 17, 315 (1931).

063846-23

M. J. WOOLLEY AND A. A. CLERKPHYSICAL REVIEW A 87, 063846 (2013)[12] J. Gough and M. R. James, IEEE Trans. Autom. Control 54,2530 (2009).
[13] K. Hammerer, M. Aspelmeyer, E. S. Polzik, and P. Zoller, Phys.
Rev. Lett. 102, 020501 (2009).
[14] F. Massel, S. U. Cho, J.-M. Pirkkaleinen, P. J. Hakonen, T. T.
Heikkilä, and M. A. Sillanpää, Nat. Commun. 3, 987 (2012).
[15] J. Suh, M. D. Shaw, H. G. LeDuc, A. J. Weinstein, and K. C.
Schwab, Nano Lett. 12, 6260 (2012).
[16] W. W. Johnson and M. Bocko, Phys. Rev. Lett. 47, 1184 (1981).
[17] M. F. Bocko and W. W. Johnson, Phys. Rev. Lett. 48, 1371
(1982).
[18] R. Onofrio, Phys. Lett. A 120, 1 (1987).
[19] R. Onofrio and A. Rioli, Phys. Rev. D 47, 2176 (1993).
[20] A. Rioli and U. Tambini, Phys. Lett. A 193, 348 (1994).
[21] C. Cinquegrana, E. Majorana, N. Pergola, P. Puppo,
P. Rapagnani, and F. Ricci, Phys. Rev. D 50, 3596 (1994).
[22] P. Bonifazi, C. Cinquegrana, E. Majorana, N. Pergola, P. Puppo,
P. Rapagnani, F. Ricci, S. Vaselli, and M. Visco, Phys. Lett. A
215, 141 (1996).
[23] J.-W. Pan, D. Bouwmeester, H. Weinfurter, and A. Zeilinger,
Phys. Rev. Lett. 80, 3891 (1998).
[24] C. W. Chou, H. de Riematten, D. Felinto, S. V. Polyakov,
S. J. van Enk, and H. J. Kimble, Nature (London) 438, 828(2005).
[25] A. Kuzmich, L. Mandel, and N. P. Bigelow, Phys. Rev. Lett. 85,
1594 (2000).
[26] K. Lalumiére, J. M. Gambetta, and A. Blais, Phys. Rev. A 81,
040301(R) (2010).
[27] A. Frisk Kockum, L. Tornberg, and G. Johansson, Phys. Rev. A
85, 052318 (2012).
[28] L.-M. Duan, G. Giedke, J. I. Cirac, and P. Zoller, Phys. Rev.
Lett. 84, 2722 (2000).
[29] J. D. Teufel, T. Donner, D. Li, J. W. Harlow, M. S. Allman,
K. Cicak, A. J. Sirois, J. D. Whittaker, K. W. Lehnert, and R. W.
Simmonds, Nature (London) 475, 359 (2011).
[30] H. J. Kimble, Phys. Scr. T 76, 127 (1998).
[31] A. Einstein, B. Podolsky, and N. Rosen, Phys. Rev. 47, 777(1935).
[32] C. Weedbrook, S. Pirandola, R. Garcia-Patrón, N. J. Cerf, T. C.
Ralph, J. H. Shapiro, and S. Lloyd, Rev. Mod. Phys. 84, 621(2012).
[33] A. D. O’Connell, M. Hofheinz, M. Ansmann, R. C. Bialczak,
M. Lenander, E. Lucero, M. Neeley, D. Sank, H. Wang,
M. Weides, J. Wenner, J. M. Martinis, and A. N. Cleland, Nature(London) 464, 697 (2010).
[34] K. C. Lee, M. R. Sprague, B. J. Sussman, J. Nunn, N. K.
Langford, X.-M. Jin, T. Champion, P. Michelberger, K. F. Reim,
D. England, D. Jaksch, and I. A. Walmsley, Science 334, 1253
(2011).
[35] B. Julsgaard, A. Kozhekin, and E. S. Polzik, Nature (London)
413, 400 (2001).

[36] J. D. Jost, J. P. Home, J. N. Amini, D. Hanneke, R. Ozeri,
C. Langer, J. J. Bollinger, D. Leibfried, and D. J. Wineland,
Nature (London) 459, 683 (2009).
[37] D. Vitali, S. Gigan, A. Ferreira, H. R. Böhm, P. Tombesi,
A. Guerreiro, V. Vedral, A. Zeilinger, and M. Aspelmeyer, Phys.
Rev. Lett. 98, 030405 (2007).
[38] S. Mancini, V. Giovannetti, D. Vitali, and P. Tombesi, Phys. Rev.
Lett. 88, 120401 (2002).
[39] M. Pinard, A. Dantan, D. Vitali, O. Arcizet, T. Briant, and
A. Heidmann, Europhys. Lett. 72, 747 (2005).
[40] H. Müller-Ebhardt, H. Rehbein, R. Schnabel, K. Danzmann, and
Y. Chen, Phys. Rev. Lett. 100, 013601 (2008).
[41] D. Vitali, S. Mancini, L. Ribichini, and P. Tombesi, J. Opt. Soc.
Am. B 20, 1054 (2003).
[42] M. J. Hartmann and M. B. Plenio, Phys. Rev. Lett. 101, 200503(2008).
[43] S. Pirandola, D. Vitali, P. Tombesi, and S. Lloyd, Phys. Rev.
Lett. 97, 150403 (2006).
[44] W. Niedenzu, R. M. Sandner, C. Genes, and H. Ritsch, J. Phys.
B 45, 245501 (2012).
[45] K. Borkje, A. Nunnenkamp, and S. M. Girvin, Phys. Rev. Lett.
107, 123601 (2011).
[46] J. Zhang, K. Peng, and S. L. Braunstein, Phys. Rev. A 68, 013808(2003).
[47] C. Joshi, J. Larson, M. Jonson, E. Andersson, and P. Öhberg,
Phys. Rev. A 85, 033805 (2012).
[48] M. Schmidt, M. Ludwig, and F. Marquardt, New J. Phys. 14,125005 (2012).
[49] H. M. Wiseman and G. J. Milburn, Quantum Measurement andControl (Cambridge University Press, Cambridge, UK, 2009).
[50] H. M. Wiseman and G. J. Milburn, Phys. Rev. A 47, 642 (1993).
[51] S. Mancini and H. M. Wiseman, Phys. Rev. A 75, 012330 (2007).
[52] S. L. Braunstein and P. van Loock, Rev. Mod. Phys. 77, 513(2005).
[53] G. Vidal and R. F. Werner, Phys. Rev. A 65, 032314 (2002).
[54] V. Giovannetti, S. Mancini, D. Vitali, and P. Tombesi, Phys. Rev.
A 67, 022320 (2003).
[55] A. C. Doherty and K. Jacobs, Phys. Rev. A 60, 2700 (1999).
[56] T. A. Palomaki, J. W. Harlow, J. D. Teufel, R. W. Simmonds,
and K. W. Lehnert, Nature (London) 495, 210 (2013).
[57] C. Lang, C. Eichler, L. Steffen, J. M. Fink, M. J. Woolley,
A. Blais, and A. Wallraff, Nat. Phys. 9, 345 (2013).
[58] R. Vijay, C. Macklin, D. H. Slichter, S. J. Weber, K. W. Murch,
R. Naik, A. N. Korotkov, and I. Siddiqi, Nature (London) 490,
77 (2012).
[59] P. Campagne-Ibarcq, E. Flurin, N. Roch, D. Darson, P. Morfin,
M. Mirrahimi, M. H. Devoret, F. Mallet, and B. Huard, Phys.
Rev. X 3, 021008 (2013).
[60] K. Zhang, P. Meystre, and W. Zhang, arXiv:1304.2459.
[61] M. J. Woolley, A. C. Doherty, G. J. Milburn, and K. C. Schwab,
Phys. Rev. A 78, 062303 (2008).

063846-24