Shelving-style QND phonon-number detection in quantum optomechanics

We propose a new method for optomechanical quantum non-demolition (QND) detection of phonon number, based on a"shelving"style measurement. The scheme uses a two-mode optomechanical system where the frequency splitting of the two photonic modes is near-resonant with the mechanical frequency. The combination of a strong optical drive and the underlying nonlinear optomechanical interaction gives rise to spin-like 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 [1]. They have enabled the exploration of a wide range of effects, including the cooling of macroscopical mechanical modes to near their quantum ground states [2][3][4] and the generation of squeezed states of light [5][6][7] and mechanical motion [8][9][10].
A key unrealized goal in optomechanics is the quantum non-demolition (QND) measurement of mechanical phonon number, something that would directly reveal the energy quantization of the mechanics [11,12]. 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 phonon-number detection both in the regime where the mechanical frequency is much smaller than the splitting between the optical modes [13][14][15][16][17], and in the regime where it is comparable to this splitting [18][19][20]. Two-cavity optomechanical systems have been realized experimentally using a variety of platforms [13,21,22].
In this work, we revisit the two-mode system of Fig. 1a, 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, many-photon 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. [18], who considered more non-resonant interactions and weak driving, such that a pertubative treatment was possible. It also contrasts from the work of Kómár et al. [19], 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. 1b, and is conceptually similar to the highly successful electron shelving technique used in trapped-ion systems [23][24][25],  1. (a) Schematic of the driven two-cavity optomechanical system described by Eqs. (1) to (4). Two photonic modes (ĉ±) are coupled nonlinearly to a mechanical modeb (amplitude g). In addition, a strong drive on the primaryĉ+ mode realizes a beam-splitter interaction betweenb andĉ− (amplitude G). κ± denote the photonic mode damping rates, γ the mechanical damping rate. (b) Schematic of measurement scheme. G causes excitations to oscillate between the mechanics and the auxiliaryĉ− mode. The nonlinear interaction converts aĉ− photon into a phonon and aĉ+ photon; this photon is then detected. The scheme allows a QND measurement ofNtot =ĉ † −ĉ − +b †b , and thus QND measurement of the initial phonon number.
as explained in Appendix A. The measurement starts by turning on the strong optical drive on the "pri-mary" mode. This induces coherent oscillations between phonons and photons in the undriven "auxiliary" mode, oscillations which necessarily conserve the total number N tot of aux-mode and mechanical excitations. The remaining nonlinear interaction then allows one to use the output field of the driven optical mode to measureN tot without changing its value -a true QND measurement. This is then effectively a measurement of mechanical phonon number, as without spurious dissipative effects, N tot is simply equal to the initial phonon number.
As with other schemes, our approach still requires a relatively strong single-photon 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. 18, or in the adiabatic regime of a position-squared measurement considered in Ref. 13-15. This faster rate could be of significant utility in new experimental designs of two-cavity optomechanical systems where one photonic mode is explicitly engineered to have extremely low dissipation [21,22]. 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 single-photon cooperativity [26], 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 Hamil-tonianĤ The first term describes the two photonc eigenmodes in the absence of optomechanical coupling, and is ( = 1 throughout) whereâ ± are annihilation operators for the two photonic modes (frequencies ω c − J and ω c + J respectively). We call these the primary (â + ) and auxiliary (â − ) modes.
The second term in Eq. (1) describes a mechanical mode with frequency ω m and annihilation operatorb: We have introduced the parameter δΩ = 2J − ω m which represents the mismatch between the photonic mode splitting and the mechanical frequency. Similar to reference [18], we will be interested in the regime where the mechanical frequency is similar in magnitude to 2J, and hence |δΩ| J. Finally, the last term in Eq. (1) describes the optomechanical coupling, which takes the form of a mechanically-mediated tunneling interaction between the two modes:Ĥ where g is the single-photon optomechanical coupling amplitude. For the specific case of a membrane-inthe-middle style optomechanical system [13], a detailed derivation of this model is given in Refs. 14 and 27; its limitations are discussed in Refs. 27 and 28. We next use a standard treatment to describe the effects of dissipation on the system [29,30]. We assume that each optical mode is coupled to an independent zerotemperature Markovian reservoir, giving rise to an amplitude damping rate of κ ± /2 forâ ± , respectively. The mechanics are similarly coupled to a thermal Markovian reservoir characterized by a thermal occupancy n th , leading to an amplitude damping rate of γ/2 forb. Note that for a membrane-in-the-middle type system, treating the optical modes as seeing independent reservoirs can be problematic [17], 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,â + , mode near its resonance frequency ω 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: Here, â + g=0 is average primary-mode amplitude induced by the drive in the absence of optomechanical coupling, andĉ + is the displaced annihilation operator for this mode in the rotating frame.
In this frame we findĤ S =Ĥ eff +Ĥ NR wherê Finally, we assume that the optical mode gap 2J is sufficiently large that we can safely make a rotating wave approximation and drop the non-resonant interactions described byĤ NR . This requires in practice J |δΩ|, G, g, κ + . 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
A. Conserved excitation number The nonlinear interaction in Eq. (6) has the form of a standard three-wave mixing Hamiltonian, describing e.g. a non-degenerate parametric amplifier with a dynamical pump mode. It has three conserved excitation numbers (typically known as Manley-Rowe constants of motion [32]) reflecting the SU(1) and SU(1,1) symmetries of the interaction. Including the beam-splitter interaction reduces the symmetry of the Hamiltonian, andĤ eff has only a single conserved excitation numberN tot : Note that asN 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 input-output port used for the measurement. This corresponds to κ + = 0, while γ = κ − 0. Starting with this limit will let us present a clear picture of how our proposed phonon number measurement operates: the conservation ofN 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 κ + γ, κ − is also relevant experimentally. While mechanical damping rates γ are typically orders of magnitude smaller than photonic damping rates, it might first seem surprising that one could achieve κ − κ + . However, recent experiments with microwave-circuit optomechanics have achieved precisely this kind of regime through clever experimental design.
We analyze the effects of non-zero γ, κ − in Section IV. As expected, they limit the maximum possible measurement time, putting constraints on whether phonon number detection is achievable.

