Optomechanical dual-beam backaction-evading measurement beyond the rotating-wave approximation
OOptomechanical dual-beam backaction-evading measurementbeyond the rotating-wave approximation
Daniel Malz and Andreas Nunnenkamp
Cavendish Laboratory, University of Cambridge, Cambridge CB3 0HE, United Kingdom (Dated: November 7, 2018)We present the exact analytical solution of the explicitly time-periodic quantum Langevin equation describingthe dual-beam backaction-evasion measurement of a single mechanical oscillator quadrature due to Braginsky,Vorontsov and Thorne beyond the commonly used rotating-wave approximation. We show that counterrotatingterms lead to extra sidebands in the optical and mechanical spectra and to a modification of the main peak.Physically, the backaction of the measurement is due to periodic coupling of the mechanical resonator to a lightfield quadrature that only contains cavity-filtered shot noise. Since this fact is independent of other degreesof freedom the resonator might be coupled to, our solution can be generalized, including to dissipatively orparametrically squeezed oscillators, as well as recent two-mode backaction evasion measurements.
I. INTRODUCTION
A continuous measurement of the position of a harmonicoscillator is subject to the “standard quantum limit” (SQL),a limit directly imposed by Heisenberg’s uncertainty rela-tion [1, 2]. An observable that can be monitored without pre-cision limit is called “quantum non-demolition” (QND) vari-able, such that its continuous measurement can avoid the mea-surement backaction (BA) [3] and thus open the way to thedetection of weak forces, such as those due to gravitationalwaves [4].There has been continued interest to implement backaction-evading (BAE) measurements [5]. Following a detailed theo-retical proposal, the first demonstration with a sensitivity be-yond the SQL was in optomechanics [6, 7] and they have sinceproven very useful [8–11]. Despite the importance of suchmeasurements, Ref. [6] discusses only lowest-order correc-tions to the rotating-wave approximation (RWA).Here, using the recently developed Floquet approach [12],we derive the exact solution to the equations describing a BAEmeasurement. Due to the presence of CR terms, this consti-tutes a solution to genuinely explicitly time-dependent quan-tum Langevin equations, as there is no frame in which theybecome stationary. The solution is possible because the me-chanical oscillator couples solely to a light quadrature that isindependent of the mechanical oscillator, and only containsfiltered shot noise. The coupling is periodic, which leads totwo-time correlators that are not time-translation invariant. Insuch a situation, the power spectrum is a time average of theFourier transform of the autocorrelator [12].In the following, after introducing our model in Sec. II, webriefly remind the reader of results in RWA in Sec. III. Then,we present all aspects of the solution, from the Floquet frame-work to the resulting spectra and backaction in Sec. IV. Fi-nally, in Sec. V, we show how to generalize our solution tothe important cases of dissipative and parametric squeezing,as well as two-mode BAE measurements (with explicit calcu-lations in Appendices A to D). We conclude in Sec. VI. ω − ω + ω cav − ω m − ω m ω m δ = ω m FIG. 1. The dual-beam backaction-evading (BAE) measurementscheme proposed in Ref. [3, 6]. Two lasers of equal strength drivethe cavity at frequencies ω ± = ω cav ± δ , where δ = ω m such thattheir sidebands overlap at the cavity frequency ω cav . Counterrotating(CR) terms cause peaks at ω cav ± ω m to appear. The output lightwill be an approximate BAE measurement of a rotating quadratureof the mechanical oscillator, the CR terms being responsible for thefinite backaction. II. MODEL
In cavity optomechanics, the photons in a cavity couple tothe motion of a mirror via radiation pressure. A BAE mea-surement of one of the quadratures of the oscillator can beimplemented by applying two drives of equal strength at fre-quencies ω cav ± ω m (see Fig. 1). Here, ω cav ( ω m ) is the fre-quency of a cavity (mechanical) mode. Originally due to Bra-ginsky et al. [3], the first detailed analysis in the context ofcavity optomechanics was given in Ref. [6].We consider the Hamiltonian of a standard cavity optome-chanical system H = H OM + H rest + H drives + H baths , (1)where ( (cid:126) = 1 ) H OM = ω cav a † a + ω m b † b − g a † a ( b † + b ) , (2a) H drives = α ( e − iω + t + e − iω − t ) a † + h.c. (2b) a, b are the bosonic annihilation operators of the cavity modeand the mechanical oscillator, respectively. The cavity modefrequency is ω cav , the mechanical frequency ω m , the couplingstrength via radiation pressure g , and the driving strength ofthe drives with frequencies ω ± = ω cav ± δ is α , where we haveleft δ unspecified for now (Fig. 1). The BAE measurement isrealized for δ = ω m . A detailed derivation of the optomechan-ical Hamiltonian can be found for instance in Ref. [13]. H baths couples the oscillator and cavity modes to baths at finite and a r X i v : . [ c ond - m a t . m e s - h a ll ] N ov zero temperature, respectively. H rest will remain unspecifiedfor now, it can contain terms that couple the harmonic oscilla-tor to other degrees of freedom. The only requirement is that [ H rest , a + a † ] = 0 , i.e., that the other degrees of freedom donot couple to this measurement cavity quadrature.To proceed, we split the light field into a coherent part andfluctuations, go into a rotating frame, a = ¯ ae − iω cav t ( e iδt + e − iδt + d ) , and linearize the Hamiltonian. Without loss ofgenerality, δ > . Under the usual assumptions of Markovianbaths and a one-sided cavity, the resulting Hamiltonian H = ω m b † b − G cos( δt )( b + b † )( d + d † ) (3)gives rise to Langevin equations [2, 14] that are periodic intime ˙ d = − κ d + √ κd in + 2 iG cos( δt )( b + b † ) , (4a) ˙ b = (cid:16) − iω m − γ (cid:17) b + √ γb in + 2 iG cos( δt )( d + d † ) + i [ H rest , b ] . (4b)Here, we have defined the enhanced optomechanical cou-pling constant G = g ¯ a , and the mechanical (optical)dissipation rate γ ( κ ). b in , d in are input noise operatorswith (cid:104) d in ( t ) d † in ( t (cid:48) ) (cid:105) = δ ( t − t (cid:48) ) , (cid:104) d † in ( t ) d in ( t (cid:48) ) (cid:105) = 0 , (cid:104) b in ( t ) b † in ( t (cid:48) ) (cid:105) = ( n th + 1) δ ( t − t (cid:48) ) , and (cid:104) b † in ( t ) b in ( t (cid:48) ) (cid:105) = n th δ ( t − t (cid:48) ) .In Sec. IV A, we solve Eqs. (4a) and (4b) without furtherapproximations, in particular without the rotating wave ap-proximation (RWA). III. BACKACTION EVASION IN RWA
As a reminder, we first consider a BAE measurement withinRWA. We define a quadrature rotating a frequency δ as X δ = be iδt + b † e − iδt . (5)The cavity equation of motion is (cf. Eq. (4a)) ˙ d = − κ d + √ κd in + iG (cid:0) be iδt + b † e − iδt (cid:1) . (6)The cavity couples only to one quadrature of the mechanicaloscillator X δ (set by the phase relation of the external drives).The equation of motion for that quadrature can be obtainedfrom Eq. (4b) (here also in RWA) ˙ X δ = − γ X δ + √ γX δ, in + i [ H rest , X δ ]+ iG (cid:104) d (cid:16) e i ( ω m − δ ) t − e i ( δ − ω m ) t (cid:17) − h.c. (cid:105) . (7)If δ = ω m , the term in square brackets in the second linevanishes, and X ω m is entirely unaffected by the cavity. Wecan readily solve the equation of motion for d , Eq. (6) d ( ω ) = χ c ( ω ) (cid:2) √ κd in ( ω ) + iGX δ ( ω ) (cid:3) , (8) with cavity response function χ c ( ω ) = ( κ/ − iω ) − . Thusthe optical spectrum is directly related to the quadrature spec-trum S d † d ( ω ) = | χ c ( ω ) | G S X δ X δ ( ω ) , (9)where for a stationary process the Wiener-Khinchin theoremensures that the Fourier transform S AB ( ω ) = (cid:90) ∞−∞ d t e iωt (cid:104) A ( t ) B (0) (cid:105) (10)is the quantum noise spectral density. The scheme corre-sponds to a true BAE measurement, if the mechanical quadra-ture rotating with frequency δ is a QND variable. Here, this isthe case for δ = ω m (or, if applicable, δ = ω m, eff , the effec-tive mechanical frequency). Independent of RWA, the input-output relation d out = d in − √ κd can be used to the obtain theoutput spectrum S d † out d out ( ω ) = κS d † d ( ω ) . For details on howto measure and interpret the output spectrum, see Ref. [15].
IV. FLOQUET FRAMEWORK
The analysis in Sec. III above breaks down if we includecounterrotating (CR) terms. Using the framework developedin Ref. [12], it is possible to find the stationary state of theexplicitly time-periodic quantum Langevin equations exactly,which gives us the opportunity to precisely determine thebackaction of the BAE measurement.First, we split system operators into Fourier components d ( t ) = ∞ (cid:88) n = −∞ e inδt d ( n ) ( t ) , (11a) d † ( t ) = ∞ (cid:88) n = −∞ e inδt d ( n ) † ( t ) , (11b)and define d ( n ) ( ω ) = (cid:90) ∞−∞ d t e iωt d ( n ) ( t ) , (12a) d ( n ) † ( ω ) = (cid:90) ∞−∞ d t e iωt d ( n ) † ( t ) . (12b)Note that this convention results in [ d ( n ) ( ω )] † = d ( − n ) † ( − ω ) .The reason for doing this is that Eqs. (4a) and (4b) can nowbe written as an infinite set of coupled, stationary Langevinequations ˙ d ( n ) = (cid:16) − inδ − κ (cid:17) d ( n ) + δ n, √ κd in + iG (cid:16) b ( n − + b ( n +1) † + λb ( n +1) + λb ( n − † (cid:17) , (13a) ˙ b ( n ) = (cid:16) − iω m − inδ − γ (cid:17) b ( n ) + δ n, √ γb in + iG (cid:16) I ( n +1) + λI ( n − (cid:17) + i [ H rest , b ] ( n ) . (13b)where x ( n ) ≡ b ( n ) + b ( n ) † , I ( n ) ≡ d ( n ) + d ( n ) † , δ n, is theKronecker delta, and [ H rest , b ] ( n ) is the n th Fourier componentof the commutator. Where feasible, here and in the following,we will mark counterrotating (i.e., off-resonance) terms by a λ (note that this only makes sense when δ is close to ω m ), suchthat RWA corresponds to λ = 0 and the full solution to λ = 1 .To further guide the intuition, we remark that the second lineof Eq. (13a) equals iG ( x ( n − + x ( n +1) ) . One can think ofthe two laser drives to result in two separate couplings to theposition x of the resonator.We define the “spectrum” S A † A ( ω, t ) ≡ (cid:90) ∞−∞ d τ e iωτ C AA ( τ, t ) , (14)with C AA ( τ, t ) ≡ (cid:10) A † ( t + τ ) A ( t ) (cid:11) . The time dependence of S A † A ( ω, t ) can be expressed as a Fourier series [12] S A † A ( ω, t ) = ∞ (cid:88) n = −∞ e iδnt S ( n ) A † A ( ω ) (15)with components S ( m ) A † A ( ω ) = (cid:88) n (cid:90) d ω (cid:48) π (cid:68) A ( n ) † ( ω + nδ ) A ( m − n ) ( ω (cid:48) ) (cid:69) . (16)It can be shown that the zeroth Fourier component, which atthe same time is the time average of S ( ω, t ) , is the measuredpower spectrum [12]. A. Solution beyond RWA
To solve Eqs. (4a) and (4b), we define the optical quadra-tures I = d + d † and Q = − i ( d − d † ) . From Eq. (4a) weobtain ˙ I = − κ I + √ κI in , (17a) ˙ Q = − κ Q + √ κQ in + 4 G cos( δt ) x. (17b)For an illustration how the mechanical and optical quadra-tures couple together, see Fig. 2. For example, the arrowfrom X δ to Q indicates that X δ appears in the equation ofmotion of Q and that the latter therefore is influenced by theformer. Since I commutes with the system Hamiltonian (3), [ H, I ] = 0 , there are no arrows pointing toward it in Fig. 2and we can solve its equation of motion (17a) directly I ( ω ) = √ κχ c ( ω ) I in ( ω ) , and I ( n ) ( ω ) = δ n, I ( ω ) . (18)Thus (cid:104) I ( ω ) I ( ω (cid:48) ) (cid:105) = 2 πδ ( ω + ω (cid:48) ) | χ c ( ω ) | and (cid:104) I ( ω ) (cid:105) = 0 .Equations (17a) and (18) imply that I is shot noise, filtered bythe cavity, and is independent of the mechanics. Furthermore,this result does not depend on H rest , as long as [ H rest , I ] = 0 .The coupling to I is periodic (cf. Eq. (4b)), a consequence ofamplitude beating of the coherent state in the cavity. There-fore, the measurement has the same effect on the mechanicalresonator as a time-periodic coupling to filtered shot noise. lightmechanics ˆ X δ ˆ P δ ˆ I ˆ Q conjugateconjugate grey: periodicinfo BA CR terms FIG. 2. Coupling of quadratures. ˆ Q records the desired information(info) about ˆ X δ . The measurement backaction (BA) is acting on ˆ P δ .Marked in gray are coupling due to counterrotating (CR) terms thatare neglected in RWA. With Eq. (18) we can solve Eq. (13b) in the case H rest = 0 by going into frequency space b ( n ) ( ω ) = χ m ( ω − nδ ) (cid:104) δ n, √ γb in + 2 √ κGf ( n ) in (cid:105) , (19a) b ( n ) † ( ω ) = χ ∗ m ( − ω + nδ ) (cid:104) δ n, √ γb † in − √ κGf ( n ) in (cid:105) , (19b)where the new bath noise operators f ( n ) in ( ω ) = iχ c ( ω )( δ n, + λδ n, − ) I in ( ω ) / . (20)They obey f ( n ) † in = − f ( n ) in , f ( − n ) in = f ( n ) in , and (cid:104) f in (cid:105) = 0 . Thetime-dependence is best seen in the time domain, where f in ( t ) = i cos( δt ) I ( t ) / √ κ. (21)This expression explicitly contains the time-dependent cou-pling and the filtered shot noise I ( t ) . The correlator is (cid:104) f in ( t ) f in ( t (cid:48) ) (cid:105) = − e − κ | t − t (cid:48) | / cos( δt ) cos( δt (cid:48) ) . (22)For stationary noise, the RHS of Eq. (22) would have to de-pend solely on the difference t − t (cid:48) .We can rewrite the equation of motion for b in terms of thenew input f in ˙ b = (cid:16) − iω m − γ (cid:17) b + √ γb in + 2 √ κGf in + i [ H rest , b ] . (23)From this equation it is clear that we are always able to passfrom the “unmeasured system” ( G = 0 ) to the “measured”one ( G (cid:54) = 0 ) by substituting √ γb in → √ γb in + 2 √ κGf in . (24)We use this to generalize the solution in Sec. V. B. Importance of Floquet framework
At this point it is useful to reflect on the advantage of theFloquet approach. It might appear as if it were unnecessary,since with Eqs. (17a) and (18), we can already write down ˙ b = (cid:16) − iω m − γ (cid:17) b + √ γb in + 2 iG cos( δt ) I ( t ) + i [ H rest , b ] , (25)which for H rest = 0 has the solution b ( ω ) = χ m ( ω ) {√ γb in ( ω ) + iG [ I ( ω − δ ) + I ( ω + δ )] } . (26)Since x ( ω ) = b ( ω ) + b † ( ω ) , this might prompt us to derivethe “power spectrum” S xx ( ω ) = (cid:90) d ω (cid:48) π (cid:104) x ( ω ) x ( ω (cid:48) ) (cid:105) . (27)This, however, is not the correct power spectrum, since x ( t ) does not describe a stationary stochastic process. It is pos-sible to remove the non-stationary terms manually (abovethey will be (cid:104) I ( ω − δ ) I ( ω (cid:48) − δ ) (cid:105) and (cid:104) I ( ω + δ ) I ( ω (cid:48) + δ ) (cid:105) )and thus obtain the stationary part, which is the power spec-trum [12]. In contrast, the systematic solution via Fouriercomponents, which all obey stationary Langevin equations,is well-defined. In addition, the Fourier components are moreversatile, and allow writing down the spectrum in arbitraryrotating frames [12]. They also provide more intuition, sincedifferent Fourier components tend to have different physicalorigins. Finally, if no exact solution is viable, they simplifythe process of approximating to desired order. C. Mechanical quadrature spectrum
We define a rotating quadrature as before (see Eq. (5)). Itsspectrum is [12] S (0) X δ X δ ( ω ) = S ( − bb ( ω + δ ) + S (0) bb † ( ω + δ )+ S (2) b † b † ( ω − δ ) + S (0) b † b ( ω − δ ) , (28)where S AB is defined through Eq. (16). We have only writ-ten down the zeroth Fourier component (the stationary part),because that is the part of the spectrum that will be measured.We can split it up into two parts S (0) X δ X δ ( ω ) = S RWA X δ X δ ( ω ) + S CR X δ X δ ( ω ) , (29)with S RWA X δ X δ being the result in RWA (dependent on H rest ), and S CR X δ X δ being due to CR terms. Note that S RWA X δ X δ is unchangedfrom the unmeasured case, since this is the essence of a BAEmeasurement. For H rest = 0 , we obtain our first major result S CR X δ X δ ( ω ) = κλ G (cid:0) | χ m ( ω + δ ) χ c ( ω + 2 δ ) | + | χ m ( − ω + δ ) χ c ( − ω + 2 δ ) | (cid:1) . (30)If H rest is nonzero, but couples weakly, such that its effect iswell approximated by introducing an effective damping con-stant and mechanical frequency, the analysis is unchanged,such that the correction takes the same form S CR X δ X δ ( ω ) = κλ G (cid:8) | χ m, eff ( ω + δ ) χ c ( ω + 2 δ ) | + | χ m, eff ( − ω + δ ) χ c ( − ω + 2 δ ) | (cid:9) , (31) − − − × frequency ω/κ FIG. 3. Exact optical output spectrum κS (0) d † d Eq. (32) (blue, solid) incomparison to a perfect BAE measurement of the modified mechan-ical spectrum (31) from Sec. IV C (yellow, dashed), and the RWAresult Eq. (9) (red, dotted). The full spectrum (blue, solid) is calcu-lated via Eq. (32), where additional terms due to the counterrotatingquadrature X − δ appear. The sidepeaks of the modified mechanicalspectrum do not show in the yellow curve, as they are suppressed by | χ c ( ω ) | . Parameters Γ eff /κ = 1 , ω m /κ = 2 , n eff = 10 , G/κ =10 , δ = ω m (thus C = 400 ). Note that the large value of Γ eff ischosen for visibility. but with an effective susceptibility χ m, eff ( ω ) = [Γ eff − i ( ω − ω m, eff )] − . Note that if δ (cid:54) = ω m, eff , the measured quadra-ture rotates at a different frequency than the natural oscillatorquadratures and thus the measurement backaction (BA) willcontaminate the measurement at later times, even in RWA.Only the case δ = ω m, eff is backaction-evading (BAE). Sincethis is our main interest, we will fix δ = ω m, eff in the follow-ing. For the general case δ (cid:54) = ω m, eff , see Eq. (F1). The changein the spectrum Eq. (31) includes a modification of the mainpeak, and new peaks corresponding to the upper sideband ofthe blue drive and the lower sideband of the red drive. D. Optical spectrum
As we have seen, CR terms modify the mechanical quadra-ture spectrum, which is reflected in the optical spectrum S (0) d † d ( ω ) = G | χ c ( ω ) | (cid:104) S (0) xx ( ω − δ ) + S (0) xx ( ω + δ )+ S (2) xx ( ω − δ ) + S ( − xx ( ω + δ ) (cid:105) = G | χ c ( ω ) | (cid:104) S (0) X δ X δ ( ω ) + λ S (0) X − δ X − δ ( ω )+ λ (cid:16) S (0) X δ X − δ ( ω ) + S (0) X − δ X δ ( ω ) (cid:17)(cid:105) . (32)Comparing the second line of Eq. (32) to Eq. (9), we noticethat there are extra terms present, namely those with at leastone λ out front. They are due to additional terms containing X δ and P δ on the RHS of the equation of motion for Q . InFig. 2, they are denoted by the two gray arrows pointing to Q .Thus, this is an imperfection of the BAE measurement.In fact, the additional contributions stem from the spectrumof the counterrotating quadrature X − δ (i.e., with frequency − δ ) and correlations of that quadrature with X δ . That promptsa remarkably literal interpretation of “counterrotating terms”.It also clarifies the origin of the measurement BA. The CRquadrature can be written X − δ ( t ) = cos(2 δt ) X δ ( t ) + sin(2 δt ) P δ ( t ) , (33)which makes it apparent that the measurement picks up someinformation about P δ , the quadrature conjugate to X δ .In contrast, the measurement backaction , as discussed inSec. IV C, arises from the gray arrows emanating from I . Thiscorrection is contained within S X δ X δ in Eq. (32).In order to distinguish the two contributions (the imperfec-tion and BA), we plot three functions in Fig. 3. The full opticalspectrum (32) in blue (solid) encompasses both contributions.A perfect BAE measurement of the modified mechanical spec-trum Eqs. (9) and (31) is shown in yellow (dashed), it picks upthe thermal contribution and the BA. Finally, the expected re-sult in RWA is red (dotted), which corresponds to a perfectBAE measurement of an otherwise undriven mechanical os-cillator. The sidepeaks of the modified mechanical spectrumdo not show in the yellow curve, as they are suppressed by | χ c ( ω ) | . E. Mechanical and optical variances
An important application of BAE measurements is deter-mining the quadrature variance, necessary for the verificationof quantum squeezing. The BA imperils this by increasing thevariance to (cid:10) X δ (cid:11) = (cid:10) X RWA (cid:11) + 2 n CR , (34)where (cid:10) X RWA (cid:11) = (cid:10) X unmeasured (cid:11) is the result calculated with-out CR terms and (cf. (31)) n CR ≡ (cid:90) d ω π S CR X δ X δ ( ω ) = 8 G ( κ/ Γ eff + 1)(Γ eff + κ ) + 16 ω m, eff . (35)Perhaps the most surprising aspect about this result is that n CR is independent of the temperature of the mechanical bath,whether it is squeezed or not, and what quadrature we mea-sure. This fact is already realized on the level of the spectrumcorrection Eq. (31). The reason, as we have noted above, isthat the optical quadrature I in the measurement cavity is in-dependent of the mechanics and that the BA is solely due tothis quadrature (see Fig. 2). Note that although we have writ-ten the change in terms of a number of phonons n CR , it is notthermal heating that causes this effect, but rather the extrac-tion of information about the conjugate quadrature.In the optical spectrum, we saw that sidepeaks appear (cf.Fig. 3). To get an approximation to the variance of the mea-sured mechanical quadrature, we integrate over the main peakof Eq. (32). Here, we calculate the error of this method, whichcan be used for underdamped oscillators Γ eff (cid:28) ω m, eff . The weight of the main peak is (cid:10) X meas (cid:11) = (cid:10) X RWA (cid:11) + 2 n CR ω m, eff (2Γ eff + κ )(Γ eff + κ )(Γ eff + 16 ω m, eff ) , (36)with n CR given by Eq. (35). It can be experimentally accessedby integrating S (0) d † out d out ( ω ) / ( κG | χ c ( ω ) | ) from − ω m, eff to ω m, eff , see Appendix E. The second term is the BA. For ex-ample, if Γ eff /κ = 10 − , G/κ = 0 . , and ω m /κ = 5 (i.e., C = 10 ), the variance increases by . (in units of thezero-point fluctuations of the mechanical oscillator). Since n CR ∝ G , stronger driving quickly makes this effect notice-able. In typical experimental regimes, where Γ eff (cid:28) ( ω m , κ ) ,the BA in Eq. (36) is well approximated by n CR , such that (cid:10) X meas (cid:11) ≈ (cid:10) X δ (cid:11) (34), as desired.Finally, we would like to gain some insights about the goodand bad cavity limits of Eqs. (35) and (36). To give some intu-ition what these limits imply for the optical spectrum, considerthe solution of Eq. (4a) d ( t ) = (cid:90) t −∞ d τ e κ ( τ − t )2 [2 iG cos( ω m τ ) x ( τ ) + d in ( τ )] . (37)In the good cavity limit ( ω m (cid:29) κ ), the exponential decay isnegligible over a period π/ω m . Thus (cid:90) πωm d τ cos( ω m τ )[cos( ω m τ ) X δ ( τ ) + sin( ω m τ ) P δ ( τ )] ≈ (1 / (cid:90) πωm d τ X δ ( τ ) , (38)since X δ ( ω ) and P δ ( ω ) are centered around ω = 0 . On theother hand, in the bad cavity limit κ > ω m the photons leavethe cavity faster than the change in coupling parameter, suchthat no averaging takes place. A simple argument then showsthat the contribution from X δ is roughly three times as big asthe contribution from P δ .These properties are reflected in the variances Eqs. (35)and (36). In the bad cavity limit, where κ (cid:29) ω m, eff , the num-ber of added phonons is n CR → G Γ eff (Γ eff + κ ) . (39)In this regime the separation of the two drives is small com-pared to the bandwidth of the cavity, and therefore the BA be-comes significant ( n CR ∼ ) at a cooperativity C ∼ , where C ≡ G /κ Γ eff . As we have seen, the reason is that resonantand CR terms couple with equal strength.In contrast, the good cavity limit is κ (cid:29) Γ eff and ω m, eff (cid:29) κ , where Eq. (35) reduces to n CR → κG eff ω m, eff . (40)This agrees with the perturbative result in Ref. [6]. Equa-tion (40) tells us that the BA depends inversely on the effec-tive mechanical dissipation rate. Physically, this is because Γ eff , the rate at which the mechanical oscillator relaxes to itssteady state, competes with the BA rate due to CR terms,which is independent of Γ eff . The BA also decreases withincreasing ω m , as CR terms become less resonant. We findthat the BA becomes significant ( n CR ∼ ) for a cooperativity C ∼ ω m /κ . Therefore, in the good cavity regime, the CRterms are suppressed by a factor ∼ ω m /κ in comparison toresonant terms, due to the mechanism described above. V. GENERALIZATION
Above, we have taken H rest = 0 for simplicity, but we canin fact take any H rest and write down the solution, as long aswe know the steady state of H OM + H rest + H baths , i.e., withoutthe BAE measurement ( G = 0 ). Furthermore, we need [ H rest , I ] = 0 . (41)If the unmeasured solution is also periodic, the considerationin [12] applies: For incommensurate periods, only the station-ary part will be picked up by the measurement, but if they arecommensurate, some other Fourier components may enter theoutput spectrum.To formulate a general theory, we collect all system (input)operators into a vector (cid:126)F ( (cid:126)F in ). The Langevin equations read ˙ (cid:126)F ( t ) = A ( t ) (cid:126)F ( t ) + L (cid:126)F in ( t ) . (42)For a time-independent A ( t ) = A , a Fourier transform yields (cid:126)F ( ω ) = χ ( ω ) L (cid:126)F in ( ω ) , (43)with susceptibility matrix χ ( ω ) = ( − iω − A ) − . We cando the same for a time-periodic A ( t ) = A ( t + T ) once wehave reformulated the problem in terms of Fourier compo-nents. Then (cid:126)F ( n ) ( ω ) = ∞ (cid:88) m = −∞ χ ( n − m ) ( ω − nδ ) L (cid:126)F ( m ) in ( ω ) . (44)If we know how to solve the unmeasured system ( G = 0 ),we know how to find χ ( n ) in that case. Importantly, Eq. (44)does not make any assumption about the noise (cid:126)F ( n ) in , apartfrom stationarity. For this reason, the replacement in Eq. (24)above will leave the susceptibility matrix χ ( n ) ( ω ) unchangedand therefore can be used to calculate the measurement BA forall systems with suitable H rest . In the end, we can write downthe scattering matrix in terms of the susceptibility matrix S ( n ) ( ω ) = δ ,n − Lχ ( n ) ( ω ) L , (45)which can be use to calculate the output fields (cid:126)F ( n ) out ( ω ) = (cid:88) m S ( n − m ) ( ω − nδ ) (cid:126)F ( m ) in ( ω ) . (46)For examples how χ can look like, see the detailed calcu-lations in Appendices A to D. Appendix A contains the setup used to produce bichromatic squeezing [9, 10, 12, 16–19], Ap-pendix B discusses a weak-coupling version of the former,which is easily adapted to other weakly coupled systems, andAppendix C contains the parametric squeezing case, which isan example with more than one independent, relevant com-ponent in χ . Last is Appendix D, which extends the methodhere to two mechanical modes, and covers a recently devel-oped two-mode BAE measurement [20, 21]. VI. CONCLUSION
Using the framework from Ref. [12], we derived the fullsolution of an optomechanical system subject to a dual-beamBAE measurement [5, 7–10]. This enables us to calculate themodification of the spectrum and quantify the measurementbackaction precisely. Furthermore, we demonstrate that ourtechnique is versatile, by showing how to generalize the cal-culation to systems where the mechanical resonator is addi-tionally coupled to other degrees of freedom and illustrate thetechnique with several examples.
ACKNOWLEDGMENTS
We are grateful to Chan Lei and John Teufel for stimulatingand insightful discussions. A.N. holds a University ResearchFellowship from the Royal Society and acknowledges addi-tional support from the Winton Programme for the Physicsof Sustainability. D.M. acknowledges support by the UK En-gineering and Physical Sciences Research Council (EPSRC)under Grant No. EP/M506485/1.
Appendix A: Dissipative bichromatic squeezing
An important application of the type of BAE measurementdiscussed in this article is the verification of quantum squeez-ing in mechanical resonators, e.g., Refs. [9, 10]. Here wegeneralize our method to this squeezing scheme, proposed in[16, 17], which has successfully produced quantum squeezingof the resonator [9, 10, 18, 19]. In this case H rest = − ∆ c † c − (cid:2) c (cid:0) G + e iδt b + G − b † (cid:1) + h.c. (cid:3) , (A1)where we have already displaced and linearized. c is the an-nihilation operator of another cavity in a frame rotating withthe lower-frequency drive. Furthermore, we have applied theRWA to H rest . The governing Langevin equations are ˙ c = (cid:16) i ∆ − κ (cid:17) c + √ κ c in + i (cid:0) G − b + G + e − iδt b † (cid:1) , (A2a) ˙ b = (cid:16) − iω m − γ (cid:17) b + √ γb in + i (cid:0) G − c + G + e − iδt c † (cid:1) , (A2b)where c in corresponds to a vacuum (zero temperature) bath,and b in to a finite temperature bath with mean occupation n th ,so (cid:10) c in ( t ) c in ( t (cid:48) ) † (cid:11) = δ ( t − t (cid:48) ) , (cid:10) b in ( t ) † b in ( t (cid:48) ) (cid:11) = δ ( t − t (cid:48) ) n th , (cid:10) b in ( t ) b in ( t (cid:48) ) † (cid:11) = δ ( t − t (cid:48) )(1+ n th ) (other correlators are zero). In the Floquet ansatz, Eq. (A2b) can be expressed in termsof a block-diagonal infinite matrix [12] with blocks χ − opt ( ω − nδ ) − iG − − iG + − iG − χ − m ( ω − nδ ) − iG + iG + χ − ∗ opt ( − ω + ( n + 2) δ ) iG − iG + iG − χ − ∗ m ( − ω + ( n + 2) δ ) c ( n ) ( ω ) b ( n ) ( ω ) c ( n +2) † ( ω ) b ( n +2) † ( ω ) = √ κc ( n ) in ( ω ) √ γb ( n ) in ( ω ) √ κc ( n +2) † in ( ω ) √ γb ( n +2) † in ( ω ) , (A3)where the cavity and mechanical response functions read χ opt ( ω ) = [ κ/ − i ( ω + ∆)] − (cid:54) = χ c and χ m ( ω ) = [ γ/ − i ( ω − ω m )] − , respectively.Blocks with zero input will vanish in the steady state. Sincethe noise resides entirely in the zeroth Fourier component (itdescribes a stationary process), c ( n ) in = b ( n ) in = 0 , ∀ n (cid:54) = 0 . (A4)Although the result for general δ and ∆ is available [12],we focus on the simple and relevant case δ = − ∆ = ω m ,where (cid:18) b (0) ( ω ) b (2) † ( ω ) (cid:19) = J ( ω ) (cid:18) iG − χ opt ( ω )0 iG + χ opt ( ω ) (cid:19) (cid:18) √ γb in √ κ c in (cid:19) , (A5)with J ( ω ) = (cid:2) χ − m ( ω ) + G χ opt ( ω ) (cid:3) − . Now we include the noise of a BAE measurement as perSec. V. The measurement requires another cavity mode cou-pled to the mechanical oscillator. We define (cid:126)F = (cid:0) b b † c c † (cid:1) , (cid:126)F in = (cid:0) b in b † in c in c † in (cid:1) . (A6)From Eq. (A5), we can directly infer the elements (cid:32) χ (0)11 ( ω ) χ (0)13 ( ω ) χ (2)21 ( ω − δ ) χ (2)23 ( ω − δ ) (cid:33) = J ( ω ) (cid:18) iG − χ opt ( ω )0 iG + χ opt ( ω ) (cid:19) . (A7)Furthermore, the elements of χ are not all independent. Ingeneral, χ ( n ) ∗ ij ( ω ) = χ ( − n )¯ i ¯ j ( − ω ) , (A8)where ¯ i is the index for the operator that is the hermitian con-jugate of the operator indexed by i (so ¯1 = 2 , ¯3 = 4 , and viceversa ).We are only interested in the elements that connect the sys-tem operators with the mechanical bath, as those will deter-mine the susceptibility to the measurement noise. Fortunately,the only non-zero ones are χ (0)11 ( ω ) = J ( ω ) = χ (0) ∗ ( − ω ) . (A9)This gives the problem exactly the same structure as the casewith H rest = 0 , except that here χ m ( ω ) is replaced by J ( ω ) . Nevertheless, we give the general formula for the addedcontribution due to the measurement b ( n ) add ( ω ) = 2 √ κG (cid:88) m K ( n − m ) ( ω − nδ ) f ( m ) in ( ω ) , (A10a) b ( n ) † add ( ω ) = 2 √ κG (cid:88) m K ( − n + m ) ∗ ( − ω + nδ ) f ( m ) in ( ω ) , (A10b)where K ( n ) ( ω ) ≡ χ ( n )11 ( ω ) − χ ( n )12 ( ω ) . (A11)This gives (note similarity with Eq. (19b)) b ( n ) add = 2 √ κGJ ( ω − nδ ) f ( n ) in ( ω ) . (A12)and hermitian conjugate. The CR correction can be calculatedas before and looks very familiar (cf. Eq. (31)) S CR X δ X δ ( ω ) = κλ G (cid:0) | J ( ω + δ ) χ c ( ω + 2 δ ) | + | J ( − ω + δ ) χ c ( − ω + 2 δ ) | (cid:1) . (A13)The reason why it is so simple is that here b (0) ( ω ) only cou-ples to b in and not to b † in . This ceases to be the case for ∆ (cid:54) = − ω m , or in parametric squeezing in Appendix C. Appendix B: Dissipative bichromatic squeezing in the absenceof strong-coupling effects
As a relevant example for how a weakly coupled H rest canlead to effective parameters in the mechanical susceptibility,we consider the weak-coupling version of Appendix A. Weplace the red-detuned drive on the red sideband ∆ = − ω m ,but allow the other to vary. To second order in G ± , the effec-tive quantum Langevin equation is (in a frame rotating withfrequency δ ) ˙ b = − (cid:20) i ( ω m, eff − δ ) + Γ eff (cid:21) b + √ γb in + 2 i G√ κ s in , (B1)with [12] Γ eff = γ + 4 κ (cid:18) G − − G ε /κ (cid:19) , (B2a) ω m, eff = ω m + G ε ( κ/ + ε , (B2b) s in = G − d in + G + d † in G . (B2c)Here, ε ≡ δ − ω m ) and G ≡ (cid:113) G − − G . s in isa Bogoliubov rotation of the original optical bath opera-tors, and is therefore a squeezed, vacuum bath with nonzeroanomalous correlators, such as (cid:104) s in ( ω ) s in ( ω (cid:48) ) (cid:105) = − πδ ( ω + ω (cid:48) ) G + G − / G .The equivalent master equation (generalization of Ref. [17]to general drive detuning δ ) is ˙ ρ = − i [( ω m, eff − δ ) b † b, ρ ] + { Γ eff ( n eff + 1) D [ b ]+Γ eff n eff D [ b † ] + 4 G κ D [ β ] (cid:27) ρ, (B3)where D [ a ] ρ ≡ aρa † − ( a † aρ + ρa † a ) is the Lindblad su-peroperator.In order to include the measurement noise, it is easier towork in the laboratory frame, where ˙ b = − (cid:18) iω m, eff + Γ eff (cid:19) b + √ γb in + 2 √ κGf in + 2 i G√ κ s in . (B4)There is a subtlety here, as the anomalous averages of s in ,will become rotating in the transition from a rotating frameinto the laboratory frame. We can ignore this difficulty, be-cause we are only interested in the correction. RewritingEq. (B4) in Fourier components, we find that the only inde-pendent nonzero element of the susceptibility matrix is χ (0)11 ( ω ) = χ m, eff ( ω ) ≡ [Γ eff / − i ( ω − ω m, eff )] − . (B5)We could have obtained this from Eq. (A9) by setting G = 0 (weak coupling) and replacing ω m and γ by their modifiedvalues. Thus we can use the result Eq. (A13) with J → χ m, eff .Therefore, the results in the article are valid in the weak cou-pling case as well, with γ → Γ eff and ω m → ω m, eff .Optionally, we can remain in a rotating frame, but thatmeans we have to rotate the added measurement noise to ˜ f in ≡ e iδt f in . This implies ˜ f ( n ) in = f ( n − in . Solving the Langevinequation gives χ (0)11 ( ω ) = [Γ eff / − i ( ω − ω m, eff + δ )] − .Consulting Eq. (A10b), we find that the new set of Fouriercomponents ˜ b ( n ) add ( ω ) = b ( n − add ( ω ) . It is straightforward tocheck that in the end ˜ b add ( t ) = e iδt b add ( t ) .This result can be adapted to a wide variety of cases, aslong as they can be approximated by coupling the harmonicoscillators to baths only (and potentially modify its effectiveparameters). Appendix C: Parametric squeezing
Squeezing is induced naturally in a degenerate parametricamplifier [22]. Whilst it is limited to 3 dB of squeezing, itssolution has other interesting aspects. Here, H rest = ω m b † b + ( µb e iω m t + h.c. ) , (C1)where µ is the parametric driving strength. Without the BAEmeasurement ( G = 0 ), the quantum Langevin equation is ˙ b = − iω m b − iµb † e − iω m t − γ b + √ γb in . (C2)In Fourier components, (cid:18) χ − m ( ω − nω m ) iµ − iµ χ − m ( ω − nω m ) (cid:19) (cid:18) b ( n ) ( ω ) b ( n +2) † ( ω ) (cid:19) = √ γ (cid:18) δ n, b in δ n, b † in (cid:19) , (C3)where χ m ( ω ) ≡ [ γ/ − i ( ω − ω m )] − . Inverting leads to (cid:18) b ( n ) ( ω ) b ( n +2) † ( ω ) (cid:19) = A ( ω − nω m ) × (cid:18) χ − m ( ω − nω m ) − iµiµ χ − m ( ω − nω m ) (cid:19) (cid:18) √ γδ n, b in √ γδ n, b † in (cid:19) , (C4)where A ( ω ) ≡ (cid:2) χ − m ( ω − nω m ) − µ (cid:3) − . (C5)It is again enough to consider only one block (the othernonzero block is the hermitian conjugate) (cid:18) b (0) ( ω ) b (2) † ( ω ) (cid:19) = A ( ω ) (cid:18) χ − m ( ω ) − iµiµ χ − m ( ω ) (cid:19) (cid:18) √ γb in (cid:19) . (C6)To make the discussion as simple as possible, we choose δ = ω m for the measurement scheme. Then the same analysisas above can be applied to obtain χ (0)11 ( ω ) = A ( ω ) χ − m ( ω ) , (C7a) χ (2)12 ( ω − ω m ) = − iµA ( ω ) , (C7b) χ ( − ∗ ( − ω + 2 ω m ) = iµA ∗ ( ω ) . (C7c)In order to evaluate the added noise due to the BAEmeasurement, we have to add the noise Eq. (20) and useEq. (A10b). We find ( λ labels terms with CR origin) b ( n ) add ( ω ) = i √ κGχ c ( ω ) I in ( ω ) × iµA ( ω − δ ) if n = 3 iµλA ( ω + δ ) + A ( ω − δ ) χ − m ( ω − δ ) if n = 1 λA ( ω + δ ) χ − m ( ω + δ ) if n = − (C8)This reverts to Eq. (19b) when we set µ = 0 . Using Eq. (28)with the spectra calculated from Eq. (C8), we obtain the cor-rection to the quadrature spectrum. The resulting spectra haveterms rotating at multiples of δ . They could be measured bycoupling to another suitably driven cavity mode [12], but tendto be very small. Appendix D: Two-mode BAE
In this section, we consider another recently experimentallydemonstrated QND scheme [20, 21]. The goal here is to mea-sure a collective quadrature of two mechanical oscillators, inorder to eventually measure both quadratures of an externalforce with the same device. Whilst the overall theory is moregeneral [23], in the specific experiment we consider here, bothmechanical resonators couple to the same cavity and the prob-lem has the Hamiltonian [20] H = ( ω m + Ω) b † b + ( ω m − Ω) b † b + ω cav a † a + g ( b + b † ) a † a + g ( b + b † ) a † a + H drive + H diss (D1)Note that in Ref. [20] the mechanical oscillators have annihi-lation operators a and b , and the cavity c , whereas here, theyare b , b , and a , respectively. As before, bichromatic drivingleads to a coherent state ˆ a = 2 cos(Ω t )¯ a + ˆ d, (D2)where ¯ a is a real number, and (cid:104) ˆ d (cid:105) = 0 .Our generic solution is applicable, because Eq. (41) is ful-filled after linearizing, no matter which of the oscillators weput into H rest . For example, we could chose H rest,1 = ( ω m + Ω) b † b + g ( b + b † )2 cos(Ω t )¯ a ( d + d † ) , (D3)but the choice with → is equally valid. Our result for thebackaction on the mechanical oscillators (all of Sec. IV, par-ticularly Eq. (31)) thus applies to both resonators individually.The correction to the cavity spectrum is slightly more trickyto find, since the fluctuations in the two resonators are corre-lated, leading to cross terms. Instead of Eq. (4a) we have ˙ d = − κ d + √ κd in + 2 i cos(Ω t ) ( G x + G x ) , (D4)where x , are the position operators for the two oscillators,and G , ≡ g , ¯ a . x , are both given through Eq. (19b), butwith the respective parameters for each resonator.To calculate the optical spectrum, first note that there is across correlation. For m, n ∈ { , − } , we have (cid:90) d ω (cid:48) π (cid:68) x ( m )1 ( ω + mδ ) x ( n )2 ( ω (cid:48) ) (cid:69) = − κG G × | χ c ( ω + mδ ) | χ x ( ω ) χ x ( − ω − ( m + n ) δ ) , (D5) with χ x , defined analogously to χ x in Appendix G. Then S (0) d † d contains two copies of Eq. (32) (one for each resonator),in addition to G | χ c ( ω ) | Re (cid:104) S (0) x x ( ω − δ ) + S (0) x x ( ω + δ )+ S (2) x x ( ω − δ ) + S ( − x x ( ω + δ ) (cid:105) , (D6)where S x x are given through Eqs. (16) and (D5). Here, weused the property [ S ( n ) AB ( ω )] † = S ( − n ) B † A † ( ω + nδ ) [12].This demonstrates that our technique can also be applied formultiple modes coupled to the cavity, as long as [ H rest , I ] = 0 . Appendix E: Integration of the main peak
One goal of a BAE measurement of a mechanical oscillatorquadrature is to extract the quadrature variance. We show thatthe weight of the main peak of the optical spectrum S (0) d † d ( ω ) G | χ c ( ω ) | = S (0) xx ( ω − δ ) + S (0) xx ( ω + δ )+ S (2) xx ( ω − δ ) + S ( − xx ( ω + δ ) (E1)is a good measure for the quadrature variance. The weight ofits middle peak at ω = 0 is W = (cid:90) ∞−∞ d ω π S (0) xx ( ω ) + 2Re[ S (2) xx ( ω )] . (E2)Equation (E2) can be evaluated with the formulae in Ap-pendix G. Given the output spectrum κS (0) d † d ( ω ) , the weightcan be approximated by integrating from − ω m to ω m W ≈ (cid:10) X meas (cid:11) ≡ (cid:90) ω m − ω m d ω π S (0) d † out d out ( ω ) κG | χ c ( ω ) | . (E3) Appendix F: Full quadrature spectrum
Here we present the full expression for the spectrum of therotating quadrature X δ S (0) X δ X δ ( ω ) = γ (cid:0) | χ m ( ω + δ ) | ( n th + 1) + | χ m ( − ω + δ ) | n th (cid:1) + κG (cid:8) − | χ c ( ω ) | Re [ χ m ( ω + δ ) χ m ( − ω + δ )]+ | χ m ( ω + δ ) | (cid:0) | χ c ( ω ) | + | χ c ( ω + 2 δ ) | (cid:1) + | χ m ( − ω + δ ) | (cid:0) | χ c ( ω − δ ) | + | χ c ( ω ) | (cid:1)(cid:9) . (F1)It contains a part due to the thermal bath of the oscillator (terms with a γ ) and an optical part (terms with κ ). If δ = ω m , χ m ( − ω + δ ) = χ ∗ m ( ω + δ ) , and the negative term in the curly brackets cancels two resonant contributions in the second line.This cancellation makes the measurement BAE in RWA. The two terms left over cause the two side peaks to appear. Invertingthe sign of the first term in curly brackets would lead to the spectrum of the conjugate quadrature S (0) P δ P δ ( ω ) , which gets all theBA in RWA. The BA due to CR terms is the same in both. If the oscillator is weakly coupled to other degrees of freedom, thisformula is still approximately correct, using effective parameters Γ eff , ω m, eff and n eff .0 Appendix G: Optical spectrum
In this section we outline the derivation of the optical spectrum. We use the expression given in the main text (32). For H rest = 0 , the necessary correlators are (cid:90) d ω (cid:48) π (cid:68) x (0) ( ω ) x (0) ( ω (cid:48) ) (cid:69) = γ (cid:2) | χ m ( − ω ) | n th + | χ m ( ω ) | ( n th + 1) (cid:3) , (G1)and for m, n ∈ { , − } (cid:90) d ω (cid:48) π (cid:68) x ( m ) ( ω + mδ ) x ( n ) ( ω (cid:48) ) (cid:69) = − κG | χ c ( ω + mδ ) | χ x ( ω ) χ x ( − ω − ( m + n ) δ ) , (G2)where χ x ( ω ) ≡ χ m ( ω ) − λχ ∗ m ( − ω ) . This gives S (0) xx ( ω ) = γ (cid:2) | χ m ( − ω ) | n th + | χ m ( ω ) | ( n th + 1) (cid:3) + κG (cid:0) | χ x ( − ω ) | | χ c ( ω + δ ) | + | χ x ( ω ) | | χ c ( ω − δ ) | (cid:1) , (G3a) S (2) xx ( ω ) = − κG | χ c ( ω + δ ) | χ ∗ x ( − ω ) χ ∗ x ( ω + 2 δ ) , (G3b) S ( − xx ( ω ) = − κG | χ c ( ω − δ ) | χ x ( ω ) χ x ( − ω + 2 δ ) , (G3c)whence S (0) d † d ( ω ) can be constructed, using (cf. Eq. (32)) S ( n ) d † d ( ω ) = G χ c ( ω ) χ ∗ c ( ω + nδ ) (cid:104) S ( n ) xx ( ω − δ ) + S ( n ) xx ( ω + δ ) + S ( n +2) xx ( ω − δ ) + S ( n − xx ( ω + δ ) (cid:105) . (G4) [1] C. M. Caves, K. S. Thorne, R. W. P. Drever, V. D. Sandberg,and M. Zimmermann, Reviews of Modern Physics , 341(1980).[2] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, andR. J. Schoelkopf, Reviews of Modern Physics , 1155 (2010).[3] V. B. Braginsky, Y. I. Vorontsov, and K. S. Thorne, Science , 547 (1980).[4] B. P. Abbott and Others, Physical Review Letters , 061102(2016), arXiv:1602.03837.[5] M. F. Bocko and R. Onofrio, Reviews of Modern Physics ,755 (1996).[6] A. A. Clerk, F. Marquardt, and K. Jacobs, New Journal ofPhysics , 095010 (2008), arXiv:0802.1842.[7] J. B. Hertzberg, T. Rocheleau, T. Ndukum, M. Savva, A. A.Clerk, and K. C. Schwab, Nature Physics , 213 (2010).[8] J. Suh, A. J. Weinstein, C. U. Lei, E. E. Wollman, S. K. Steinke,P. Meystre, A. A. Clerk, and K. C. Schwab, Science , 1262(2014), arXiv:1312.4084.[9] F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado,and J. D. Teufel, Physical Review X , 041037 (2015),arXiv:1509.01629v1.[10] C. U. Lei, A. J. Weinstein, J. Suh, E. E. Wollman, A. Kronwald,F. Marquardt, A. A. Clerk, and K. C. Schwab, Physical ReviewLetters , 100801 (2016), arXiv:1605.08148.[11] E. S. Polzik and K. Hammerer, Annalen der Physik , A15(2014), arXiv:1405.3067.[12] D. Malz and A. Nunnenkamp, Physical Review A , 023803(2016), arXiv:1605.04749. [13] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Reviewsof Modern Physics , 1391 (2014).[14] C. Gardiner and P. Zoller, Quantum Noise: A Handbook ofMarkovian and Non-Markovian Quantum Stochastic Methodswith Applications to Quantum Optics , Springer Series in Syn-ergetics (Springer, 2004).[15] A. J. Weinstein, C. U. Lei, E. E. Wollman, J. Suh, A. Metel-mann, A. A. Clerk, and K. C. Schwab, Physical Review X ,041003 (2014), arXiv:1404.3242.[16] A. Mari and J. Eisert, Physical Review Letters , 213603(2009), arXiv:0911.0433.[17] A. Kronwald, F. Marquardt, and A. A. Clerk, Physical ReviewA , 063833 (2013).[18] J.-M. Pirkkalainen, E. Damsk¨agg, M. Brandt, F. Massel, andM. A. Sillanp¨a¨a, Physical Review Letters , 243601 (2015),arXiv:1507.04209.[19] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald,F. Marquardt, A. A. Clerk, and K. C. Schwab, Science , 952(2015), arXiv:1507.01662.[20] M. J. Woolley and A. A. Clerk, Physical Review A , 063846(2013), arXiv:1304.4059.[21] C. F. Ockeloen-Korppi, E. Damsk¨agg, J.-M. Pirkkalainen, A. A.Clerk, M. J. Woolley, and M. A. Sillanp¨a¨a, Physical ReviewLetters , 140401 (2016).[22] D. Walls and G. J. Milburn, Quantum Optics , edited by D. Wallsand G. J. Milburn (Springer Berlin Heidelberg, Berlin, Heidel-berg, 2008).[23] M. Tsang and C. M. Caves, Physical Review X2