Shelvingstyle QND phononnumber detection in quantum optomechanics
Abstract
We propose a new method for optomechanical quantum nondemolition (QND) detection of phonon number, based on a “shelving” style measurement. The scheme uses a twomode optomechanical system where the frequency splitting of the two photonic modes is nearresonant with the mechanical frequency. The combination of a strong optical drive and the underlying nonlinear optomechanical interaction gives rise to spinlike dynamics which facilitate the measurement. This approach allows phonon number measurement to be accomplished parametrically faster than in other schemes which are restricted to weak driving. The ultimate power of the scheme is controlled by the size of the single photon optomechanical cooperativity.
I Introduction
Quantum optomechanical systems, where photons interact with mechanical motion, have been the recent subject of intense experimental and theoretical activity Aspelmeyer2014a (). They have enabled the exploration of a wide range of effects, including the cooling of macroscopical mechanical modes to near their quantum ground states Teufel2011 (); Chan2011 (); Peterson2016 () and the generation of squeezed states of light Brooks2012 (); SafaviNaeini2013 (); Purdy2013 () and mechanical motion Wollman2015 (); Pirkkalainen2015 (); Lecocq2015 ().
A key unrealized goal in optomechanics is the quantum nondemolition (QND) measurement of mechanical phonon number, something that would directly reveal the energy quantization of the mechanics Santamore2004 (); Martin2007 (). The most promising proposals involve optomechanical systems where a single mechanical resonator couples to two photonic modes, as depicted schematically in Fig. 1a. Such systems can in principle be used for phononnumber detection both in the regime where the mechanical frequency is much smaller than the splitting between the optical modes Thompson2008 (); Jayich2008 (); Miao2009 (); Gangat2011 (); Yanay2016a (), and in the regime where it is comparable to this splitting Ludwig2012 (); Komar2013 (); BasiriEsfahani2012 (). Twocavity optomechanical systems have been realized experimentally using a variety of platforms Thompson2008 (); Paraiso2015 (); Toth2016 ().
In this work, we revisit the twomode system of Fig. (a)a, and focus on the regime where the mechanics are almost resonant with the frequency difference of the photonic modes, and where one photonic mode is strongly driven. This strong drive generates a coherent, manyphoton optomechanical interaction which cannot be treated perturbatively. We show that this regime allows a new approach to phonon number measurement, one which offers advantages over previous proposals (namely a much faster measurement). The regime we study differs from that of Ludwig et al. Ludwig2012 (), who considered more nonresonant interactions and weak driving, such that a pertubative treatment was possible. It also contrasts from the work of Kómár et al. Komar2013 (), who considered a perfectly resonant system and weak optical driving, but did not consider phonon number measurement physics.
Our measurement scheme is sketched in Fig. (b)b, and is conceptually similar to the highly successful electron shelving technique used in trappedion systems Dehmelt1968 (); Nagourney1986 (); Javanainen1986 (), as explained in Appendix A. The measurement starts by turning on the strong optical drive on the “primary” mode. This induces coherent oscillations between phonons and photons in the undriven “auxiliary” mode, oscillations which necessarily conserve the total number \hat{N}_{\rm tot} of auxmode and mechanical excitations. The remaining nonlinear interaction then allows one to use the output field of the driven optical mode to measure \hat{N}_{\rm tot} without changing its value  a true QND measurement. This is then effectively a measurement of mechanical phonon number, as without spurious dissipative effects, \hat{N}_{\rm tot} is simply equal to the initial phonon number.
As with other schemes, our approach still requires a relatively strong singlephoton optomechanical coupling (i.e. it must be large compared to the geometric mean of the two photonic damping rates) to resolve mechanical Fock states. However, it allows measurement rates that are parametrically larger than the perturbative regime considered in Ref. Ludwig2012, , or in the adiabatic regime of a positionsquared measurement considered in Ref. Thompson2008, ; Jayich2008, ; Miao2009, . This faster rate could be of significant utility in new experimental designs of twocavity optomechanical systems where one photonic mode is explicitly engineered to have extremely low dissipation Paraiso2015 (); Toth2016 (). In such systems, the ultimate limit to Fock state detection comes from the requirement that the measurement rate be faster than the heating rate of the mechanics due to their intrinsic dissipation. As we show, the relevant parameter controlling the utility of our scheme in this limit is the singlephoton cooperativity Aspelmeyer2014 (),
The remainder of the paper is arranged as follows. In Section II we introduce the theoretical model of our twomode optomechanical system, and discuss the parameter regime of interest. In Section III we introduce a mapping to an effective spin system that provides a convenient way to describe the system dynamics; we focus here on the ideal case where the only dissipation is due to the coupling to the input/output port used for the measurement. In Section IV we discuss the measurement protocol. Finally, in Section V we study the effects of the unwanted dissipation both in the aux mode and in the mechanics, and the limits that they impose on QND measurement.
II Model
The optomechanical system of interest (c.f. Fig. 1a) consists of two photonic modes and a single mechanical mode. Its coherent dynamics is described by the Hamiltonian
\begin{split}\displaystyle\hat{H}_{S}&\displaystyle=\hat{H}_{\rm opt}+\hat{H}_% {\rm m}+\hat{H}_{\rm int}.\end{split}  (1) 
The first term describes the two photonc eigenmodes in the absence of optomechanical coupling, and is (\hbar=1 throughout)
\begin{split}\displaystyle\hat{H}_{\rm opt}&\displaystyle=\mathopen{}% \mathclose{{}\left(\omega_{c}J}\right)\hat{a}_{+}^{\dagger}\hat{a}_{+}+% \mathopen{}\mathclose{{}\left({\omega_{c}+J}}\right)\hat{a}_{}^{\dagger}\hat{% a}_{},\end{split}  (2) 
where \hat{a}_{\pm} are annihilation operators for the two photonic modes (frequencies \omega_{c}J and \omega_{c}+J respectively). We call these the primary (\hat{a}_{+}) and auxiliary (\hat{a}_{}) modes.
The second term in Eq. (1) describes a mechanical mode with frequency \omega_{\rm m} and annihilation operator \hat{b}:
\hat{H}_{\rm m}=\omega_{\rm m}\hat{b}^{\dagger}\hat{b}=\mathopen{}\mathclose{{% }\left({2J\delta\Omega}}\right)\hat{b}^{\dagger}\hat{b}.  (3) 
We have introduced the parameter \delta\Omega=2J\omega_{m} which represents the mismatch between the photonic mode splitting and the mechanical frequency. Similar to reference Ludwig2012 (), we will be interested in the regime where the mechanical frequency is similar in magnitude to 2J, and hence \delta\Omega\ll J.
Finally, the last term in Eq. 1 describes the optomechanical coupling, which takes the form of a mechanicallymediated tunneling interaction between the two modes:
\begin{split}\displaystyle\hat{H}_{\rm int}&\displaystyle=g\mathopen{}% \mathclose{{}\left({\hat{b}+\hat{b}^{\dagger}}}\right)\mathopen{}\mathclose{{}% \left({\hat{a}_{}^{\dagger}\hat{a}_{+}+\hat{a}_{+}^{\dagger}\hat{a}_{}}}% \right),\end{split}  (4) 
where g is the singlephoton optomechanical coupling amplitude. For the specific case of a membraneinthemiddle style optomechanical system Thompson2008 (), a detailed derivation of this model is given in Refs. Jayich2008, ; Cheung2011, ; its limitations are discussed in Refs. Cheung2011, ; Xuereb2012, .
We next use a standard treatment to describe the effects of dissipation on the system Gardiner2004 (); Clerk2010 (). We assume that each optical mode is coupled to an independent zerotemperature Markovian reservoir, giving rise to an amplitude damping rate of \kappa_{\pm}/2 for \hat{a}_{\pm}, respectively. The mechanics are similarly coupled to a thermal Markovian reservoir characterized by a thermal occupancy n_{\rm th}, leading to an amplitude damping rate of \gamma/2 for \hat{b}. Note that for a membraneinthemiddle type system, treating the optical modes as seeing independent reservoirs can be problematic Yanay2016a (), as it neglects a possible dissipationinduced coupling of the optical normal modes. Such a coupling will not play a role in our scheme, as it is suppressed by the optical mode splitting 2J, while the main optomechanical coupling is resonant.
Our measurement scheme will involve applying a strong coherent drive to the primary, \hat{a}_{+}, mode near its resonance frequency \omega_{c}J. For simplicity, we assume the case of a perfectly resonant drive, as qualitatively similar results are obtained for a detuned drive. In this case, it is useful to move to a rotating frame, and also make a standard displacement transformation on the optical modes:
\begin{gathered}\displaystyle\hat{a}_{}\to e^{i\mathopen{}\mathclose{{}\left% ({\omega_{c}+J\delta\Omega/2}}\right)t}\hat{c}_{},\qquad\hat{b}\to e^{i% \mathopen{}\mathclose{{}\left({\omega_{\rm m}+\delta\Omega/2}}\right)t}\hat{b}% ,\\ \displaystyle\hat{a}_{+}\to e^{i\mathopen{}\mathclose{{}\left({\omega_{c}J}}% \right)t}\mathopen{}\mathclose{{}\left({\langle{\hat{a}_{+}}\rangle_{g=0}+\hat% {c}_{+}}}\right),\quad G=g\langle{\hat{a}_{+}}\rangle_{g=0}.\end{gathered}  (5) 
Here, \langle{\hat{a}_{+}}\rangle_{g=0} is average primarymode amplitude induced by the drive in the absence of optomechanical coupling, and \hat{c}_{+} is the displaced annihilation operator for this mode in the rotating frame.
In this frame we find \hat{H}_{S}=\hat{H}_{\rm eff}+\hat{H}_{\rm NR} where
\displaystyle\begin{split}\displaystyle\hat{H}_{\rm eff}&\displaystyle=\tfrac{% \delta\Omega}{2}\begin{pmatrix}\hat{c}_{}^{\dagger}\hat{c}_{}\hat{b}^{% \dagger}\hat{b}\end{pmatrix}G\begin{pmatrix}\hat{c}_{}^{\dagger}\hat{b}+\hat% {b}^{\dagger}\hat{c}_{}\end{pmatrix}\\ &\displaystyle\quadg\begin{pmatrix}\hat{b}^{\dagger}\hat{c}_{+}^{\dagger}\hat% {c}_{}+\hat{c}_{}^{\dagger}\hat{c}_{+}\hat{b}\end{pmatrix}\end{split}  (6)  
\displaystyle\hat{H}_{\rm{NR}}  \displaystyle=e^{i\mathopen{}\mathclose{{}\left(4J\delta\Omega}\right)t}\hat% {b}^{\dagger}\hat{c}_{}^{\dagger}\begin{pmatrix}G+g\hat{c}_{+}\end{pmatrix}+% \operatorname{h.c.}  (7) 
Finally, we assume that the optical mode gap 2J is sufficiently large that we can safely make a rotating wave approximation and drop the nonresonant interactions described by \hat{H}_{\rm{NR}}. This requires in practice J\gg\delta\Omega,G,g,\kappa_{+}. We are left with the interactions depicted in Fig. 1a: a beam splitter interaction between the mechanics and the auxiliary mode, and a nonlinear interaction involving all three modes.
Note that in the absence of driving, i.e. G=0, the photonic modes are only driven by vacuum noise, and the nonlinear optomechanical interaction in Eq. 6 effectively vanishes. Therefore, without driving the mechanics are completely decoupled from the photonic modes. We imagine that before the measurement is turned on by driving the primary photonic mode, the mechanics are first prepared close to the ground state, e.g. via standard cavity cooling techniques. The scheme we describe in what follows then allows one to measure the mechanical phonon number.
III Dynamics and effective spin description
III.1 Conserved excitation number
The nonlinear interaction in Eq. 6 has the form of a standard threewave mixing Hamiltonian, describing e.g. a nondegenerate parametric amplifier with a dynamical pump mode. It has three conserved excitation numbers (typically known as ManleyRowe constants of motion Manley1956 ()) reflecting the SU(1) and SU(1,1) symmetries of the interaction. Including the beamsplitter interaction reduces the symmetry of the Hamiltonian, and \hat{H}_{\rm eff} has only a single conserved excitation number \hat{N}_{\rm tot}:
\hat{N}_{\rm tot}=\hat{b}^{\dagger}\hat{b}+\hat{c}_{}^{\dagger}\hat{c}_{},% \hskip 17.071654pt\mathopen{}\mathclose{{}\left[\hat{N}_{\rm tot},\hat{H}_{\rm eff% }}\right]=0  (8) 
Note that as \hat{N}_{\rm tot} is independent of the primary mode, it continues to be conserved even if this mode is damped.
In what follows, we first analyze the ideal case where the only dissipation in the system is due to its coupling to the inputoutput port used for the measurement. This corresponds to \kappa_{+}\neq 0, while \gamma=\kappa_{}\simeq 0. Starting with this limit will let us present a clear picture of how our proposed phonon number measurement operates: the conservation of \hat{N}_{\rm tot} in this limit is at the heart of the shelving scheme depicted in Fig. 1b. In this ideal case, we find it is always possible to determine the initial mechanical phonon number by simply measuring for long enough; this is akin to standard shelving measurements.
The regime \kappa_{+}\gg\gamma,\kappa_{} is also relevant experimentally. While mechanical damping rates \gamma are typically orders of magnitude smaller than photonic damping rates, it might first seem surprising that one could achieve \kappa_{}\ll\kappa_{+}. However, recent experiments with microwavecircuit optomechanics have achieved precisely this kind of regime through clever experimental design.
We analyze the effects of nonzero \gamma,\kappa_{} in Section IV. As expected, they limit the maximum possible measurement time, putting constraints on whether phonon number detection is achievable.
III.2 Mapping to an effective spin system
To understand the dynamics of our system, it is helpful to make the conservation of \hat{N}_{\rm tot} explicit by representing the mechanical and auxmodes by an effective spin (a Schwinger boson representation in reverse). We thus introduce the effective (dimensionless) angular momentum operators \boldsymbol{\hat{J}} via
\begin{gathered}\displaystyle\hat{J}_{+}=\hat{b}^{\dagger}\hat{c}_{}\qquad% \hat{J}_{}=\hat{c}_{}^{\dagger}\hat{b}\\ \displaystyle\hat{J}_{\pm}\equiv\hat{J}_{x}\pm i\hat{J}_{y}\qquad\hat{J}_{z}=% \tfrac{1}{2}\mathopen{}\mathclose{{}\left({\hat{b}^{\dagger}\hat{b}\hat{c}_{% }^{\dagger}\hat{c}_{}}}\right).\end{gathered}  (9) 
These operators satisfy standard angular momentum commutation relations. A direct calculation shows:
\hat{J^{2}}=\frac{\hat{N}_{\rm tot}}{2}\mathopen{}\mathclose{{}\left(\frac{% \hat{N}_{\rm tot}}{2}+1}\right).  (10) 
The conservation of \hat{N}_{\rm tot} now just corresponds to our effective spin having a fixed length.
With these operators, the effective Hamiltonian is
\begin{split}\displaystyle\hat{H}_{\rm eff}&\displaystyle=\boldsymbol{B}\cdot% \boldsymbol{\hat{J}}g\begin{pmatrix}\hat{J}_{}\hat{c}_{+}+\hat{c}_{+}^{% \dagger}\hat{J}_{+}\end{pmatrix}.\end{split}  (11) 
where the effective magnetic field is
\boldsymbol{B}=2G\mathbf{e}_{x}+\delta\Omega\mathbf{e}_{z}\equiv B\mathbf{e}_{B}  (12) 
with \mathbf{e}_{x},\mathbf{e}_{y},\mathbf{e}_{z} the unit vectors in the respective directions. As the \hat{c}_{+} mode is damped, we see that \hat{H}_{\rm eff} takes the form of a generalized spinboson model. The quadratic terms in the original version of \hat{H}_{\rm eff} now correspond to a Zeeman field on our spin. The order g nonlinear term in \hat{H}_{\rm eff} now corresponds to a systembath coupling where a bath boson can be created while adding angular momentum in the z direction to the spin.
The HeisenbergLangevin equations of motion take the form
\displaystyle\dot{\hat{c}}_{+}=\tfrac{\kappa_{+}}{2}\hat{c}_{+}+\sqrt{\kappa_% {+}}\hat{c}_{\rm in}+ig\hat{J}_{+}  (13)  
\displaystyle\boldsymbol{\dot{\hat{J}}}=\boldsymbol{\hat{J}}\times\mathopen{}% \mathclose{{}\left({\boldsymbol{B}+g\mathopen{}\mathclose{{}\left[{\mathbf{e}_% {}\hat{c}_{+}+\mathbf{e}_{+}\hat{c}_{+}^{\dagger}}}\right]}}\right)  (14) 
where \mathbf{e}_{\pm}=\mathbf{e}_{x}\pm i\mathbf{e}_{y}. The operator \hat{c}_{\rm in} describes the vacuum fluctuations entering the system from the inputoutput line coupled to the principal + mode (i.e. from the inputoutput waveguide that will be used for measurement); it has zero mean and correlation functions
\langle{\hat{c}_{\rm in}^{\dagger}\mathopen{}\mathclose{{}\left({t}}\right)% \hat{c}_{\rm in}\mathopen{}\mathclose{{}\left({t^{\prime}}}\right)}\rangle=0% \qquad\langle{\hat{c}_{\rm in}\mathopen{}\mathclose{{}\left({t}}\right)\hat{c}% _{\rm in}^{\dagger}\mathopen{}\mathclose{{}\left({t^{\prime}}}\right)}\rangle=% \delta\mathopen{}\mathclose{{}\left({tt^{\prime}}}\right).  (15) 
When the optomechanical drive is initially turned on, \hat{N}_{\rm tot}=\hat{n}_{b}, where \hat{n}_{b}\equiv\hat{b}^{\dagger}\hat{b} is the initial mechanical phonon number. From Eq. (10), we see that the phonon number manifests itself in the dynamics by directly setting the size of the effective spin \boldsymbol{\hat{J}}. The nonlinear interaction creates photons in the principal mode \hat{c}_{+} in a way that depends on \hat{J}_{\pm}, and hence on the magnitude of \boldsymbol{\hat{J}}. Thus, by monitoring the output field from this mode, one can determine the initial phonon number. We stress that the “backaction” from the this measurement only alters the direction our effective spin points in, but not its size. The measurement is thus QND.
For easy reference, a synopsis of the mapping from bosons to spins is provided in Table 1.
Spin term  Optomechanical Equivalent  