B. Mapping to an effective spin system
To understand the dynamics of our system, it is helpful to make the conservation ofN tot explicit by representing the mechanical and aux-modes by an effective spin (a Schwinger boson representation in reverse). We thus introduce the effective (dimensionless) angular momentum operatorsĴ viaĴ These operators satisfy standard angular momentum commutation relations. A direct calculation shows: The conservation ofN tot now just corresponds to our effective spin having a fixed length. With these operators, the effective Hamiltonian iŝ where the effective magnetic field is with e x , e y , e z the unit vectors in the respective directions. As theĉ + mode is damped, we see thatĤ eff takes the form of a generalized spin-boson model. The quadratic terms in the original version ofĤ eff now correspond to a Zeeman field on our spin. The order g nonlinear term inĤ eff now corresponds to a system-bath coupling where a bath boson can be created while adding angular momentum in the z direction to the spin. The Heisenberg-Langevin equations of motion take the formċ where e ± = e x ± ie y . The operatorĉ in describes the vacuum fluctuations entering the system from the inputoutput line coupled to the principal + mode (i.e. from the input-output waveguide that will be used for measurement); it has zero mean and correlation functions When the optomechanical drive is initially turned on, N tot =n b , wheren b ≡b †b is the initial mechanical The z-projection of the effective spin represents the difference between the number of phonons and number of aux-mode quanta.
The many-photon optomechanical interaction G exchanges photons and phonons, and hence acts like a magnetic field in ex direction.

Bz δΩ
The mismatch between the photonic mode splitting and the mechanical frequency acts as an ez field.
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Ĵ . The nonlinear interaction creates photons in the principal modeĉ + in a way that depends onĴ ± , and hence on the magnitude of 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 I.

C. Dynamics in the limit g κ+
We consider the system dynamics in the typical to experiment, κ + g. We begin by solving Eq. (13) to find . We see thatĉ + (t) depends on the behaviour ofĴ at earlier times, and that the range of this memory effect is κ −1 + . We can solve for the spin dynamics during this short interval to a good approximation by ignoring the effects of g, aŝ As long as τ is at most ∼ κ + , we can drop the O(gτ ) term as it is proportional to the small parameter With this approximation, the solution forĉ + (t) takes the simpler form: Here, the non-Hermitian operatorĴ c is linear in the components ofĴ , and given bŷ The initial transient behaviour ofĉ + is described bŷ Equation (19) implies that the primary mode (and its output optical field) will be sensitive to components of J in a manner that depends on both B and κ + . In the extreme Markovian limit κ + B we findĴ c →Ĵ + , implying that theĉ + is able to measure both transverse components of our effective spin. In contrast, in the limit B κ + , theĉ + mode is only able to respond to the non-rotating component of the spin, and henceĴ c → 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ödinger-picture density matrixρ J = Trĉ We again consider the g κ + limit, and derive a weakcoupling master equation, retaining terms to leading order in η. For simplicity, we also consider the typical limit where B ηg = 2g 2 /κ + , allowing us to drop rapidly oscillating terms (i.e. make a secular approximation).
We useĴ = e B ·Ĵ to denote the projection of the spin along the axis of the effective magnetic field B , and useĴ ± to denote raising and lowering operators forĴ . Defining D[â]·ρ =âρâ † − 1 2 â †âρ +ρâ †â as the standard Lindblad superoperator, we finḋ 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 bŷ The remaining terms describe dephasing (rate Γ ϕ ), and spin raising / lowering transitions along the axis e B (rates Γ ± ). The dissipative rates are Note that unless B z = B, the vacuum noise driving theâ + mode is able to cause spin flips upwards in energy: this is possible because we are strongly driving thê a + mode, and can be viewed as being analogous to the phenomenon of "quantum heating" [33][34][35]. Note also that we always assume that the drive frequency ω c − J is much larger than G, |δΩ|. From Eq. (22) we find that the steady state of the spin is simply a thermal state, with an effective temperature 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 δΩ 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).

D. QND nature of the measurement
AsN tot is a conserved quantity, it is clear that our protocol allows for a QND measurement ofN tot : this quantity commutes with the system Hamiltonian, even when we include the coupling κ + to the input-output 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 ofN tot corresponds to the initial mechanical phonon numbern 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 untiln b =N tot , i.e. until all the original phonons are returned to the mechanical resonator.

A. 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ĉ + . Eq. (19) tells us that the output field will be sensitive to particular components of the spin vectorĴ . Recall that in our original lab frame, the primary mode is driven at its resonance frequency ω + = ω 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 γ = κ − = 0.
From standard input-output theory (and working in the interaction picture defined by Eq. (5)), the cavity output field leaving the + mode is given bŷ and the phonon number can be obtained from measuring an output quadrature, where the angle α parameterizes the choice of quadrature.
Using the result in Eq. (19), we find that the average value of this output quadrature is given by: In the steady state, our spin relaxes to a thermal state, and in general the value of Ĵ c (t) will be non-zero 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. (20) and (25) that in the steadystate Ĵ c (t) is purely real, implying that the optimal choice of quadrature corresponds to an angle α = π/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 toX α out ) 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 [14,30] indicate that the mechanical Fock staten b = n+1 can be distinguished from Fock staten b = n when the measurement time is sufficiently long, namely where and the prefactor is determined by As one can see from Fig. 3, for a broad range of G/δΩ, the prefactor 1/j 2 n is of order unity, and the timescale of the measurement is set by τ 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. [18]. Their scheme requires a measurement time which scales as t ∼ (δΩ/G) 2 τ meas to resolve two adjacent Fock states. However, the setup requires an interaction that is sufficiently non-resonant 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 |δΩ| and G |δΩ|, 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 ∼ |δΩ|. 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κ + to resolve phonon number states after t = 50κ −1 + (see Ref. 18 Fig. 3), Eq. (31) implies our shelving protocol could do the same in a system where g = 0.3κ + .  29), it is proportional to the amplitude of the measurement output. G is the many-photon optomechanical coupling, and δΩ is the mismatch between the mechanical frequency and the photonic mode splitting. In the spin mapping, 2G/δΩ = Bx/Bz. Different curves correspond to different initial phonon numbers n b . The steady-state value of Ĵ c is sensitive to phonon number for a wide range of G/δΩ. All curves correspond to the case where there is no spurious dissipation, γ = κ− = 0, and the limit 4G 2 + (δΩ) 2 g 2 /κ+. Note that while Ĵ c(t) is generally complex, it is real in the steady state.