\hat{J}_{z}  \tfrac{\hat{b}^{\dagger}\hat{b}\hat{c}_{}^{\dagger}\hat{c}_{}}{2}  The zprojection of the effective spin represents the difference between the number of phonons and number of auxmode quanta. 
\hat{J}_{+}, \hat{J}_{}  \hat{b}^{\dagger}\hat{c}_{}, \hat{c}_{}^{\dagger}\hat{b}  
B_{x}  2G  The manyphoton optomechanical interaction G exchanges photons and phonons, and hence acts like a magnetic field in \mathbf{e}_{x} direction. 
B_{z}  \delta\Omega  The mismatch between the photonic mode splitting and the mechanical frequency acts as an \mathbf{e}_{z} field. 
III.3 Dynamics in the limit g\ll\kappa_{+}
We consider the system dynamics in the typical to experiment, {\kappa_{+}\gg g}. We begin by solving Eq. 13 to find
\begin{split}\displaystyle\hat{c}_{+}\mathopen{}\mathclose{{}\left({t}}\right)% &\displaystyle=e^{\tfrac{1}{2}\kappa_{+}t}\hat{c}\mathopen{}\mathclose{{}% \left({0}}\right)+\hat{\zeta}\mathopen{}\mathclose{{}\left({t}}\right)\\ &\displaystyle\quad+ig\int_{0}^{t}{d\tau}\;e^{\tfrac{1}{2}\kappa_{+}\tau}\hat% {J}_{+}\mathopen{}\mathclose{{}\left({t\tau}}\right).\end{split}  (16) 
where \hat{\zeta}\mathopen{}\mathclose{{}\left({t}}\right)=\sqrt{\kappa}\int_{0}^{t}% {d\tau}\;e^{\tfrac{1}{2}\kappa_{+}\mathopen{}\mathclose{{}\left({t\tau}}% \right)}\hat{c}_{\rm in}\mathopen{}\mathclose{{}\left({\tau}}\right).
We see that \hat{c}_{+}(t) depends on the behaviour of \boldsymbol{\hat{J}} at earlier times, and that the range of this memory effect is \kappa_{+}^{1}. We can solve for the spin dynamics during this short interval to a good approximation by ignoring the effects of g, as
\boldsymbol{\hat{J}}\mathopen{}\mathclose{{}\left({t\tau}}\right)=e^{i% \boldsymbol{B}\cdot\boldsymbol{\hat{J}}\tau}\mathopen{}\mathclose{{}\left({% \boldsymbol{\hat{J}}\mathopen{}\mathclose{{}\left({t}}\right)+O\mathopen{}% \mathclose{{}\left({g\tau}}\right)}}\right)e^{i\boldsymbol{B}\cdot\boldsymbol% {\hat{J}}\tau}.  (17) 
As long as \tau is at most \sim\kappa_{+}, we can drop the O\mathopen{}\mathclose{{}\left({g\tau}}\right) term as it is proportional to the small parameter
\eta=2g/\kappa_{+}.  (18) 
With this approximation, the solution for \hat{c}_{+}(t) takes the simpler form:
\hat{c}_{+}\mathopen{}\mathclose{{}\left({t}}\right)=e^{\tfrac{1}{2}\kappa_{+% }t}\hat{C}_{0}+\hat{\zeta}\mathopen{}\mathclose{{}\left({t}}\right)+i\eta\hat{% J}_{c}\mathopen{}\mathclose{{}\left({t}}\right)+O\mathopen{}\mathclose{{}\left% ({\eta}}\right)^{2}.  (19) 
Here, the nonHermitian operator \hat{J}_{c} is linear in the components of \boldsymbol{\hat{J}}, and given by
\begin{gathered}\displaystyle\hat{J}_{c}=\mathopen{}\mathclose{{}\left[{\frac{% B_{x}\boldsymbol{B}+\tfrac{\kappa_{+}}{2}\mathbf{e}_{+}\times\boldsymbol{B}+% \mathopen{}\mathclose{{}\left({\tfrac{\kappa_{+}}{2}}}\right)^{2}\mathbf{e}_{+% }}{\mathopen{}\mathclose{{}\left({\tfrac{\kappa_{+}}{2}}}\right)^{2}+B^{2}}}}% \right]\cdot\boldsymbol{\hat{J}}.\end{gathered}  (20) 
The initial transient behaviour of \hat{c}_{+} is described by {\hat{C}_{0}=\hat{c}_{+}\mathopen{}\mathclose{{}\left({0}}\right)i\eta\hat{J}% _{c}\mathopen{}\mathclose{{}\left({0}}\right)}.
Equation 19 implies that the primary mode (and its output optical field) will be sensitive to components of \boldsymbol{\hat{J}} in a manner that depends on both \boldsymbol{B} and \kappa_{+}. In the extreme Markovian limit {\kappa_{+}\gg B} we find \hat{J}_{c}\to\hat{J}_{+}, implying that the \hat{c}_{+} is able to measure both transverse components of our effective spin. In contrast, in the limit B\gg\kappa_{+}, the \hat{c}_{+} mode is only able to respond to the nonrotating component of the spin, and hence \hat{J}_{c}\to(B_{x}/B)\mathbf{e}_{B}\cdot\hat{J}.
It is also instructive to consider the dynamics in terms of the reduced density matrix describing the spin. The interaction picture density matrix is defined in terms of the Schrödingerpicture density matrix \hat{\rho}_{J}=\operatorname{Tr}_{\hat{c}_{+}}\mathopen{}\mathclose{{}\left[{% \hat{\rho}}}\right] via
\tilde{\rho}_{J}(t)=e^{i\boldsymbol{B}\cdot{\boldsymbol{\hat{J}}}t}\hat{\rho}% _{J}(t)e^{i\boldsymbol{B}\cdot{\boldsymbol{\hat{J}}}t}.  (21) 
We again consider the g\ll\kappa_{+} limit, and derive a weakcoupling master equation, retaining terms to leading order in \eta. For simplicity, we also consider the typical limit where B\gg\mathopen{}\mathclose{{}\left({\eta g=2g^{2}/\kappa_{+}}}\right), allowing us to drop rapidly oscillating terms (i.e. make a secular approximation).
We use \hat{J}_{\parallel}=\mathbf{e}_{B}\cdot\boldsymbol{\hat{J}} to denote the projection of the spin along the axis of the effective magnetic field \boldsymbol{B} , and use \hat{J}_{\pm}^{\prime} to denote raising and lowering operators for \hat{J}_{\parallel}. Defining \mathcal{D}\mathopen{}\mathclose{{}\left[{\hat{a}}}\right]\cdot{\hat{\rho}}=% \hat{a}\hat{\rho}\hat{a}^{\dagger}\tfrac{1}{2}\mathopen{}\mathclose{{}\left({% \hat{a}^{\dagger}\hat{a}\hat{\rho}+\hat{\rho}\hat{a}^{\dagger}\hat{a}}}\right) as the standard Lindblad superoperator, we find
\begin{split}\displaystyle\dot{\tilde{\rho}}_{J}&\displaystyle=i\mathopen{}% \mathclose{{}\left[{\tilde{\rho}_{J},\hat{H}_{1}}}\right]+\\ &\displaystyle\mathopen{}\mathclose{{}\left(\Gamma_{\varphi}\mathcal{D}% \mathopen{}\mathclose{{}\left[{\hat{J}_{\parallel}}}\right]+\Gamma_{+}\mathcal% {D}\mathopen{}\mathclose{{}\left[{\hat{J}_{+}^{\prime}}}\right]+\Gamma_{}% \mathcal{D}\mathopen{}\mathclose{{}\left[{\hat{J}_{}^{\prime}}}\right]}\right% )\cdot{\tilde{\rho}_{J}}.\end{split}  (22) 
This has the form of a standard master equation for a spin coupled to a bath. The first term describes a modification of the coherent dynamics, described by
\hat{H}_{1}=\tfrac{g^{2}B}{B^{2}+\mathopen{}\mathclose{{}\left({\tfrac{\kappa_% {+}}{2}}}\right)^{2}}\mathopen{}\mathclose{{}\left[{\tfrac{B_{z}}{B}\hat{J}_{% \parallel}\tfrac{B^{2}+B_{z}^{2}}{2B^{2}}\mathopen{}\mathclose{{}\left({% \boldsymbol{\hat{J}}^{2}\hat{J}_{\parallel}^{2}}}\right)}}\right].  (23) 
The remaining terms describe dephasing (rate \Gamma_{\varphi}), and spin raising / lowering transitions along the axis \mathbf{e}_{B} (rates \Gamma_{\pm}). The dissipative rates are
\begin{gathered}\displaystyle\Gamma_{\varphi}=\frac{B_{x}^{2}}{B^{2}}\frac{4g^% {2}}{\kappa_{+}},\qquad\Gamma_{\pm}=\frac{\mathopen{}\mathclose{{}\left(B\pm B% _{z}}\right)^{2}}{2B^{2}}\frac{g^{2}\kappa_{+}}{B^{2}+\mathopen{}\mathclose{{}% \left({\tfrac{\kappa_{+}}{2}}}\right)^{2}}.\end{gathered}  (24) 
Note that unless B_{z}=B, the vacuum noise driving the \hat{a}_{+} mode is able to cause spin flips upwards in energy: this is possible because we are strongly driving the \hat{a}_{+} mode, and can be viewed as being analogous to the phenomenon of “quantum heating” Dykman2011 (); Peano2010 (); Lemonde2015 (). Note also that we always assume that the drive frequency \omega_{c}J is much larger than G,\delta\Omega.
From Eq. 22 we find that the steady state of the spin is simply a thermal state,
\tilde{\rho}_{J}\to e^{B\hat{J}_{\parallel}/T_{\rm eff}}/\operatorname{Tr}% \mathopen{}\mathclose{{}\left[{e^{B\hat{J}_{\parallel}/T_{\rm eff}}}}\right]  (25) 
with an effective temperature
B/T_{\rm eff}=2\ln\mathopen{}\mathclose{{}\left({\frac{B+B_{z}}{BB_{z}}}}% \right).  (26) 
Again, unless B_{z}=B, the steady state will not correspond to the spin being in its ground state (due to the previously mentioned quantum heating physics). Note that in the case where the effective field is completely in the x direction (i.e. the detuning \delta\Omega vanishes), our master equation describes an effective infinite temperature, and the spin will be completely depolarized in the steady state.
Figure 2 illustrates the dissipative spin dynamics described above (see Appendix C for details on numerics).
III.4 QND nature of the measurement
As \hat{N}_{\rm tot} is a conserved quantity, it is clear that our protocol allows for a QND measurement of \hat{N}_{\rm tot}: this quantity commutes with the system Hamiltonian, even when we include the coupling \kappa_{+} to the inputoutput port used for the measurement.
Perhaps less obvious is that this also corresponds to a QND measurement of the phonon number in the mechanics. As already discussed, when the control drive is first turned on, the value of \hat{N}_{\rm tot} corresponds to the initial mechanical phonon number \hat{n}_{b}. For a true QND measurement of phonon number, one would like the final phonon number to return to this value after the measurement is complete. This is easily achieved by turning off the large control drive applied to the primary mode. Once this drive is off, G=0, the nonlinear interaction (last term in Eq. 6) keeps acting until \hat{n}_{b}=\hat{N}_{\rm tot}, i.e. until all the original phonons are returned to the mechanical resonator.
IV Description of the measurement output
IV.1 Homodyne measurement and measurement time in the absence of spurious dissipation
Having discussed its dynamics, we now turn to how information on the spin system (and hence the initial photon number) emerges in the output field leaving the primary phononic mode \hat{c}_{+}. Eq. 19 tells us that the output field will be sensitive to particular components of the spin vector \boldsymbol{\hat{J}}. Recall that in our original lab frame, the primary mode is driven at its resonance frequency \omega_{+}=\omega_{c}J. Information on our effective spin will be optimally encoded in one quadrature of the primary mode output field. We thus consider a homodyne measurement which directly probes this optimal output field quadrature, again starting with the ideal case where \gamma=\kappa_{}=0.
From standard inputoutput theory (and working in the interaction picture defined by Eq. (5)), the cavity output field leaving the + mode is given by
\hat{c}_{\rm out}(t)=\hat{c}_{\rm in}(t)\sqrt{\kappa_{+}}\hat{c}_{+}(t),  (27) 
and the phonon number can be obtained from measuring an output quadrature,
\hat{X}^{\alpha}_{\rm out}(t)=\mathopen{}\mathclose{{}\left(e^{i\alpha}\hat{c}% _{\rm out}(t)+e^{i\alpha}\hat{c}_{\rm out}^{\dagger}(t)}\right)/\sqrt{2},  (28) 
where the angle \alpha parameterizes the choice of quadrature.
Using the result in Eq. 19, we find that the average value of this output quadrature is given by:
\langle{\hat{X}_{\rm out}^{\alpha}(t)}\rangle=\sqrt{\frac{8g^{2}}{\kappa_{+}}}% \mathrm{Im}\mathopen{}\mathclose{{}\left[{e^{i\alpha}\langle{\hat{J}_{c}(t)}% \rangle}}\right].  (29) 
In the steady state, our spin relaxes to a thermal state, and in general the value of \langle{\hat{J}_{c}(t)}\rangle will be nonzero and reflect the overall size of the spin. It will thus be related to the initial mechanical phonon number, as illustrated in Fig. 3.
One finds from Eqs. 25 and 20 that in the steadystate \langle{\hat{J}_{c}(t)}\rangle is purely real, implying that the optimal choice of quadrature corresponds to an angle \alpha=\pi/2. This corresponds to measuring the phase quadrature of the output light leaving the primary mode.
The simplest way to extract the phonon number is to first turn on the measurement by driving the primary photonic mode, and then integrating the output homodyne current (which is proportional to \hat{X}_{\rm out}^{\alpha}) for a time t. For sufficiently long integration times, the signal of phonon number in the average homodyne current will be resolvable above the noise. For weak couplings, the noise in the homodyne current is simply due to vacuum noise and is Gaussian. As a result, standard calculations Clerk2010 (); Jayich2008 () indicate that the mechanical Fock state \hat{n}_{b}=n+1 can be distinguished from Fock state \hat{n}_{b}=n when the measurement time is sufficiently long, namely
t>\epsilon^{1}\mathopen{}\mathclose{{}\left(\tfrac{1}{j_{n}}}\right)^{2}\tau_% {\rm meas}  (30) 
where
\tau_{\rm meas}=\kappa_{+}/g^{2}.  (31) 
and the prefactor is determined by
\begin{split}\displaystyle j_{n}=4\mathrm{Im}\Big{[}e^{i\alpha}\Big{(}&% \displaystyle\langle{\hat{J}_{c}\mathopen{}\mathclose{{}\left({t\to\infty}}% \right)}\rangle_{n_{b}=n+1}\\ &\displaystyle\langle{\hat{J}_{c}\mathopen{}\mathclose{{}\left({t\to\infty}}% \right)}\rangle_{n_{b}=n}\Big{)}\Big{]}.\end{split}  (32) 
As one can see from Fig. 3, for a broad range of G/\delta\Omega, the prefactor 1/j_{n}^{2} is of order unity, and the timescale of the measurement is set by \tau_{\rm meas} as given in Eq. 31.
We stress that the time scale required to distinguish between different mechanical Fock states in our setup is parametrically shorter than that obtained in the proposal of Ludwig et al. Ludwig2012 (). Their scheme requires a measurement time which scales as t\sim\mathopen{}\mathclose{{}\left({\delta\Omega/G}}\right)^{2}\tau_{\rm meas} to resolve two adjacent Fock states. However, the setup requires an interaction that is sufficiently nonresonant such that hybridization of phonons and photons is weak, and such that a perturbative treatment is applicable. In practice, this requires both the conditions g\ll\delta\Omega and G\ll\delta\Omega, yielding a measurement time which scales inversely with a small parameter. Our approach is not constrained in the same way, and can operate much closer to perfect resonance, where G\sim\delta\Omega. As a result, for the same value of g, our approach yields greatly enhanced measurement speed. This in turn implies that an equally rapid measurement can be achieved with a smaller optomechanical coupling: where Ludwig et al. use g=3\kappa_{+} to resolve phonon number states after t=50\kappa_{+}^{1} (see Ref. Ludwig2012, Fig. 3), Eq. 31 implies our shelving protocol could do the same in a system where g=0.3\kappa_{+}.
IV.2 Optimal parameters
Our system has a fairly large number of tunable parameters. In particular, both nonzero components of the effective field \boldsymbol{B}=(G,0,\delta\Omega) can be selected for optimal performance. We remind the reader that G is the manyphoton coupling strength and \delta\Omega the mismatch between the mechanical frequency and photonic mode splitting.
A key parameter to optimize the size of the signal in the longtime homodyne current which encodes the initial phonon number. From Eq. 29, this is directly determined by the steadystate value of \langle{\hat{J}_{c}}\rangle. In steady state, we obtain from Eqs. 25 and 20,
\langle{\hat{J}_{c}}\rangle\to\frac{B_{x}}{B}\langle{\langle{\hat{J}_{% \parallel}}\rangle}\rangle_{T_{\rm eff}}  (33) 
where \langle{\langle{\hat{J}_{\parallel}}\rangle}\rangle_{T_{\rm eff}} is the thermal average of \hat{J}_{\parallel}=\boldsymbol{e}_{B}\cdot\boldsymbol{\hat{J}} taken at a temperature T_{\rm eff} given by Eq. 26.
It follows from these expressions that to have an appreciable signature of the phonon number in the steadystate value of \langle{\hat{J}_{c}}\rangle we need B_{x} and B_{z} to be comparable. A finite value of B_{x} (i.e. manyphoton coupling G) is needed to ensure that the steady state has a transverse component that can be detected by the primary photonic mode, as this mode couples to \hat{J}_{+}. A nonzero value of B_{z} (i.e. mismatch between mechanical frequency and photonic mode splitting) is needed so that the effective temperature describing the steadystate is finite and the spin is not in a completely depolarized state.
To make the above discussion more concrete, consider the simple case where we have one phonon intially, and our effective spin corresponds to a spin 1/2. There,
\langle{\hat{J}_{c}\mathopen{}\mathclose{{}\left({t\to\infty}}\right)}\rangle_% {n_{b}=1}=\frac{B_{x}B_{z}}{B^{2}+B_{z}^{2}}=\frac{G\delta\Omega}{2G^{2}+% \delta\Omega^{2}}.  (34) 
Figure 3 shows the steady state value of \langle J_{c}\rangle for other values of initial phonon number. We find the optimal ratio, for a system with one or a few initial phonons, is
G/\delta\Omega\approx 3/4.  (35) 
While the steadystate homodyne signal is not sensitive to the size of \mathbf{B} (but rather only to the ratio B_{z}/B_{x}), the overall magnitude of the effective field is nonetheless important to our scheme. If B\gg\kappa_{+}, then it follows from Eqs. 24 and 22 that the relaxation of the spin to its steady state is greatly suppressed. In this limit, there is simply very little density of states in the primary + mode for transitions that involve changing the energy of our spin. We find that it is generally optimal to have the spin reach the steady state on a time scale much shorter than \tau_{\rm meas}, requiring B\ll\kappa_{+}. That is, one should avoid the manyphoton optomechanical strong coupling regime.
Finally, to have the simple shelving dynamics depicted in Fig. 1b, the coherent oscillations between c_{} and b should be underdamped, leading to the requirement (G\sim B)\gg g^{2}/\kappa_{+}. This condition is easy to meet in experiment given the typical weakness of singlephoton optomechanical coupling strengths. Putting these conditions together, we find that the optimal regime for phonon number measurement requires:
g^{2}/\kappa_{+}\ll\mathopen{}\mathclose{{}\left(B=\sqrt{G^{2}+(\delta\Omega)^% {2}}}\right)\ll\kappa_{+}.  (36) 
An alternative approach to the measurement which results in faster measurement rates but is more complicated to implement experimentally is discussed in Appendix B.
IV.3 Resolving Fock states
We next present numerical results for the measurement output of our scheme based on simulations of the full master equation dynamics, including \hat{c}_{+}. We construct an estimator for phonon number from the integrated homodyne current:
\begin{split}&\displaystyle\hat{n}_{b}^{\rm meas}\mathopen{}\mathclose{{}\left% ({t}}\right)\equiv\\ &\displaystyle\qquad\mathopen{}\mathclose{{}\left[{\sqrt{\tfrac{8g^{2}}{\kappa% _{+}}}\tfrac{G\delta\Omega}{2G^{2}+\delta\Omega^{2}}}}\right]^{1}\cdot% \mathopen{}\mathclose{{}\left(\frac{1}{t}\int_{0}^{t}d\tau\hat{X}_{\rm out}^{% \pi/2}(\tau)}\right),\end{split}  (37) 
This estimator has been constructed so that at long times
\langle{\hat{n}_{b}^{\rm meas}(t)}\rangle\to\mathopen{}\mathclose{{}\left\{% \begin{array}[]{lc}n_{b,{\rm init}}&n_{b,{\rm init}}=0,1\\ f[n_{b,{\rm init}}]\,n_{b,{\rm init}}&n_{b,{\rm init}}\geq 2\end{array}}\right.  (38) 
where n_{b,{\rm init}} is the initial phonon number, and f[n] can be calculated from the effective temperature above, and goes to \mathopen{}\mathclose{{}\left[{\tfrac{2G\delta\Omega}{2G^{2}+\delta\Omega^{2}}% }}\right]^{1} for n\to\infty.
The behavior of the measurement output, expressed in terms of the phonon number estimator \hat{n}_{b}^{\rm meas}(t), are presented in Fig. 4, again for the ideal case of \gamma=\kappa_{}=0. As expected from our analysis, different phonon numbers are clearly resolvable after an integration time on the order of \tau_{\rm meas}. The details of these numerical calculations are given in Appendix C.
Note that in this ideal situation where there is no spurious dissipation, the resolving power of our measurement continues to increase indefinitely as the integration time t is increased. This reflects the shelving nature of the dynamics. We now go on to address the additional limitations placed by unwanted dissipation.
V Spurious dissipation: mechanical and auxiliary mode losses
Our analysis so far has treated the ideal case, where the only dissipation is the necessary coupling between the primary photonic mode and the waveguide (or transmission line) used to collect the output field which serves as the measurement record. This idealized analysis gives a useful picture for understanding the more realistic experimental case where there is also mechanical dissipation (damping rate \gamma) and damping of the auxiliary mode (damping rate \kappa_{}).
All three forms of spurious dissipation (mechanical loss, mechanical heating, auxiliary mode loss) cause the quantity \hat{N}_{\rm tot} to change, disrupting the shelving physics that is at the heart of our scheme. If spurious dissipation changes \hat{N}_{\rm tot}, the simple correspondence between the measurement record and the initial mechanical phonon number is lost. Focusing on small initial phonon numbers, the fastest rate at which such spurious processes occur is given by \max\mathopen{}\mathclose{{}\left[{\kappa_{},\gamma\mathopen{}\mathclose{{}% \left({2n_{\rm th}+1}}\right)}}\right]. Thus, one requires that the measurement essentially occur on a timescale much shorter than the inverse of this rate. This translates to the condition
\mathopen{}\mathclose{{}\left(\tau_{\rm meas}^{1}=\frac{g^{2}}{\kappa_{+}}}% \right)\gg\max\mathopen{}\mathclose{{}\left[{\kappa_{},\gamma\mathopen{}% \mathclose{{}\left({2n_{\rm th}+1}}\right)}}\right].  (39) 
In Fig. 5 we show the results of detailed master equation simulations including all forms of spurious dissipation which show this intuitive understand is correct.
It is interesting to consider two relevant limits of the condition in Eq. (39). In the case where \kappa_{} is much faster than the mechanical heating \gamma(2n_{\rm th}+1), the condition reduces to g^{2}\gg\kappa_{+}\kappa_{}. This is the same constraint that limits other proposals for measuring phonon number in twocavity optomechanical systems Miao2009 (); Ludwig2012 ().
given recent progress in constructing twomode optomechanical systems with extremely low auxiliarycavity dissipation Paraiso2015 (), it is worth considering the opposite limit, where \gamma\mathopen{}\mathclose{{}\left({2n_{\rm th}+1}}\right) dominates. In this case, our condition reduces to
\mathopen{}\mathclose{{}\left(C_{1}=\frac{4g^{2}}{\kappa_{+}\gamma}}\right)\gg 2% n_{\rm th}+1.  (40) 
where we have introduced the singlephoton cooperativity C_{1}. This regime is illustrated in Fig. 5. Numerically, we find the requirement is {C_{1}\gtrsim 100\mathopen{}\mathclose{{}\left({2n_{\rm th}+1}}\right)} to distinguish between the vacuum state and a single phonon, if the mechanics are the dominant source of spurious dissipation.
VI Outlook
We have presented and analyzed a new approach for QND phonon number measurement in quantum optomechanics, one which adapts the successful shelving strategy used in trappedion systems. By considering regimes where the driving of the photonic system is strong (and not perturbative), we are able to obtain measurement rates that are parametrically larger than previously considered approaches. At the level of theory, our approach of mapping a bosonic problem onto a spin system could be useful in more complex optomechanical systems where there are again conserved quantities akin to the excitation number \hat{N}_{\rm tot} in our system.
Our approach is most suited to systems where the auxiliary (unmeasured) photonic mode has an extremely low dissipation rate; recent experiments on twocavity microwavefrequency optomechanical systems Toth2016 () present a promising route to this regime.
Acknowledgements
This research was undertaken in part thanks to funding from NSERC and the Canada Research Chairs program.
Appendix A Analogy with electron shelving
The shelving protocol, pioneered by Hans Dehmelt Dehmelt1968 (); Nagourney1986 (); Javanainen1986 (), is a form of QND measurement of the state of twolevel atomic system. Given an ion in its ground state (\mathopen{}\mathclose{{}\left\lvert g}\right\rangle) or a metastable excited state (\mathopen{}\mathclose{{}\left\lvert e}\right\rangle), the protocol makes use of a shortlived auxiliary state (\mathopen{}\mathclose{{}\left\lvert a}\right\rangle). By applying a resonant laser drive we induce Rabi oscillations between \mathopen{}\mathclose{{}\left\lvert g}\right\rangle and \mathopen{}\mathclose{{}\left\lvert a}\right\rangle; this is accompanied by the spontaneous decay process \mathopen{}\mathclose{{}\left\lvert a}\right\rangle\to\mathopen{}\mathclose{{}% \left\lvert g}\right\rangle. The detection of photoemission in the appropriate wavelength is then confirmation that the atom is not in the excited state \mathopen{}\mathclose{{}\left\lvert e}\right\rangle. However, throughout the process, the atom cannot be said to be in the original ground state \mathopen{}\mathclose{{}\left\lvert g}\right\rangle. Instead it is in some hybrid state of \mathopen{}\mathclose{{}\left\lvert g}\right\rangle and \mathopen{}\mathclose{{}\left\lvert a}\right\rangle. This process is sketched out in figure Fig. (a)a.
Notice that the rate of spontaneous emission is low, and distinguishing the signal from the noise requires a prolonged measurement similar to what is discussed in Section IV. Thus the original shelving protocol also constitutes a weak measurement.
The analogy to our measurement scheme is most direct in the case where we are trying to distinguish an initial mechanical phonon number of \hat{n}_{b}=0 or \hat{n}_{b}=1. Here, the system is originally in the state \mathopen{}\mathclose{{}\left\lvert 1,0}\right\rangle or \mathopen{}\mathclose{{}\left\lvert 0,0}\right\rangle (where \mathopen{}\mathclose{{}\left\lvert n_{b},\hat{n}_{}}\right\rangle denotes a Fock state with the respective number of phonons and auxiliary photons). The application of the drive G generates Rabi oscillations between \mathopen{}\mathclose{{}\left\lvert 1,0}\right\rangle and \mathopen{}\mathclose{{}\left\lvert 0,1}\right\rangle, while a decay process allows the relaxation \mathopen{}\mathclose{{}\left\lvert 0,1}\right\rangle\to\mathopen{}\mathclose{% {}\left\lvert 1,0}\right\rangle with emission of a primary photon, first into the cavity and then out of it through the dissipation channel \kappa_{+}. By detecting these emitted photons, we can determine the initial phonon number. This process is sketched out in figure Fig. (b)b.
Our scheme takes the basic shelving concept further, allowing one to distinguish between multiple possible values of \hat{n}_{b} while remaining nondestructive. As discussed in the text, information on the phonon number is encoded in the rate at which emitted photons are produced.
Appendix B Optimized measurement prototocol
As discussed in Section IV.2, phonon number measurement in the steady state is subject to two conflicting constraints: the magnitude of \langle{\boldsymbol{\hat{J}}}\rangle is proportional to B_{z}/B while the measurable cavity field is coupled to the transverse components, \hat{c}_{+}\sim\langle{\hat{J}_{+}}\rangle\sim B_{x}/B\mathopen{}\mathclose{{}% \left\lvert\langle{\boldsymbol{\hat{J}}}\rangle}\right\rvert.
These competing requirements on the effective field can be avoided by arranging parameters so that the system has a long relaxation time, and performing the measurement in the transient period before a steady state is reached. We can then align the spin entirely along the axis of measurement and increase the magnitude of the cavity field. As we have seen in Eq. 24, a long relaxation time is achieved in the regime B\gg\kappa_{+}. We therefore suggest the following protocol:

The system is set up with some large energy mismatch \delta\Omega\gtrsim\kappa_{+}. This ensures a slow relaxation to the steady state once the primary mode drive is turned on, the relevant rates being {\Gamma_{\pm}\sim\mathopen{}\mathclose{{}\left({\kappa_{+}/B}}\right)^{2}\tau_% {m}^{1}}.

With the driving turned off, G=0, the mechanics are prepared in their initial state. As discussed in Section II, the auxiliary mode remains at zero population in the absence of driving. This corresponds to the spin pointing along \mathbf{e}_{z}.

We then ramp up the primary mode drive semiadiabatically, until G\gg\delta\Omega. In the spin language, we increase the B_{x} component until the effective field is nearly parallel with \mathbf{e}_{x}.

If the driving is increased sufficiently slowly, the spin itself will track the direction of the field and finally point along \mathbf{e}_{x}. We make use of the separation of scales to ramp over a time t_{\rm ramp} such that
B^{1}\ll\kappa_{+}^{1}\ll t_{\rm ramp}\ll\tau_{\rm meas}. (41) 
Finally, we perform a homodyne measurement as described in Section IV. Here we use the estimator
\hat{n}_{b}^{\rm meas}\mathopen{}\mathclose{{}\left({t}}\right)\equiv\mathopen% {}\mathclose{{}\left[{\sqrt{\tfrac{2g^{2}}{\kappa}}}}\right]^{1}\cdot% \mathopen{}\mathclose{{}\left(\frac{1}{t}\int_{0}^{t}d\tau\hat{X}_{\rm out}^{% \pi/2}(\tau)}\right). (42)
This procedure is illustrated in Fig. 7. We see that the required measurement time is about three times shorter for this set of parameters.
Appendix C Details of Numerical Simulations
The numerical calculations presented throughout this work are the results of full master equation treating each of the three bosonic modes in our system:
\begin{split}\displaystyle\dot{\hat{\rho}}=i\mathopen{}\mathclose{{}\left[{% \hat{\rho},\hat{H}_{\rm eff}}}\right]&\displaystyle +\kappa_{+}\mathcal{D}% \mathopen{}\mathclose{{}\left[{\hat{c}_{+}}}\right]\cdot\hat{\rho}+\kappa_{}% \mathcal{D}\mathopen{}\mathclose{{}\left[{\hat{c}_{}}}\right]\cdot\hat{\rho}% \\ &\displaystyle+\mathopen{}\mathclose{{}\left[{\gamma\mathopen{}\mathclose{{}% \left({n_{\rm th}+1}}\right)\mathcal{D}[\hat{b}]+\gamma n_{\rm th}\mathcal{D}[% \hat{b}^{\dagger}]}}\right]\cdot\hat{\rho},\end{split}  (43) 
where \hat{H}_{\rm eff} is given in Eq. 6 and
\mathcal{D}\mathopen{}\mathclose{{}\left[{\hat{a}}}\right]\cdot{\hat{\rho}}=% \hat{a}\hat{\rho}\hat{a}^{\dagger}\tfrac{1}{2}\mathopen{}\mathclose{{}\left({% \hat{a}^{\dagger}\hat{a}\hat{\rho}+\hat{\rho}\hat{a}^{\dagger}\hat{a}}}\right)  (44) 
is the standard Lindblad superoperator.
The numerical simulations were performed using QuTiP Johansson2013 (). We truncated the Hilbert space to {\hat{c}_{+}^{\dagger}\hat{c}_{+}=0,1,2}, and to 0\leq\hat{N}_{\rm tot}\leq 6 for simulations involving spurious dissipation (truncation is not needed if \kappa_{}=\gamma=0). We then solved the eigenvalue problem defined by Eq. 43 and calculated \hat{\rho}\mathopen{}\mathclose{{}\left({t}}\right) at any time as a sum of exponentials.
The initial state in all cases was \hat{\rho}\mathopen{}\mathclose{{}\left({0}}\right)=\mathopen{}\mathclose{{}% \left\lvert\psi_{0}}\right\rangle\mathopen{}\mathclose{{}\left\langle\psi_{0}}\right\rvert where
\hat{b}^{\dagger}\hat{b}\mathopen{}\mathclose{{}\left\lvert\psi_{0}}\right% \rangle=n_{b}\mathopen{}\mathclose{{}\left\lvert\psi_{0}}\right\rangle\qquad% \hat{c}_{\pm}\mathopen{}\mathclose{{}\left\lvert\psi_{0}}\right\rangle=0.  (45) 
To simulate the driving of the primary mode, the parameter G is initially ramped from zero to its final value. For Figs. 5 and 4 we ramp G linearly over a time span of \kappa_{+}^{1}, while in Fig. 7 we increase it exponentially as explained in the figure caption. Experimentally these behaviors can be achieved using pulseshaping techniques.
The standard deviations plotted are given by
\Delta n_{b}^{\rm meas}\mathopen{}\mathclose{{}\left({t}}\right)=\sqrt{\langle% {\mathopen{}\mathclose{{}\left({\hat{n}_{b}^{\rm meas}\mathopen{}\mathclose{{}% \left({t}}\right)}}\right)^{2}}\rangle\langle{\hat{n}_{b}^{\rm meas}\mathopen% {}\mathclose{{}\left({t}}\right)}\rangle^{2}}.  (46) 
They were calculated using the quantum regression theorem Gardiner2004 (), taking into account the full dynamics of the system and not just cavity vacuum noise.
References
 (1) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt. Cavity Optomechanics. Springer, Berlin (2014).
 (2) J. D. Teufel et al. Nature 475 359, (2011).
 (3) J. Chan et al. Nature 478 89, (2011).
 (4) R. W. Peterson et al. Phys. Rev. Lett. 116 063601, (2016).
 (5) D. W. C. Brooks et al. Nature 488 476, (2012).
 (6) A. H. SafaviNaeini et al. Nature 500 185, (2013).
 (7) T. P. Purdy, P.L. Yu, R. W. Peterson, N. S. Kampel, and C. A. Regal. Phys. Rev. X 3 031012, (2013).
 (8) E. E. Wollman et al. Science 349, (2015).
 (9) J.M. Pirkkalainen, E. Damskägg, M. Brandt, F. Massel, and M. Sillanpää. Phys. Rev. Lett. 115 243601, (2015).
 (10) F. Lecocq, J. B. Clark, R. W. Simmonds, J. Aumentado, and J. D. Teufel. Phys. Rev. X 5 041037, (2015).
 (11) D. H. Santamore, A. C. Doherty, and M. C. Cross. Phys. Rev. B 70 144301, (2004).
 (12) I. Martin and W. H. Zurek. Phys. Rev. Lett. 98 120401, (2007).
 (13) J. D. Thompson et al. Nature 452 72, (2008).
 (14) A. M. Jayich et al. New J. Phys. 10 095008, (2008).
 (15) H. Miao, S. Danilishin, T. Corbitt, and Y. Chen. Phys. Rev. Lett. 103 100402, (2009).
 (16) A. A. Gangat, T. M. Stace, and G. J. Milburn. New J. Phys. 13 043024, (2011).
 (17) Y. Yanay, J. C. Sankey, and A. A. Clerk. Phys. Rev. A 93 063809, (2016).
 (18) M. Ludwig, A. H. SafaviNaeini, O. Painter, and F. Marquardt. Phys. Rev. Lett. 109 063601, (2012).
 (19) P. Kómár et al. Phys. Rev. A 87 013839, (2013).
 (20) S. BasiriEsfahani, U. Akram, and G. J. Milburn. New J. Phys. 14 085017, (2012).
 (21) T. K. Paraïso et al. Phys. Rev. X 5 041024, (2015).
 (22) L. D. Tóth, N. R. Bernier, A. Nunnenkamp, A. K. Feofanov, and T. J. Kippenberg. arXiv:1602.05180, (2016).
 (23) H. G. Dehmelt. Adv. At. Mol. Phys. 3 53, (1968).
 (24) W. Nagourney, J. Sandberg, and H. Dehmelt. Phys. Rev. Lett. 56 2797, (1986).
 (25) J. Javanainen. Phys. Rev. A 33 2121, (1986).
 (26) M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt. Rev. Mod. Phys. 86 1391, (2014).
 (27) H. K. Cheung and C. K. Law. Phys. Rev. A 84 023812, (2011).
 (28) A. Xuereb and P. Domokos. New J. Phys. 14 095027, (2012).
 (29) C. W. Gardiner and P. Zoller. Quantum Noise. Springer, Berlin/Heidelberg (2004).
 (30) A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf. Rev. Mod. Phys. 82 1155, (2010).
 (31) D. F. Walls and G. J. Milburn. Quantum optics. Springer, Berlin/Heidelberg, 2nd edition (1997).
 (32) J. M. Manley and H. E. Rowe. Proc. Inst. Radio Eng. 44 904, (1956).
 (33) M. I. Dykman, M. Marthaler, and V. Peano. Phys. Rev. A 83 052115, (2011).
 (34) V. Peano and M. Thorwart. Phys. Rev. B 82 155129, (2010).
 (35) M.A. Lemonde and A. A. Clerk. Phys. Rev. A 91 033836, (2015).
 (36) J. R. Johansson, P. D. Nation, and F. Nori. Comput. Phys. Commun. 184 1234, (2013).