B. Optimal parameters
Our system has a fairly large number of tunable parameters. In particular, both non-zero components of the effective field B = (G, 0, δΩ) can be selected for optimal performance. We remind the reader that G is the manyphoton coupling strength and δΩ the mismatch between the mechanical frequency and photonic mode splitting.
A key parameter to optimize the size of the signal in the long-time homodyne current which encodes the initial phonon number. From Eq. (29), this is directly determined by the steady-state value of Ĵ c . In steady state, we obtain from Eqs. (20) and (25), where Ĵ T eff is the thermal average ofĴ = e B ·Ĵ taken at a temperature T 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 Ĵ c we need B x and B z to be comparable. A finite value of B x (i.e. many-photon 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Ĵ + . A non-zero value of B z (i.e. mismatch between mechanical frequency and photonic mode splitting) is needed so that the effective temperature describing the steady-state 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, (34) Figure 3 shows the steady state value of J c for other values of initial phonon number. We find the optimal ratio, for a system with one or a few initial phonons, is While the steady-state homodyne signal is not sensitive to the size of 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 κ + , then it follows from Eqs. (22) and (24) 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 τ meas , requiring B κ + . That is, one should avoid the many-photon 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 ∼ B) g 2 /κ + . This condition is easy to meet in experiment given the typical weakness of single-photon optomechanical coupling strengths. Putting these conditions together, we find that the optimal regime for phonon number measurement requires: An alternative approach to the measurement which results in faster measurement rates but is more complicated to implement experimentally is discussed in Appendix B.

C. 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ĉ + . We construct an estimator for phonon number from the integrated homodyne current: This estimator has been constructed so that at long times where n b,init is the initial phonon number, and f [n] can be calculated from the effective temperature above, and goes to 2GδΩ 2G 2 +δΩ 2 −1 for n → ∞. (t) as a function of integration time t. The estimator, obtained from the time-integrated homodyne current, is described in Eq. (37), and time is given in units of τmeas = κ+/g 2 . Curves correspond to the mean value of the estimator, shaded regions correspond to the standard deviation. Different curves correspond to different initial mechanical phonon numbers. One sees that different initial mechanical phonon numbers become distinguishable on a time scale proportional to τmeas. Parameters correspond to g = 0.01κ+, G = 0.1κ+, δΩ = 0.13κ+, and to the ideal case where γ = κ− = 0.
The behavior of the measurement output, expressed in terms of the phonon number estimatorn meas b (t), are presented in Fig. 4, again for the ideal case of γ = κ − = 0. As expected from our analysis, different phonon numbers are clearly resolvable after an integration time on the order of τ 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 γ) and damping of the auxiliary mode (damping rate κ − ).
All three forms of spurious dissipation (mechanical loss, mechanical heating, auxiliary mode loss) cause the quantityN tot to change, disrupting the shelving physics that is at the heart of our scheme. If spurious dissipation changesN tot , the simple correspondence between the measurement record and the initial mechanical phonon (t) as a function of integration time t when spurious dissipation is taken into account. As in Fig. 4, the estimator is described in Eq. (37), and time is in units of τmeas = κ+/g 2 . Curves correspond to the mean value of the estimator and shaded regions correspond to the standard deviation. Insets show the behavior of N tot over the same range of times. We consider parameters where mechanical dissipation is the key limit to the scheme; as such, the relevant parameter is the single-photon cooperativity C1 = 4g 2 /κ+γ. When spurious dissipation is negligible (left), N tot is nearly conserved and different initial mechanical phonon numbers n b are easily resolved if we integrate on a time scale longer than τmeas. When spurious dissipation is relevant but sufficiently small, the integrated homodyne current can still be used to determine the initial phonon number n b at intermediate times. At longer times, the mechanical system begins to relax into a thermal equilibrium independent of its initial state. If the rate of spurious dissipation is too fast (right) the system relaxes before differences in n b are resolvable. Parameters used are g = 0.1κ+, G = 0.1κ+, δΩ = 0.13κ+, κ− = 10 −4 κ+, n th = 100 and γ varies. The details of the numerical integration are given in Appendix C.
number is lost. Focusing on small initial phonon numbers, the fastest rate at which such spurious processes occur is given by max[κ − , γ(2n th + 1)]. Thus, one requires that the measurement essentially occur on a timescale much shorter than the inverse of this rate. This translates to the condition 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 κ − is much faster than the mechanical heating γ(2n th + 1), the condition reduces to g 2 κ + κ − . This is the same constraint that limits other proposals for measuring phonon number in two-cavity optomechanical systems [15,18].
given recent progress in constructing two-mode optomechanical systems with extremely low auxiliarycavity dissipation [21], it is worth considering the opposite limit, where γ(2n th + 1) dominates. In this case, our condition reduces to where we have introduced the single-photon cooperativity C 1 . This regime is illustrated in Fig. 5. Numerically, we find the requirement is C 1 100(2n th + 1) 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 trapped-ion 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 numberN 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 microwave-frequency optomechanical systems [22] present a promising route to this regime. (a) The measurement protocol. At time t = 0 the system is initiated in a pure state. From time 0 < t < (t f = 0.1τmeas), the driving is ramped up, with G(t) = G f /100 1−t/t f . We then keep G fixed and the measurement protocol outlined in Section IV is performed. (b-c) Behavior of the spin system. Initially, with δΩ G, the spin is aligned along ez. When the driving is ramped up semi-adiabatically, the direction of the spin rotates with it to point along ez. At a later time, t ∼ (Γ+ + Γ−) −1 (see Eq. (24)), the magnitude of the spin decays into its steady state value. (d) The phonon number estimator as a function of time during the measurement. This plot corresponds to the gray shaded area on the left. Curves correspond to the mean value of the estimator, and shaded regions correspond to the standard deviation. The different phonon numbers can be differentiated about three times faster than in Fig. 4. The time is in units of τmeas = κ+/g 2 . The parameters used here are g = 0.01κ+, δΩ = κ+, G f = 5κ+ and with no spurious dissipation, κ− = γ = 0. Details about the numerical calculation are given in Appendix C.

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: whereĤ eff is given in Eq. (6) and is the standard Lindblad superoperator. The numerical simulations were performed using QuTiP [36].
To simulate the driving of the primary mode, the parameter G is initially ramped from zero to its final value. For Figs. 4 and 5 we ramp G linearly over a time span of κ −1 + , while in Fig. 7 we increase it exponentially as explained in the figure caption. Experimentally these behaviors can be achieved using pulse-shaping techniques.
The standard deviations plotted are given by They were calculated using the quantum regression theorem [29], taking into account the full dynamics of the system and not just cavity vacuum noise.