text:
# A Quantum Non-demolition measurement of Fock states of mesoscopic mechanical oscillators

## I Introduction

Quantum mechanics tells us that the energy of an oscillator is quantized. However, an observation of quantum limited mechanical motion in macroscopic objects has not been possible because the energy associated with individual phonons is typically much smaller than the thermal energy C03 ; R01 .

Advances in nanotechnology have enabled experimenters to build ever smaller mechanical oscillators with high resonance frequencies and quality factors HZMR03 . As an individual phonon energy becomes comparable to or greater than $`k_\mathrm{B}T`$, quantum effects begin to appear and it should be possible to realize various quantum phenomena.

In this paper, we investigate the possibility of observing transitions amongst the Fock states of a mesoscopic mechanical oscillator. To do this requires the coupling of the system oscillator to a measurement device that sensitively detects the phonon number of the system oscillator but does not itself change the excitation level of the oscillator. In the quantum regime it becomes very important to model the precise way that a quantum system interacts with any measuring apparatus, as well as with the environment. Specifically, it is necessary to take into account the measurement backaction and to design the system-readout interaction so as to allow the best possible measurement of the desired observable. We will show that it is possible in principle to take advantage of the non-linear interaction between modes of oscillation of an elastic beam or beams to track the state of the oscillator as it jumps between number states due to its coupling to the surrounding thermal environment.

The laws of quantum mechanics tell us that, even in the absence of instrumental or thermal noise, a measurement will tend to disturb the state of the measured system. The interaction between the system and the measurement apparatus means that while information about the measured observable may be read out from the state of the meter after interacting with the system, the quantum mechanical uncertainty in the initial state of the meter leads to random changes in the conjugate observable of the system. This backaction noise on the conjugate observable is an inevitable result of the very interaction that allows the measurement to take place. It has long been recognized that such backaction noise places a fundamental limit on the sensitivity of physical measurements BK . However, the class of measurements known as quantum non-demolition (QND) measurements partially circumvents this problem by guaranteeing that the backaction noise does not affect the results of future measurements of the same quantity. The idea of a QND measurement is widely discussed in the literature (for example, see BK , CTDSZ80 , and WM ). In a QND measurement the interaction Hamiltonian between system and meter commutes with the internal Hamiltonian of the system: an ideal QND measurement is *repeatable* since the backaction noise does not affect the dynamics of the measured observable. In this paper, we are interested in a QND measurement of phonon number. The conjugate observable of number is phase, thus, the measurement backaction in our case will result in diffusion of the phase of the mechanical oscillations. However, the scheme allows the complete determination of the oscillator excitation level and thus projects the system onto a number state in an idealized limit.

The scheme for the QND measurement of phonon number that we consider uses two anharmonically coupled modes of oscillation of a mesoscopic elastic structure. The resonant frequencies of these two modes are different. The higher frequency mode is the system to be measured, while the lower frequency oscillator serves as the meter (we refer this oscillator as ancilla). The key idea of the scheme is that, from the point of view of the readout oscillator, the interaction with the system constitutes a shift in resonance frequency that is proportional to the time-averaged excitation of the system oscillator. This frequency shift may be detected as a change in the phase of the ancilla oscillations when driven on resonance. We show that this scheme realizes an ideal QND measurement of phonon number in an appropriate limit. To measure the phase of the ancilla oscillator we imagine a magnetomotive detection scheme so that the actual physical quantity measured is an electric current that couples to the ancilla displacement. Thus our task is to understand how the strong measurement of the current, represented by the von Neumann projection scheme on the current operator, yields information on the system phonon number, and in turn affects the dynamics of the system via the indirect coupling through the ancilla oscillator, and in the presence of the inevitable coupling of the ancilla and system oscillators to the environment. This QND measurement scheme where a nonlinear potential provides a phase shift to one oscillator that reflects the excitation of the other is analogous to the experiment of Peil and Gabrielse PG99 , which demonstrated a QND measurement of the excitation of a single trapped electron. Theoretical discussions of such approximate QND measurements of the excitation of an oscillator date back at least to Unruh U78 .

In Sec. II we introduce our model and construct the Hamiltonian describing the two oscillators, the magnetomotive coupling, and the coupling to the environment represented by baths of harmonic oscillators. For the ancilla oscillator displacement to directly indicate the system phonon number the time scale of the ancilla dynamics must be much shorter than that of the system. This actually allows us to adiabatically eliminate the ancilla operators to obtain dynamical equations for the system alone. Thus, in Sec. III, we obtain a reduced master equation for the density matrix of the system oscillator, which allows us to focus on the physics of the system dynamics. This adiabatic elimination of ancilla degrees of freedom is often considered in quantum optics. However, the adiabatic elimination used in quantum optics is at temperature zero, and we need to reformulate the method for finite temperatures.

Once we know the system dynamics, we next focus on obtaining the experimental outcome. Quantum mechanics allows us to determine the state of the system *conditioned* on the measured current $`I(t)`$. The von Neumann projection postulate says that after a measurement a quantum system in some possibly mixed initial state is projected onto the eigenstate corresponding to the measurement outcome. The continual measurement and projection of the current $`I(t)`$ provides accumulating information on the system phonon number, and correspondingly a projection onto phonon number states. The theory of quantum trajectories B86 ; B87 ; C93 ; Wth  has been developed to deal with such continuous measurements. In Secs. IV we discuss such quantum trajectory equations for our system. The method leads to a stochastic master equation for the system density matrix, where the stochastic component comes from the particular value of the measured current at each time, which itself is a stochastic variable since it is the outcome of a quantum measurement. These equations of motion for the system conditioned on a particular sequence of measurement results allow us to investigate the possibility of tracking the evolution of the system as it jumps between number states due to its interaction with the thermal bath. Some details of the formulation of the operators describing the measurement current that are needed to derive the stochastic master equation are given in Appendix A.

The main discussion of the physical implications of the model is in Sec. V where we consider the parameters that are necessary to observe the oscillator jump between number states. A reader not interested in the details of the derivation could read Sec. V and following sections after the description of the model in Sec. II. Finally, in Sec. VI we conclude with a discussion of the feasibility of the scheme based on current technology and future enhancements.

## II Constructing the Hamiltonian

### II.1 The model

In this section we introduce the model system and show how the coupling between the system and ancilla approximates a QND coupling in an appropriate limit. We then derive equations of motion that take into account the couplings to the environment and the interactions that drive and monitor the oscillations of the ancilla.

Consider a mesoscopic beam with rectangular cross section. There are two orthogonal flexing modes that are not coupled in the linear elasticity theory, but are coupled anharmonically. This coupling exists in nature between the two orthogonal flexing modes of single mechanical beam. However, the coupling can also be controlled and engineered: a similar coupling of bending modes in two elastic beams has been proposed by Yurke SP:Y . In this scheme two mesoscopic elastic beams with rectangular cross section are connected by a series of mechanical coupling devices. These devices have the effect of allowing only one type of strain (the longitudinal stretch) to pass to the other beam. In this paper we focus on the extent to which such mechanical devices are able to realize a QND measurement and the constraints this places on the specifications of the device, and the temperature at which the experiment is performed.

For convenience, we refer to the system of interest as oscillator $`0`$ and the ancilla as oscillator $`1`$, and the corresponding resonant frequencies of the two modes as $`\omega _0`$ and $`\omega _1`$, respectively. The ancilla is driven at its resonant frequency, and a measurement apparatus is attached to the ancilla. The whole structure is kept at a low temperature $`T`$ such that $`\mathrm{}\omega _0k_BT`$, where $`\mathrm{}=h/2\pi `$ is Plank’s constant, and $`k_B`$ is Boltzmann’s constant. The oscillators are weakly coupled to the environment. Figure 1 shows a schematic of our model.

### II.2 System Hamiltonian

Converting the schematic model in the last section into a tractable mathematical model and obtaining the system dynamics requires some assumptions and simplifications.

Firstly, we focus on the anharmonic coupling and the limit in which it satisfies the QND condition. In linear elasticity theory, the two flexing modes, which are perpendicular to each other, propagate independently without interacting. Beyond the linear approximation these modes are coupled. Expansion of the elastic energy with respect to the strain tensor is taken up to second order in the harmonic approximation. By symmetry the coupling between the modes first occurs at fourth order, proportional to $`x_0^2x_1^2`$. So we expand the anharmonic terms up to quartic order to give

$$H_{\mathrm{anh}}=\mathrm{}\left(\stackrel{~}{\lambda }_{00}x_0^4+\stackrel{~}{\lambda }_{11}x_1^4\right)+\mathrm{}\stackrel{~}{\lambda }_{01}x_0^2x_1^2,$$
(1)
where $`\stackrel{~}{\lambda }_{ij}`$ give the strengths of the nonlinear terms. The first two terms in Eq. (1) are internal anharmonic terms. Under the rotating wave approximation (see below), these terms cause a shift in the oscillator resonant frequencies and a non-linear phase shift that depends on intensity (a Kerr non-linearity) resulting in rotational shear of the state in the phase space of the two oscillators. For the system oscillator, the small shifts in the energy level spacings are not important, and can be ignored. The ancilla oscillator is externally driven, and so the nonlinearity of this oscillator may become large, for example leading to multistability for large enough drive strengths. We will assume that the drive strength is kept smaller than this range, so that again the $`x_1^4`$ nonlinearity does not play an essential role. However the $`x_0^2x_1^2`$ coupling plays an essential role in coupling the system and ancilla. Therefore, in the interests of a straightforward presentation we retain the non-linear coupling given by $`\lambda _{01}`$ and disregard the nonlinearities of the system and ancilla internal Hamiltonians. A detailed analysis including non-interacting nonlinearities and detuning in Ref. (qph04 ) has shown that in the regime of strong damping of the ancilla that we will mostly consider the effect of these anharmonic terms will be negligible for small detuning.

In terms of creation and annihilation operators, the Hamiltonian is now $`H=H_0+V_{\mathrm{anh}}`$, with $`H_0`$ the harmonic part

$$H_0=\mathrm{}\omega _0a_0^{}a_0+\mathrm{}\omega _1a_1^{}a_1,$$
(2)
and $`V_{\mathrm{anh}}`$ the anharmonic coupling

$$V_{\mathrm{anh}}=\frac{1}{4}\mathrm{}\lambda _{01}\left(a_0^{}+a_0\right)^2\left(a_1^{}+a_1\right)^2.$$
(3)
We have defined the standard raising and lowering operators for the oscillators

$$a_i=\sqrt{m_i\omega _i/2\mathrm{}}x_i+i\sqrt{1/2\mathrm{}m_i\omega _i}p_i,$$
(4)
and $`a_i^{}`$ is the Hermitian conjugate of $`a_i`$. So far we have ignored any coupling of the two oscillators to the environment so as to focus on the interaction of the two oscillators.

In order to perform a QND measurement of $`a_0^{}a_0`$ the Hamiltonian of the oscillators $`H`$ should satisfy the QND condition

$$[a_0^{}a_0,H_0+V_{\mathrm{anh}}]=0.$$
(5)
To show that this condition is satisfied in an appropriate limit it is useful to move into an interaction picture with respect to $`H_0`$. If the frequencies of the two oscillators satisfy $`\omega _0\omega _1\lambda _{01}`$ and $`\omega _i\lambda _{01}`$, then the time-dependent terms in $`V_{\mathrm{anh}}(t)`$ lead to rapid, small amplitude oscillations of $`a_i`$ that essentially average to zero over the time scales for which the non-linearity $`\lambda _{01}`$ is relevant. If we admit a time coarse-graining over times longer than the mechanical oscillation period we may ignore the rapidly oscillating terms, an approximation known as the rotating wave approximation (RWA). Another intuitive explanation for the rotating wave approximation is that the condition $`\omega _0\omega _1\lambda _{01}`$ means that the differences in energy are so large that the energy non-conserving transitions are strongly suppressed.

Disregarding the energy non-conserving terms in the Hamiltonian and then absorbing constant corrections to the system and ancilla oscillation frequency into the definition of $`\omega _0`$ and $`\omega _1`$, we obtain

$$V_{\mathrm{anh}}^{\mathrm{RWA}}(t)=\mathrm{}\lambda _{01}a_0^{}a_0a_1^{}a_1.$$
(6)
The constant term has been disregarded since it merely provides an overall phase. Note that having made the rotating wave approximation the anharmonic coupling term commutes with the observable $`a_0^{}a_0`$, and so a QND measurement can be achieved under the condition $`\omega _0\omega _1\lambda _{01}`$ 111Actually, it is not necessary to make the rotating wave approximation with respect to both oscillators as the interaction $`V_{\mathrm{anh}}`$ becomes QND after the rotating wave approximation made in an interaction picture with respect to $`\mathrm{}\omega _0a_0^{}a_0`$ only. However, we make the further rotating wave approximation to simplify the later development.. Returning to the Schrödinger picture the Hamiltonian $`H`$ now can be written as

$$H^{\mathrm{RWA}}=\mathrm{}\omega _0a_0^{}a_0+\mathrm{}\omega _1a_1^{}a_1+\mathrm{}\lambda _{01}a_0^{}a_0a_1^{}a_1.$$
(7)
In the above rotating wave Hamiltonian an excitation of the system oscillator leads to a frequency shift of the ancilla oscillator. To detect the system excitation level, the ancilla is driven on resonance and the phase shift of the oscillations is measured. The driving of the ancilla may be written in terms of term in the Hamiltonian in the Schrödinger picture

$$H_{\text{drive}}=2\mathrm{}E\left(a_1+a_1^{}\right)\mathrm{cos}\omega _1t,$$
(8)
where the parameter $`E`$ is used to characterize the strength of the drive. In the interaction picture using the RWA for $`\omega _1>E`$, we get

$$H_{\text{drive}}^{\mathrm{RWA}}=\mathrm{}E\left(a_1+a_1^{}\right).$$
(9)
Now we add the coupling of thermal baths to the system and ancilla. We employ a standard technique and model the thermal baths (the surrounding environment) as an infinite number of harmonic oscillators. The thermal baths are linearly coupled to the system or ancilla by coordinate-coordinate coupling, i.e., $`_jA_{ij}x_ix_j`$ where $`x_i`$ is the system or ancilla coordinate, $`x_j`$ is the coordinate of an oscillator in the bath, with the index $`j`$ corresponding to different bath oscillators. We will again use the rotating wave approximation for the coupling since the couplings are weak.

The nature of the coupling with the measurement instrument depends on the measurement scheme. Here we adopt a magnetomotive detection scheme suggested by Yurke et al. YGPB94 ; CR96 ; CR99 . A metallized conducting surface on the ancilla oscillator develops an electromotive force across it due to a perpendicular magnetic field and the oscillation of the beam. The voltage developed depends on

$$V=lB\frac{dx_1}{dt},$$
(10)
where $`V`$ is the voltage, $`B`$ is the magnetic field, the conductor is of length $`l`$ and $`x_1`$ is the displacement of the beam from its equilibrium position. Depending on the resistance $`R`$ of the conducting strip and remainder of the circuit this will result in a current that is then measured.

In order to quantize this measuring device we follow the standard practice in quantum electronics and model this resistance by a semi-infinite transmission lineL . This model has been used in the context of mechanical measurements (see YD ; PM01 ) and is in fact mathematically the same as the ’Rubin model’ R63 ; W . This is certainly a simplified model of an actual detection circuit which essentially assumes that the noise in the circuit is broad-band and Gaussian. The transmission line will be considered to be at a temperature corresponding to the effective noise temperature of the detection circuit and this noise will affect both the sensitivity and the heating of the ancilla and system oscillators. More realistic quantum mechanical models of amplifier circuits can be found in C98  for example. Our final model for the QND set-up will be fairly robust to the precise detection circuit. The resulting current operator is

$$I\underset{n}{}b_{\mathrm{d},n}+b_{\mathrm{d},n}^{},$$
(11)
where $`b_{\mathrm{d},n}`$ are the lowering operators for the modes of the transmission line and the proportionality constant, which is not important for our results, depends on the circuit parameters. For a linearly coupled system-bath measurement within the rotating wave approximation, the Hamiltonian describing the coupling between each measurement current mode and the ancilla is proportional to $`b_{\mathrm{d},n}^{}a_1+b_{\mathrm{d},n}a_1^{}`$. The coupling to the thermal bath modes has the same mathematical structure. In the rotating wave approximation, the difference between a coordinate-coordinate coupling, and a momentum-coordinate coupling can be absorbed into the definition of the phase of the various raising and lowering operators. As is usually done, for later convenience we will include a phase factor of $`\pi /2`$ so that the coupling to the baths and measurement current takes the latter form.

The final Hamiltonian for our model is then

$`H`$ $`=\mathrm{}\omega _0a_0^{}a_0+\mathrm{}E\left(a_1+a_1^{}\right)+\mathrm{}{\displaystyle \underset{\mathrm{s}}{}}{\displaystyle \underset{n}{\overset{\mathrm{}}{}}}\omega _{\mathrm{s},n}b_{\mathrm{s},n}^{}b_{\mathrm{s},n}`$
$`+\mathrm{}\lambda _{01}a_0^{}a_0a_1^{}a_1+i\mathrm{}\left(\mathrm{\Theta }^{}a_0\mathrm{\Theta }a_0^{}\right)+i\mathrm{}\left(\mathrm{\Gamma }^{}a_1\mathrm{\Gamma }a_1^{}\right)+i\mathrm{}\left(D^{}a_1D\left(t\right)a_1^{}\right),`$ (12)
where $`\mathrm{\Gamma },D,\mathrm{\Theta }`$ have the form, $`_\mathrm{n}^{\mathrm{}}g_\mathrm{s}\left(\omega _n\right)b_{\mathrm{s},n},`$ and the index $`\mathrm{s}`$ denotes the three different baths: the thermal bath coupled to the system (0), the thermal bath coupled to the ancilla (1), and the measurement bath coupled to the ancilla (d). The strength of the coupling to the bath modes is given by the coefficients $`g_\mathrm{s}\left(\omega _n\right)`$. Later, we will derive the relationship between these coefficients and the corresponding oscillator damping rates or quality factors.

## III Dynamics of the system

We use the dynamics described by Eq. (12) to understand the measurement process. In this section we first find a master equation that describes the evolution of the system and ancilla alone without explicitly describing the state of the environment. Secondly we further simplify this equation by making use of the difference in time-scales between the system and ancilla to obtain a master equation for the system oscillator alone by means of adiabatic elimination. This allows us to study the effect of the QND measurement coupling on the system.

### III.1 Master equation

We develop a master equation for the density operator of the system alone by integrating out the bath degrees of freedom. Because we are interested in high-Q oscillators weakly coupled to the baths we employ a rotating wave approximation and the Markov approximation that the memory time scale of the bath is short. In this regime the rotating wave master equation accurately describes the dynamics on time-scales longer than an oscillation period, and the resulting master equation preserves the positivity of the density matrix. The derivation of such master equations is widely discussed in the literature, see for example Carmichael C , Walls and Milburn WM  and Caldeira and Leggett CL83 . Note that Caldeira and Leggett do not make the rotating wave approximation (which we adopt from the quantum optics literature) but instead make a high-temperature approximation; the two equations agree in the overlap of their domain of validity (high temperature and weak coupling to the bath). However, the Caldeira-Leggett master equation can only be guaranteed to preserve the positivity of the density operator in the limit of high temperature.

Assuming that the environment and measurement baths are in thermal equilibrium, the master equation for the reduced density operator $`\rho `$ describing the state of the system and ancilla in the interaction picture takes what is known as Lindblad form

$`{\displaystyle \frac{d\rho }{dt}}`$ $`={\displaystyle \frac{i}{\mathrm{}}}[\mathrm{}E\left(a_1^{}+a_1\right)+\mathrm{}\lambda _{01}a_0^{}a_0a_1^{}a_1,\rho ]`$
$`+\nu (N_0+1)𝒟\left[a_0\right]\rho +\nu N_0𝒟\left[a_0^{}\right]\rho +\kappa (N_1+1)𝒟\left[a_1\right]\rho +\kappa N_1𝒟\left[a_1^{}\right]\rho ,`$ (13)
where

$$𝒟\left[x\right]\rho =2x\rho x^{}x^{}x\rho \rho x^{}x,$$

and $`N_\mathrm{i}`$ are the Bose-Einstein factors at frequency $`\omega _0`$ for $`N_0`$ and $`\omega _1`$ for $`N_1`$. The first term in Eq. (13) involving the commutator describes the coherent driving of the ancilla oscillator and the non-linear coupling between the two oscillators in the rotating wave approximation. The remaining terms describe the dissipative interactions with the various baths. The coefficient $`\nu `$ is

$$\nu \pi \varrho _{\mathrm{b0}}\left(\omega _0\right)\left|g_{\mathrm{b0}}\left(\omega _0\right)\right|^2,$$

where $`\varrho _{\mathrm{b0}}\left(\omega _0\right)`$ is the density of states of the bath coupled to the system at frequency $`\omega _0`$. It can be experimentally obtained from the quality factor of the system oscillator $`Q_0`$ as $`\nu =\omega _0/2Q_0`$. The coefficient $`\kappa `$ is the corresponding damping rate of the ancilla, with contributions $`\eta `$ from the coupling to the environment and $`\mu `$ from the measurement bath. Both these rates can be expressed in terms of the bath density of states at frequency $`\omega _1`$ in exactly the same way as for $`\nu `$, and $`\kappa =\eta +\mu `$. The terms containing a factor $`(N_i+1)`$ describe the spontaneous and stimulated emission of phonons into the thermal bath while the ones with $`N_i`$ correspond to absorption of phonons.

The master equation Eq. (13) can in principle be numerically integrated. However we will make some further approximations in order to derive a master equation for the system dynamics alone and show that in some limit the readout system coupling results in the phase diffusion that is required as the backaction for the QND measurement, with no extra noise above this quantum limit. To do this we assume that the ancilla is strongly damped. In this limit the ancilla relaxes rapidly to a state consistent with the instantaneous system state. As a result its dynamics are slaved to the system oscillator and can in fact be eliminated from the equations of motion. Experimentally this is the limit in which the displacement of the ancilla directly reflects the system behavior. This *adiabatic elimination* is described in the following subsection. The final result of this analysis is Eq. (25) below.

Note that $`\nu ,\kappa `$ are the widths of the oscillator resonances, and these should be taken into account when assessing the validity of the rotating wave approximation. In the presence of the coupling to the baths, the rotating wave approximation is only valid if $`\omega _0\omega _1`$ is much greater than the line width of the oscillators, i.e., $`\omega _0\omega _1\nu ,\kappa `$. This condition can be understood as not allowing non-energy-conserving transfers of a phonon between the two oscillators.

### III.2 Adiabatic elimination

For a strongly damped ancilla ($`\kappa \nu `$) the driven ancilla rapidly relaxes to a state that oscillates with a phase determined by the current system phonon number. In the interaction picture this is a displaced thermal state: a state with variance of position and momentum equal to those of a thermal state but with non-zero expectation values of position and momentum consistent with the driving and damping of the oscillator. It will be useful to transform the equations of motion in such a way as to make a perturbative expansion around this steady state. The basic idea is to transform the origin of phase space such that the ancilla steady state for the transformed master equation is a thermal state. This transformation will essentially remove the driving term in the master equation. This is the approach of Wiseman and Milburn WM00  who study adiabatic elimination in a similar model. While they assume zero temperature and therefore end up with a perturbation expansion about the displaced ancilla ground state we must generalize their techniques to finite temperature.

Following Wiseman and Milburn, we use the displacement operator, $`D(\alpha )=\mathrm{exp}\left[\alpha a_1^{}\alpha ^{}a_1\right]`$ with $`\alpha =iE/\kappa `$. The transformed system state is $`\stackrel{~}{\rho }D(\alpha )\rho D(\alpha )^{}`$ and we may write the master equation for $`\stackrel{~}{\rho }`$ as

$`\dot{\stackrel{~}{\rho }}`$ $`=D(\alpha )\dot{\rho }D(\alpha )^{}`$
$`=i\left|\alpha \right|^2\lambda _{01}[a_0^{}a_0,\stackrel{~}{\rho }]`$
$`i\lambda _{01}[a_0^{}a_0a_1^{}a_1,\stackrel{~}{\rho }]i\lambda _{01}[a_0^{}a_0\left(\alpha a_1^{}+\alpha ^{}a_1\right),\stackrel{~}{\rho }]`$
$`+\kappa (N_1+1)\left(2a_1\stackrel{~}{\rho }a_1^{}a_1^{}a_1\stackrel{~}{\rho }\stackrel{~}{\rho }a_1^{}a_1\right)+\kappa N_1\left(2a_1^{}\stackrel{~}{\rho }a_1a_1a_1^{}\stackrel{~}{\rho }\stackrel{~}{\rho }a_1a_1^{}\right)`$
$`+\nu (N_0+1)\left(2a_0\stackrel{~}{\rho }a_0^{}a_0^{}a_0\stackrel{~}{\rho }\stackrel{~}{\rho }a_0^{}a_0\right)+\nu N_0\left(2a_0^{}\stackrel{~}{\rho }a_0\stackrel{~}{\rho }a_0a_0^{}a_0a_0^{}\stackrel{~}{\rho }\right).`$ (14)
In this master equation the excitation of the ancilla oscillations leads to a frequency shift of the system oscillator described by the first three terms on the right hand side of this equation. The first term is due to the classical mean value of the ancilla oscillator energy and is just a constant shift in the system oscillation frequency. We may move to an interaction picture at this shifted frequency, the most convenient interaction picture in which to perform the adiabatic elimination. The next two terms describe the effect of the fluctuations in the ancilla excitation. The thermal bath coupling terms (the last four groups of terms) are the same as before.

The adiabatic elimination will hold when the terms proportional to $`\kappa `$ in Eq. (14) dominate in the ancilla dynamics. Thus, the adiabatic elimination is valid in a strongly damped regime such that

$$\frac{\lambda _{01}\left|\alpha \right|}{\kappa },\frac{\nu }{\kappa }ϵ1.$$
(15)
We are assuming that the ancilla oscillator relaxes faster than the system oscillator as well as that the non-linear dynamics are weak compared to the damping of the ancilla oscillator. For the consistency of the following treatment it will also be necessary to have $`\lambda _{01}N_1/\kappa ϵ^2`$. This requirement follows from the second term (the non-linear coupling term) on the right hand side of the master equation. This constraint can be achieved consistent with Eq. (15), for example, by leaving $`\alpha `$ finite and choosing $`N_1,\lambda _{01}/\kappa ϵ`$, a regime of low temperature and moderate non-linearity. The approximations are also valid at arbitrary temperature in the limit of strong driving and weak non-linearity such that $`\lambda _{01}/\kappa ϵ^2`$ and $`\alpha ϵ^1`$ hold 222However, the driving strength is chosen in such a way to avoid the bi-stability region caused by the anharmonicity of the oscillator, see VI.. Here the scaling of the driving strength is chosen to preserve the measurement sensitivity which will scale with $`\lambda _{01}\alpha /\kappa `$. In this regime the frequency shift of the system oscillator becomes large.

As mentioned above, in this displaced frame the state of the readout oscillator is close to a thermal state and we expand $`\stackrel{~}{\rho }`$ in the form

$`\stackrel{~}{\rho }`$ $`=\rho _0\rho _{N_1}+\rho _1a_1^{}\rho _{N_1}+\rho _1^{}\rho _{N_1}a_1+\rho _2a_1^{}\rho _{N_1}a_1`$
$`+\rho _2^{}a_1^2\rho _{N_1}+\rho _2^{}^{}\rho _{N_1}a_1^2+O\left(ϵ^3\right).`$ (16)
Here the $`\rho _i,i=0,1,2\mathrm{}`$ act on the system oscillator and the subscripts indicate orders of magnitude in $`ϵ`$. The scalings of the different parameters with $`ϵ`$ have been chosen to guarantee the consistency of the expansion. The quantity $`\rho _{N_1}`$ is the thermal density matrix for the ancilla, which in terms of the average excitation number $`N_1`$ is

$$\rho _{N_1}=\underset{n=0}{}\frac{1}{N_1+1}\left(\frac{N_1}{N_1+1}\right)^n|nn|.$$
(17)
This is the steady state of the master equation for an oscillator coupled to a thermal bath with temperature given by $`N_1`$ We have restricted Eq. (16) to normal ordered terms using the following identities which can be proved from this expression for $`\rho _{N_1}`$:

$`\rho _{N_1}a_1^{}`$ $`={\displaystyle \frac{N_1}{N_1+1}}a_1^{}\rho _{N_1},`$ (18)
$`a_1\rho _{N_1}`$ $`={\displaystyle \frac{N_1}{N_1+1}}\rho _{N_1}a_1.`$ (19)
These normal ordering identities are the key to generalizing the arguments of Wiseman and Milburn to the finite temperature case. Using $`\mathrm{Tr}_1(\rho _{N_1})=1`$ and $`\mathrm{Tr}_1(a_1^{}\rho _{N_1}a_1)=N_1+1`$, it can be seen that the system density matrix after tracing out the ancilla state is

$$\rho _\mathrm{s}=\mathrm{Tr}_1\left\{\stackrel{~}{\rho }\right\}=\rho _0+(N_1+1)\rho _2.$$
(20)
Now we substitute Eq. (16) into Eq. (14) and, using the normal ordering identities, derive equations of motion for the operators $`\rho _i`$, retaining terms in the evolution of $`\rho _0`$ and $`\rho _2`$ up to second order in $`ϵ`$:

$`\dot{\rho }_0=`$ $`i\lambda _{01}\left[\alpha ^{}a_0^{}a_0\rho _1\alpha \rho _1^{}a_0^{}a_0\right]+2\kappa (N_1+1)\rho _2+_0\rho _0+\kappa 𝒪\left(ϵ^3\right),`$ (21)
$`\dot{\rho }_1=`$ $`i\lambda _{01}a_0^{}a_0\alpha \rho _0+i\lambda _{01}\alpha {\displaystyle \frac{N_1}{N_1+1}}\rho _0a_0^{}a_0\kappa \rho _1+\kappa 𝒪\left(ϵ^2\right),`$ (22)
$`\dot{\rho }_2=`$ $`i\lambda _{01}a_0^{}a_0\left[\alpha \rho _1^{}+\alpha ^{}{\displaystyle \frac{N_1}{N_1+1}}\rho _1\right]+i\lambda _{01}\left[\alpha {\displaystyle \frac{N_1}{N_1+1}}\rho _1^{}+\alpha ^{}\rho _1\right]a_0^{}a_0`$
$`2\kappa \rho _2i\lambda _{01}{\displaystyle \frac{N_1}{N_1+1}}[a_0^{}a_0,\rho _0]+\kappa 𝒪\left(ϵ^3\right),`$ (23)
where $`_0\rho _0`$ refers to the damping of the system oscillator described by the last two terms in Eq. (14). When $`\kappa `$ is large, the equation for $`\rho _1`$ is strongly damped and quickly decays to the steady state. So we perform adiabatic elimination by setting $`\dot{\rho }_1=0`$ and obtaining an expression for $`\rho _1`$

$$\rho _1=i\frac{\lambda _{01}}{\kappa }\left[\alpha a_0^{}a_0\rho _0\alpha \frac{N_1}{N_1+1}\rho _0a_0^{}a_0\right]+𝒪\left(ϵ^2\right).$$
(24)
Substituting Eq. (24) into Eqs. (21, 23) and using the definition of the reduced density matrix Eq. (20) we find, up to second order in $`ϵ`$, the master equation for the reduced density matrix

$`\dot{\rho }_\mathrm{s}`$ $`={\displaystyle \frac{\lambda _{01}^2\left|\alpha \right|^2(2N_1+1)}{\kappa }}[a_0^{}a_0,[a_0^{}a_0,\rho _\mathrm{s}]]i\left\{\omega _0+\lambda _{01}(\left|\alpha \right|^2+N_1)\right\}[a_0^{}a_0,\rho _\mathrm{s}]`$
$`+\nu \left(N_0+1\right)\left(2a_0\rho _\mathrm{s}a_0^{}a_0^{}a_0\rho _\mathrm{s}\rho _\mathrm{s}a_0^{}a_0\right)+\nu N_0\left(2a_0^{}\rho _\mathrm{s}a_0\rho _\mathrm{s}a_0a_0^{}a_0a_0^{}\rho _\mathrm{s}\right).`$ (25)
This is the main result of this section. Note that the effect of the adiabatic elimination has essentially been to replace $`a_1`$ by $`\lambda _{01}|\alpha |a_0^{}a_0/\kappa `$, an indication that by measuring the ancilla oscillations it will be possible to obtain information about the system phonon number.

## IV Quantum trajectories

Equation (25) describes the statistical behavior of the system due to the coupling to the thermal bath and the indirect coupling to the ancilla thermal bath and the measurement bath, but does not tell us how the measured current reflects the system state, or about the correlations of the system dynamics with particular measurement outcomes. In this section we derive an equation of motion for the state of the system conditioned on a particular sequence of measurement outcomes. This equation is termed a quantum trajectoryB86 ; B87 ; C93 ; Wth , and results from continually projecting onto eigenstates of the current. Since the current effectively measures phonon number, this measurement process will tend to force the system towards a pure number state that is consistent with the measurement current. The time scale for this to occur will depend on the coupling of the system to the measurement apparatus, which is in turn connected to the sensitivity of the measurement. On the other hand, the coupling of the system to a thermal bath will lead it to absorb and emit energy from the bath. Thus, in order to determine which number state the system is in and track its evolution, it must be possible to distinguish between one number state and the next in a time that is short compared to the time scale over which phonons are absorbed from and emitted into the thermal bath.

A quantum trajectory is constructed as follows. Over each infinitesimal time interval, the system and the measurement bath states become weakly entangled via the interaction Hamiltonian. As a consequence, at each time instant, the state of the system influences the distribution of the possible values of the current $`I`$ that may be obtained in the measurement. In turn, von Neumann projections of the entangled states allow us to calculate the effect of the measurement of the current on the system state. The appropriate projection is onto the current eigenvector corresponding to the measured current value. This results in a stochastic master equation for the state evolution. To implement the quantum trajectory approach we perform a simulation by picking the measurements $`I(t)`$ from the correct probability distribution and following the corresponding evolution of the system state. The $`I(t)`$ curve produced by such a simulation is representative of a single experimental run, and is a useful predictor of what the experimentalist might see. There will be a signal contribution that reflects the system state, as well as a white noise background arising from both thermal and quantum noise.

### IV.1 Description of the measurement

While quantum trajectories are discussed in general at zero temperature in the quantum optics literature, Wiseman has discussed the quantum trajectory equations for homodyne detection at finite temperature (Ref. Wth , § 4.4.1). The demodulated current that reflects changes in the phase of the ancilla oscillation in our setup is mathematically analogous to homodyne detection at finite temperature, and so we can adopt these results here.

The measurement bath is described by the Boson operators $`b_{d,n}`$. Since the measurement bath is assumed to be large, the finely spaced modes with a smooth density of states leads to a short memory time, a result known as the Markov limit. To exploit this is is useful to introduce a global bath operator which captures the combination of bath modes that interact with the ancilla oscillating at frequency $`\omega _1`$ at time $`t`$ (see Appendix A for the derivation of these results)

$$B_t=\frac{1}{\sqrt{2\pi \rho _\mathrm{d}(\omega _1)}g_\mathrm{d}\left(\omega _1\right)}\underset{n}{}g_\mathrm{d}\left(\omega _n\right)b_{\mathrm{d},n}e^{\iota \left(\omega _n\omega _1\right)t},$$
(26)
and has time-local commutation rules in the Markov approximation

$$[B_t,B_t^{}^{}]=\delta \left(tt^{}\right).$$
(27)
The operator $`B_t`$ should be considered to be a linear combination of Schrödinger picture operators, with the phase factors of the coefficients depending on the parameter $`t`$ GCZ92 . In quantum optics this is termed the input field operator and, roughly speaking, describes the combination of bath modes interacting with the system at time $`t`$. In terms of these operators the current Eq. (11) (appropriately scaled to remove proportionality constants) is

$$I(t)=B_t+B_t^{},$$
(28)
and the interaction Hamiltonian between the ancilla and the measurement bath in the interaction picture is

$$H_{\mathrm{int}}^\mathrm{I}(t)=i\mathrm{}\sqrt{2\mu }(B_ta_1^{}B_t^{}a_1).$$
(29)
with $`\mu =\pi \varrho _\mathrm{d}\left(\omega _1\right)\left|g_\mathrm{d}\left(\omega \right)\right|^2`$ the ancilla damping rate coming from the measurement bath coupling as before.

The idea of the calculation is to consider the interaction of the ancilla with the bath at time $`t`$, represented by the operator $`B_t`$, over a small time interval $`\mathrm{\Delta }t`$. It is supposed that each “element” in the time sequence of the bath $`B_t`$ is initially described by a thermal state. Over the interval $`\mathrm{\Delta }t`$ the ancilla and bath states become weakly entangled. Measurement of the current (i.e., the bath operator $`B_t+B_t^{}`$) then finds a value of the current equal to an eigenvalue $`I`$ of the current operator, with the corresponding eigenstate $`|I`$, with a probability distribution $`P(I)`$ given by the density matrix of the entangled state in the usual way.

$$P(I)=I\left|\rho \left(t+\mathrm{\Delta }t\right)\right|I.$$
(30)
The measurement also projects the density matrix onto the eigenstate $`|I`$

$$\rho \frac{|II\left|\rho \left(t+\mathrm{\Delta }t\right)\right|II|}{I\left|\rho \left(t+\mathrm{\Delta }t\right)\right|I}.$$
(31)
Since the value of the current measured is a stochastic variable, this projection adds a stochastic component to the evolution of the density matrix.

To follow the evolution over a time $`\mathrm{\Delta }t`$ it is useful to introduce the normalized operator

$$\mathrm{\Delta }B=\left[_0^{\mathrm{\Delta }t}B_t𝑑t\right]/\sqrt{\mathrm{\Delta }t}B_t\sqrt{\mathrm{\Delta }t},$$
(32)
which satisfies the commutation rule

$$[\mathrm{\Delta }B\left(t\right),\mathrm{\Delta }B^{}\left(t\right)]=1.$$
(33)
At time $`t`$ the density matrix representing the ancilla and the segment of the measurement bath represented by $`\mathrm{\Delta }B\left(t\right)`$ can be written as a direct product of the system plus ancilla $`\rho \left(t\right)`$ and bath $`\rho _\mathrm{b}\left(t\right)`$ density matrices

$$\overline{\rho }\left(t\right)=\rho \left(t\right)\rho _\mathrm{b}\left(t\right),$$
(34)
and $`\rho _\mathrm{b}\left(t\right)`$ is a thermal state. To lowest order the evolution under the interaction Eq. (29) gives

$$\overline{\rho }\left(t+dt\right)=\rho \left(t\right)\rho _\mathrm{b}\left(t\right)+\sqrt{2\mu }\sqrt{\mathrm{\Delta }t}[\mathrm{\Delta }B^{}a_1a_1^{}\mathrm{\Delta }B,\rho \left(t\right)\rho _\mathrm{b}\left(t\right)]+𝒪\left(\mathrm{\Delta }t\right).$$
(35)
The second term on the right hand side of this equation is the leading order term in the weak entangling of the state, and will lead, after projection, to the stochastic part of the density matrix evolution. Using the $`\mathrm{\Delta }B`$ notation has made the $`𝒪\left(\sqrt{\mathrm{\Delta }t}\right)`$ size of this term explicit. To derive the deterministic part of the evolution equation we would need to keep the $`𝒪\left(\mathrm{\Delta }t\right)`$ terms, but since these are already known (the master equation in Lindblad form) we will not do this here.

The scheme is now to project this density matrix onto an eigenstate of $`\mathrm{\Delta }B+\mathrm{\Delta }B^{}`$ chosen with a probability given by $`\overline{\rho }\left(t+dt\right)`$. Because of the weak coupling of the bath with the system, this will give a small additional contribution (actually proportional to $`\sqrt{\mathrm{\Delta }t}`$) to the system density matrix depending on the value of the current measured. Since the combination of operators $`\mathrm{\Delta }B+\mathrm{\Delta }B^{}`$ is just the displacement $`X`$ of the harmonic oscillator represented by the operator $`\mathrm{\Delta }B`$, this projection is most easily done by first transforming the state of the bath into a Wigner function representation (see for example Ga ). At time $`t`$, the bath oscillator described by $`\mathrm{\Delta }B(t)`$ is in a thermal state and the distribution of $`X`$ is a Gaussian centered at $`X=0`$ and with width $`2N_1+1=\mathrm{coth}\left(\mathrm{}\omega _1/2k_\mathrm{B}T\right)`$. Following the evolution of the state shows that at time $`t+\mathrm{\Delta }t`$ and to $`𝒪\left(\sqrt{\mathrm{\Delta }t}\right)`$ the distribution of $`X`$ after the evolution corresponding to the operation Eqs. (31) and (35) remains Gaussian and with the same width, but now centered around $`\sqrt{2\mu }\mathrm{Tr}_{\rho \left(t\right)}\left\{a_1+a_1^{}\right\}\mathrm{\Delta }t`$. This means that the variable $`\sqrt{\mathrm{\Delta }t}X`$ is a Gaussian random variable given by

$$\sqrt{\mathrm{\Delta }t}X=\sqrt{2\mu }a_1+a_1^{}\mathrm{\Delta }t+\sqrt{2N_1+1}dW,$$
(36)
with $`dW`$ a Wiener increment with $`dW^2=\mathrm{\Delta }t`$. These results give us expressions for the measured current and the effect of the measurement on the system density matrix.

Since the current is $`X/\sqrt{\mathrm{\Delta }t}`$, the first important result is that the measured current integrated over time $`\mathrm{\Delta }t`$ is

$$I\left(t\right)\mathrm{\Delta }t=\sqrt{2\mu }a_1+a_1^{}\mathrm{\Delta }t+\sqrt{2N_1+1}dW,$$
(37)
or in differential form

$$I\left(t\right)=\sqrt{2\mu }a_1+a_1^{}\left(t\right)+\sqrt{2N_1+1}\xi \left(t\right),$$
(38)
where $`\xi \left(t\right)=dW/dt`$ represents white noise with correlations

$`\xi \left(t\right)`$ $`=0,`$ (39)
$`\xi \left(t\right)\xi \left(t^{}\right)`$ $`=\delta \left(tt^{}\right).`$ (40)
The second result is for the increment of the system density matrix after evolution through $`\mathrm{\Delta }t`$ and projection by the measurement (cf. Wth  Eq. (4.113)),

$`d\rho ^{\mathrm{st}}\left(t\right)`$ $`=X|\overline{\rho }(t+\mathrm{\Delta }t)|X/p(X)`$
$`=\sqrt{\mathrm{\Delta }t}X{\displaystyle \frac{\sqrt{2\mu }}{2N_1+1}}\left[\left(N_1+1\right)\left(a_1\rho ^{\mathrm{st}}+\rho ^{\mathrm{st}}a_1^{}\right)N_1\left(a_1^{}\rho ^{\mathrm{st}}+\rho ^{\mathrm{st}}a_1^{}\right)\mathrm{Tr}\left\{a_1\rho ^{\mathrm{st}}+\rho ^{\mathrm{st}}a_1^{}\right\}\rho ^{\mathrm{st}}\right]+𝒪\left(\mathrm{\Delta }t\right).`$ (41)
Replacing the stochastic variable $`X`$ by the expression Eq. (36), and retaining only the $`𝒪\left(\sqrt{\mathrm{\Delta }t}\right)`$ term gives

$$d\rho ^{\mathrm{st}}=\sqrt{\frac{2\mu }{1+2N_1}}\left[(N_1+1)(a_1\rho ^{\mathrm{st}}+\rho ^{\mathrm{st}}a_1^{})N_1(a_1^{}\rho ^{\mathrm{st}}+\rho ^{\mathrm{st}}a_1)a_1+a_1^{}\rho ^{\mathrm{st}}\right]dW+𝒪\left(\mathrm{\Delta }t\right).$$
(42)
Equation (42) is the stochastic term that must be added to the density matrix evolution of equation (13) to give us the stochastic master equation for the density matrix conditioned on the measurement outcome. Note that the noise term $`dW`$ appearing in Eq. (42) is the same as that appearing in Eq. (37), so that it is related to the actual current measured $`I(t)`$ by Eq. (37)

$$dW=\left(I(t)\sqrt{2\mu }a_1+a_1^{}\right)\mathrm{\Delta }t/\sqrt{1+2N_1}.$$
(43)
### IV.2 Adiabatic elimination on the stochastic master equation

Just as we did for the master equation it is possible to adiabatically eliminate the ancilla coordinates and find a stochastic master equation for the system alone. Using the same expansion for the conditioned density matrix Eq. (16) we can determine stochastic equations for $`\rho _i^{\mathrm{st}}`$ from Eq. (42). We obtain the set of differential equations

$`d\rho _0^{\mathrm{st}}`$ $`=\sqrt{{\displaystyle \frac{2\mu }{2N_1+1}}}\left\{(N_1+1)\left(\rho _1^{\mathrm{st}}+\rho _1^{\mathrm{st}}\right)\rho _1^{\mathrm{st}}+\rho _1^{\mathrm{st}}\rho _0\right\}dW+𝒪(dt),`$ (44)
$`d\rho _1^{\mathrm{st}}`$ $`=\sqrt{{\displaystyle \frac{2\mu }{2N_1+1}}}\left\{\left(N_1+1\right)\left(2\rho _2^{}^{\mathrm{st}}+\rho _2^{\mathrm{st}}\right)\rho _1^{\mathrm{st}}+\rho _1^{\mathrm{st}}\rho _1^{\mathrm{st}}\right\}dW,`$ (45)
$`d\rho _2^{\mathrm{st}}`$ $`=\sqrt{{\displaystyle \frac{2\mu }{2N_1+1}}}\left\{\rho _1^{\mathrm{st}}+\rho _1^{\mathrm{st}}\rho _2^{\mathrm{st}}\right\}dW.`$ (46)
We have written only the stochastic contributions; the terms proportional to $`dt`$ are the same as in the adiabatic elimination on the ordinary master equation.

Now we will do the adiabatic elimination as we did for the deterministic master equation. As before we wish to set $`d\rho _1^{\mathrm{st}}`$ to zero. However, since $`d\rho _1^{\mathrm{st}}`$ is stochastically driven this will not be precisely true even in the steady state. In order to replace $`\rho _1^{\mathrm{st}}`$ by its mean value at long times it is necessary to estimate the size of the fluctuations resulting from the $`dW`$-term. Following the analysis of Doherty and Jacobs DJ99  we integrate the stochastic term over the decay time of the ancilla oscillator and compare its root mean square magnitude with the deterministic terms. Consider the full equation for $`d\rho _1`$

$`d\rho _1^{\mathrm{st}}=`$ $`i\lambda _{01}\alpha \left[a_0^{}a_0\rho _0^{\mathrm{st}}{\displaystyle \frac{N_1}{N_1+1}}\rho _0^{\mathrm{st}}a_0^{}a_0\right]dt`$
$`\kappa \rho _1dt+\sqrt{{\displaystyle \frac{2\mu }{2N_1+1}}}\left\{\left(N_1+1\right)\left(2\rho _2^{}^{\mathrm{st}}+\rho _2\right)\rho _1^{\mathrm{st}}+\rho _1^{\mathrm{st}}\rho _1^{\mathrm{st}}\right\}dW+\kappa 𝒪\left(ϵ^3\right)dt.`$ (47)
We integrate this over a time $`\mathrm{\Delta }t1/\kappa `$ and use the fact that the mean values of $`\rho _0,\rho _1`$ must be slowly varying over this time to obtain

$`\mathrm{\Delta }\rho _1`$ $`={\displaystyle _0^{\mathrm{\Delta }t}}𝑑\rho _1`$ (48)
$``$ $`i\lambda _{01}\alpha \left[a_0^{}a_0\rho _0{\displaystyle \frac{N_1}{N_1+1}}\rho _0a_0^{}a_0\right]\mathrm{\Delta }t\kappa \rho _1\mathrm{\Delta }t`$
$`+\sqrt{{\displaystyle \frac{2\mu }{2N_1+1}}}\left\{\left(N_1+1\right)\left(2\rho _2^{}+\rho _2\right)\rho _1+\rho _1^{}\rho _1\right\}{\displaystyle _0^{\mathrm{\Delta }t}}𝑑W(t^{})+𝒪\left(ϵ^3\right).`$ (49)
The random number $`\mathrm{\Delta }W=_0^{\mathrm{\Delta }t}𝑑W(t^{})`$ is Gaussian distributed with mean zero and variance $`\mathrm{\Delta }t1/\kappa `$ Ga2 , thus the root mean square size of $`\mathrm{\Delta }W`$ is $`1/\sqrt{\kappa }`$. As a result the stochastic term in the update of $`\rho _1`$ scales like $`ϵ^{5/2}`$ and is negligible in comparison to the deterministic terms which scale like $`ϵ`$. As a result Eq. (24) holds exactly as before. Using Eq. (20) and $`i\alpha =i\alpha ^{}=|\alpha |`$ we finally obtain the stochastic master equation (SME) for the system

$`d\rho _\mathrm{s}`$ $`=\left\{{\displaystyle \frac{\lambda _{01}^2\left|\alpha \right|^2(2N_1+1)}{\kappa }}[a_0^{}a_0,[a_0^{}a_0,\rho _\mathrm{s}]]\right\}dti\left\{\omega _0+\lambda _{01}(\left|\alpha \right|^2+N_1)\right\}[a_0^{}a_0,\rho _\mathrm{s}]dt`$
$`+\nu \left(N_0+1\right)\left(2a_0\rho _\mathrm{s}a_0^{}a_0^{}a_0\rho _\mathrm{s}\rho _\mathrm{s}a_0^{}a_0\right)dt+\nu N_0\left(2a_0^{}\rho _\mathrm{s}a_0\rho _\mathrm{s}a_0a_0^{}a_0a_0^{}\rho _\mathrm{s}\right)dt`$
$`\sqrt{2k}\left[a_0^{}a_0\rho _\mathrm{s}+\rho _\mathrm{s}a_0^{}a_02a_0^{}a_0\rho _\mathrm{s}\right]dW,`$ (50)
where

$$k\mu \lambda _{01}^2|\alpha |^2/\left(2N_1+1\right)\kappa ^2.$$
(51)
Again the noise $`dW`$ is related to the measured current, which using Eq. (24) can be written in terms of the system phonon number

$$I(t)\mathrm{\Delta }t=\sqrt{2N_1+1}\left(2\sqrt{2k}a_0^{}a_0\mathrm{\Delta }t+dW\right).$$
(52)
## V Results

To further understand the consequences of the stochastic density matrix equation Eq. (50) we consider a case in which the initial state is a mixture of number states, $`\rho _s=_np_n|nn|`$. (A thermal state is an example of such a state.) The solution of the stochastic master equation, Eq. (50), remains a mixture of number states if the initial state is a mixture of number states. For such an initial state the stochastic master equation can be reduced to an equation for the weights $`p_n`$, which takes the form

$`dp_n`$ $`=2\nu (N_0+1)[np_n(n+1)p_{n+1}]dt`$
$`2\nu N_0[(n+1)p_nnp_{n1}]dt2\sqrt{2k}\left(n{\displaystyle \underset{n^{}}{}}n^{}p_n^{}\right)p_ndW.`$ (53)
Since mixtures of number states are invariant under changes of phase and the number states are eigenstates of the Hamiltonian, neither the phase diffusion term nor the Hamiltonian terms in the stochastic master equation contribute to the evolution of the phonon number distribution. It can also be shown that this system of equations also describes the evolution of the phonon number distribution $`p_n=n|\rho |n`$ for an arbitrary (not necessarily diagonal) initial state $`\rho (0)`$.

Equation (53) is our central result for analyzing the behavior of the measurement protocol. The first two terms of Eq. (53) containing $`dt`$ describe emission into and absorption from the thermal bath coupled to the system. We will see that the second, stochastic, term tends to concentrate the distribution $`p(n)`$ onto a single value of $`n`$. The two effects are characterized by the two time scales, $`\nu ^1`$ and $`k^1`$, and the resulting behavior depends on the ratio of these two times. We will discuss the competition by calculating the occupation number of the system, $`a_0^{}a_0\left(t\right)`$, which is given in terms of $`p_n`$ from numerical simulations of Eq. (53) by

$$a_0^{}a_0\left(t\right)=\underset{n}{}np_n\left(t\right).$$
(54)
Firstly we turn off the stochastic component and consider the solutions given by the deterministic part of Eq. (53). Figure 2 is a plot for $`k=0`$ starting from two different initial states $`|1`$ and $`|2`$ and for a bath temperature corresponding to an average occupation number $`N_0=1.62`$. The plot shows that the deterministic terms in Eq. (53) drive the system towards a mixed (thermal) state, so that the ensemble average of $`a_0^{}a_0(t)`$ gradually reaches the thermal average at the bath temperature. This is true regardless of the initial state. Note that the deterministic part of Eq. (53) also describes the average over all measurement outcomes even for nonzero $`k,`$since the stochastic term averages to zero. We can define the characteristic time the system resides in a given number state, which we call the dwell time $`t_{\mathrm{dwell}}`$, as the reciprocal of the initial transition rate given by Eq. (53) with $`k=0`$

$$t_{\mathrm{dwell}}=\frac{dp_n}{dt}|_{p_n(0)=1}=\frac{1}{2\nu \left[N_0(n+1)+\left(N_0+1\right)n\right]}.$$
(55)
Note that the dwell time depends on the initial state $`n`$, and also on the temperature of the bath through $`N_0`$.

We turn now to the dynamics resulting from the measurement process in the absence of coupling to the thermal bath. Figure 3 shows results for $`a_0^{}a_0\left(t\right)`$ for a simulation of Eq. (53) with $`\nu =0`$ and an initial condition of a thermal state. Figure 4 shows the individual probabilities $`p_n`$ for $`n=0,1,2,3`$ for the same simulation. All number states are present initially, but eventually the system is projected onto state $`|1`$ in this simulation. In other runs, with different random numbers for the stochastic term, different final states result, as expected. The plots show that the stochastic term tends to project the system state onto a pure number state on a time scale of order $`k^1`$. We call this time the collapse time $`t_{\text{coll}}`$. Since no coupling to the thermal bath is present in these simulations, once projected onto a number state, the state is stationary. The collapse onto a number state can actually be shown analytically using the solution of the system of equations (53) due to Jacobs and Knight Ja1998 .

For the phonon number $`a_0^{}a_0\left(t\right)`$ to take on discrete values with both the thermalization by the coupling to the bath and the projection by the measurement process present, we need $`t_{\mathrm{dwell}}t_{\mathrm{coll}}`$. This is illustrated in Figs. 5, which show results for the cases $`t_{\mathrm{dwell}}t_{\mathrm{coll}}`$ and the $`t_{\mathrm{dwell}}t_{\mathrm{coll}}`$ with the fixed value of $`N_0=1.62`$. We use values of $`k/\nu =250`$ giving $`t_{\mathrm{dwell}}/t_{\mathrm{coll}}=153`$ for state $`|0`$ and $`t_{\mathrm{dwell}}/t_{\mathrm{coll}}=42.4`$ for state $`|1`$, and $`k/\nu =5`$ giving $`t_{\mathrm{dwell}}/t_{\mathrm{coll}}=3.06`$ for state $`|0`$ and $`t_{\mathrm{dwell}}/t_{\mathrm{coll}}=0.85`$ for state $`|1`$. The jumps in the occupation number are clearly evident in the former case, but are not seen in the latter case. The discreteness in the phonon number is shown more clearly by plotting histograms of $`a_0^{}a_0(t)`$, Figs. 6, again using a fixed value of $`N_0=1.62`$ but with different values of $`k/\nu `$ equal to $`150`$ and $`15`$. (A bin width $`\mathrm{\Delta }a_0^{}a_0=0.1`$ is used.) The clustering of the $`a_0^{}a_0(t)`$ values around integral values is clearly evident for $`k/\nu =150`$, is still identifiable for $`k/\nu =15`$, and completely absent for $`k/\nu =3`$. The increasing sharpness of the jumps with larger $`k/\nu `$ can be seen in a more quantitative manner by plotting the standard deviation of the phonon number from integer values, the time and ensemble average of $`\left|a_0^{}a_0(t)\mathrm{Int}a_0^{}a_0(t)\right|^2`$, as a function of $`k/\nu `$ (see Fig. 7).

Since $`t_{\mathrm{dwell}}`$ is dependent on the temperature, the condition $`t_{\mathrm{dwell}}t_{\mathrm{coll}}`$ effectively places a limit on the temperature of the system oscillator even for large $`k`$. Setting $`n=N_0`$ in Eq. (55) for $`t_{\mathrm{dwell}}`$, this inequality gives the condition on the temperature for jumps in the number to be seen

$$\nu N_0(N_0+1)/2k1.$$
(56)
In order to keep the same resolution for observing clear jumps as at low temperature, $`k/\nu `$ must be increased as temperature increases. This is not an easy task for the experimenters: for example, an oscillator with $`1`$ GHz resonant frequency, which is the highest frequency currently reported for a mesoscopic oscillator HZMR03 , at $`T=0.1`$ K the average occupation number is $`N_0=1.62`$. When the temperature is raised to $`T=1`$ K, the value rises to $`N_0=20`$. Thus if we demand the same resolution for jumping in both cases, the sensitivity of the measurement at the higher temperature must be increased by a large factor. This is illustrated by Figs. 8 which show $`a_0^{}a_0(t)`$ over time for different temperatures corresponding to $`N_0=1.62`$ for $`k/\nu =150`$, and $`N_0=20`$ for $`k`$/$`\nu =1850`$. The product $`\nu N_0/k`$ has been kept constant at $`0.0108`$ in order to provide the same resolution for the jumps. Also notice from Eq. (55) that $`t_{\mathrm{dwell}}`$ decreases with the system state $`n`$ making it difficult to recognize the discrete jumps when the system state is at higher $`n`$.

We have so far considered the possibility of observing discrete occupation numbers in terms of the behavior of the variable $`a_0^{}a_0(t)`$. In actual practice the occupation number must be inferred from the measured current $`I(t)`$, and is obscured by the noise in this variable. A simple scheme to reduce the effect of the noise is to average the signal over a sliding window. We can define the measurement time $`t_\mathrm{m}`$ as the averaging time required to give unit signal to noise ratio. Thus we equate the signal $`S`$, given by averaging the current for unit phonon number $`a_0^{}a_0=1`$ over the measurement time

$$S=_0^{t_\mathrm{m}}𝑑\stackrel{~}{Q}=2\sqrt{2N_1+1}\sqrt{2k}t_\mathrm{m},$$
(57)
with the noise $`N`$ over this averaging time

$$N=\sqrt{_0^{t_\mathrm{m}}𝑑\stackrel{~}{Q}^2}=\sqrt{2N_1+1}t_\mathrm{m}^{1/2}.$$
(58)
Setting $`S/N=1`$ gives the measurement time

$$t_\mathrm{m}=\frac{1}{8k}.$$
(59)
For jumps in the measured phonon number to be detected in the current we would need $`t_{\mathrm{dwell}}t_\mathrm{m}`$. Notice that the measurement time and the collapse time are comparable. This means that if the experimenter can infer the system number state through the measurement current, then the system is actually projected to that state on the same time scale. The results for different averaging times $`\mathrm{\Delta }t`$ is shown in Fig. 9. For $`k\mathrm{\Delta }t`$ equal to 4.5 or 7.5 the averaging is sufficient to display the steps in $`a_0^{}a_0\left(t\right)`$ without too much rounding of the transitions.

The simple averaging is not actually the optimal way to extract $`a_0^{}a_0\left(t\right)`$ from $`I(t)`$. In principle a better approach is to use the stochastic master equation to reconstruct the dynamics of the system given the initial ensemble $`\rho _\mathrm{s}(t_0)`$ and the measured current $`I(t)`$. This can be seen more readily if we rewrite Eq. (52) as

$$dW=\frac{1}{\sqrt{2N_1+1}}\left\{I\left(t\right)2\sqrt{2k}a_0^{}a_0\left(t\right)\right\}dt.$$
(60)
In our simulations we draw $`I(t)`$ at random from the appropriate distribution and find the stochastic density matrix $`\rho _\mathrm{s}`$. However, using $`I(t)`$ from experimental data, the experimenter can in principle propagate Eq. (53) using Eq. (60) and then can estimate the phonon number at each time from Eq. (54). This procedure is itself a low pass filtering that reduces the noise on the measurement current, corresponding to the optimal filtering for our model of the systemb99 ; d99

## VI Parameters and constraints

The results of the previous section show that a large value of the ratio $`k/\nu `$ is crucial. To analyze the interplay of the parameters of an experimental realization, we simplify the expression for $`k/\nu `$ by assuming that most of the damping of the ancilla comes from the necessary coupling to the measurement device, rather than from the extra thermal bath, i.e. $`\mu \kappa `$. Then using the expression of $`k`$ from Eq. (51) and $`\nu =\omega _0/2Q_0`$ we get

$$\frac{k}{\nu }4(2N_1+1)^1Q_0Q_1\frac{\omega _1}{\omega _0}\left(\frac{\lambda _{01}}{\omega _1}\right)^2|\alpha |^2.$$
(61)
Equation (61) shows that the success of the measurement procedure is favored by large oscillator quality factors, large driven response $`\left|\alpha \right|`$, and a large value of the anharmonicity coupling factor $`\lambda _{01}/\omega _1`$. In addition, as we have seen, detecting individual jumps becomes harder as the temperature increases. Increasing the quality factors of mesoscopic oscillators is an active area of research. Currently, values of order $`10^3`$to $`10^4`$ seem possible. If these could be raised to the values characteristic of more macroscopic oscillators of the same material, of order $`10^6`$ or even higher, the detection of individual phonons would become correspondingly easier. The frequency ratio $`\omega _1/\omega _0`$ appearing in Eq. (61) must be less than unity for our detection scheme, but will probably not be too small because of geometry constraints. Thus the main parameters available to optimize the experimental geometry are the anharmonicity factor $`\lambda _{01}/\omega _1`$ and the dimensionless measure of the driven displacement of the ancilla, $`\left|\alpha \right|^2`$ (the number of phonons in the driven state). We now consider these factors in more detail.

### VI.1 Anharmonicity Coefficient

The interaction Hamiltonian for the system and ancilla oscillators Eq. (7) can be written

$$H^{\mathrm{RWA}}=\mathrm{}\omega _0a_0^{}a_0+\mathrm{}\left[\omega _1+\lambda _{01}n_0\right]a_1^{}a_1,$$
(62)
with $`n_0`$ the system phonon number. This equation implies that $`\lambda _{01}`$ can be estimated as the frequency shift of the ancilla oscillator for a single quantum ($`n_0=1`$) of the system oscillator.

For the prototype geometries using the two orthogonal flexing modes of a single beam, or parallel flexing modes of two longitudinally coupled beams, the nonlinear coupling arises from geometrical effects. At second order, the transverse displacement in one mode gives a longitudinal strain, which then changes the frequency of the second mode. The strain generated by the flexing motion and the frequency shift associated with this strain can be derived using elasticity theory, and have been calculated by Harrington and Roukes HR94 . A demonstration of such frequency shift detection, and direct measurement of $`\lambda _{01}`$ between two coupled beam has recently been reportedHPYR03 . The longitudinal strain produced by a single quantum in the fundamental flexing motion is

$$\chi \frac{\mathrm{}}{m_0\omega _0}\frac{1}{L_0^2},$$
(63)
where $`m_0`$ is the mass and $`L_0`$ is the length of the system beam. Then the ancilla frequency shift caused by this strain is

$$\lambda _{01}=\omega _1\frac{\zeta }{2\pi ^2}\chi \frac{L_1^2}{d_1^2},$$
(64)
where $`\zeta `$ is a geometric factor ($`\zeta =3`$ for clamped beam boundary conditions) and $`L_1,d_1`$ are the length and thickness of the ancilla beam, respectively. Introducing a dimensionless quantity,

$$R\frac{\mathrm{}^2}{m_1d_1^2}\frac{1}{\mathrm{}\omega _1},$$
(65)
then the scaled coupling coefficient can be expressed as

$$\frac{\lambda _{01}}{\omega _1}=\frac{\zeta }{2\pi ^2}\frac{m_1\omega _1L_1^2}{m_0\omega _0L_0^2}R.$$
(66)
Since the factor of the ratio of the two mode parameters will not be too large or small, the most important quantity determining the anharmonicity factor, which must not be too small for the success of the measurement scheme, is the dimensionless ratio $`R`$. This will typically be a small number. The need for small devices is seen from the scaling of this parameter with the dimensions.

### VI.2 Driving strength

The detection scheme we have considered is to measure the phase of the driven response of the ancilla oscillator. Since the detection scheme is magnetomotive, it is natural to consider the use of magnetic driving in estimating the size of the displacement parameter $`\left|\alpha \right|`$. For magnetic driving using a current $`I_{\mathrm{drive}}`$ in a magnetic field $`B`$ the dimensionless displacement can be estimated as

$$\left|\alpha \right|=Q_1\frac{BI_{\mathrm{drive}}L_1d_1}{\sqrt{2}\mathrm{}\omega _1}\sqrt{R}.$$
(67)
Again the important role of $`R`$ in limiting the size of $`\left|\alpha \right|`$ in this analysis is apparent.

We must also recognize that the size of $`\left|\alpha \right|`$ might be limited by other physical constraints, rather than by the available drive strength. One constraint might be to avoid undesired nonlinear effects in the driven beam itself. For a classical oscillator, at sufficiently large drive amplitudes nonlinear frequency pulling leads to a multiplicity of solutions and instability. This occurs when the nonlinear frequency shift is comparable with the width of the resonance $`\omega _1/Q_1`$. Using the same type of estimate for the nonlinear frequency shift as in Eq. (64) shows that this occurs for

$$\left|\alpha \right|\frac{1}{\sqrt{Q_1R}}.$$
(68)
A more detailed, quantum mechanical analysis of the driven nonlinear oscillator will be presented elsewhere qph04 .

### VI.3 Example Configuration

As a first estimate of the order of magnitude of the quantities introduced above we will construct an example configuration using parameters that seem plausible with current technology.

Recently, oscillators with resonant frequencies as high as $`1`$ GHz have been fabricated HZMR03  using silicon carbide. Thus we consider two flexing modes with resonant frequencies of $`\omega _0=2.3`$ GHz and $`\omega _1=0.36`$ GHz so that $`\omega _0\omega _1\lambda _{01},\nu ,\kappa `$ are satisfied. For this value of the system oscillator frequency, $`\mathrm{}\omega _0/k_\mathrm{B}T`$ is unity at a temperature of about $`0.1`$ K. The oscillators in Ref. HZMR03  were not very small, but it is expected that the structure can be scaled down while maintaining the high oscillation frequency. We therefore suppose smaller dimensions consistent with these frequencies, namely dimensions are $`0.6`$ $`\mu `$m $`\times `$ $`0.04`$ $`\mu `$m $`\times `$ $`0.07`$ $`\mu `$m for the system beam and $`0.6`$ $`\mu `$m $`\times `$ $`0.04`$ $`\mu `$m $`\times `$ $`0.01`$ $`\mu `$m for the ancilla beam. With these parameters we obtain

$$R=4.26\times 10^9.$$
(69)
The factor $`R`$ occurs squared in $`k/\nu `$ (via Eq. (66), which is required to be large, and so this small factor must be mitigated by the other quantities in Eq. (66), i.e. large values of Q and a large driven amplitude $`\left|\alpha \right|`$. Suppose the Q of the system oscillator is $`Q_0=10000`$ and $`Q_1=1000`$. For the size of $`\left|\alpha \right|`$ first consider the magnetic driving. A magnetic field $`B=10`$ Tesla and $`I_{\mathrm{drive}}=1`$ $`\mu `$A can raise the driven response to $`\left|\alpha \right|10^5`$. To reach $`k/\nu 1`$ the nonlinear coupling $`\lambda _{01}/\omega _1`$ that is required is then

$$\lambda _{01}/\omega _1=4.9\times 10^8.$$
(70)
With the given beam dimensions and the geometric nonlinearity, the anharmonic coupling coefficient is actually $`\lambda _{01}/\omega _1=1.3\times 10^{11}`$, about three orders of magnitude smaller. One possible way to increase this value might be to engineer the geometry of the oscillator so that the anharmonic coupling is larger than in the simple geometric nonlinearity we have consideredSP:R . Another way to increase $`k/\nu `$ is to increase $`\left|\alpha \right|`$ using a different driving scheme, although for the value or $`R`$ in Eq. (69) the limit in Eq. (68) is already exceeded for $`\left|\alpha \right|10^3`$ , so that engineering the geometry to reduce the self-nonlinearity might be necessary. An obvious way to increase $`k/\nu `$ to values greater than unity is to use oscillators with smaller dimensions, for example carbon nanotubes.

## VII Conclusion

We have analyzed a scheme to observe quantum transitions of a mesoscopic mechanical oscillator. The non-linear coupling shifts the frequency of a second (ancilla) oscillator proportionally to the excitation level of the first (system) oscillator. This frequency shift may be detected as a phase shift of the ancilla oscillation when driven on resonance. In principle, a QND measurement is possible if the coupling constant between the two oscillators $`\lambda _{01}`$ is much smaller than the resonance frequencies of the oscillators, as will usually be the case. We have derived the master equation for the system density matrix first integrating out the environment and measurement degrees of freedom, and then by removing the ancilla operator using the fact that the time scale of the system and ancilla dynamics are quite different. The master equation has three components: phase diffusion as a result of the measurement backaction; a constant energy shift due to the excitation of the ancilla oscillator, and number state transitions due to the interaction with the thermal bath (the environment).

The measurement process introduces a stochastic component into the system dynamics and we have obtained the stochastic master equation corresponding to our measurement scheme. From the stochastic master equation we identify two competing tendencies that can be characterized by two parameters. One is the coupling strength $`\nu `$ of the system and thermal bath, which is associated with the dwell time $`t_{\mathrm{dwell}}`$ between transitions. The other is the coefficient $`k`$, associated with measurements, which includes not only the coupling strength of the system to the measurement bath but also the anharmonic coupling strength between the oscillators, the driving amplitude. This coefficient is related to the measurement time $`t_\mathrm{m}`$ that is needed for a measurement to be able to produce an outcome with certainty. To observe clear quantum jumps we would need $`t_{\mathrm{dwell}}t_\mathrm{m}`$. If this condition is not satisfied, then the experimenter cannot infer the energy eigenstate of the system from the observed current.

Although our simple estimates based on plausible lithographically prepared oscillators yield values for the ratio $`t_{\mathrm{dwell}}/t_\mathrm{m}`$ too small for the observation of individual phonons, enhancements to the geometry and the trend to smaller device sizes should improve the outlook. The basic scheme and theoretical techniques developed here are fairly general, and in particular are not restricted to zero temperature, and so can be also used for other applications such as single spin detection and noise analysis for a solid state based quantum computer. Such possibilities might open up a new stage for observing quantum dynamics in mesoscopic systems.

###### Acknowledgements.

We thank Michael Roukes for providing information about the experimental progress in his group. DHS is grateful to Gerard Milburn and Tony Leggett for insightful discussions and thanks for the Institute for Quantum Informaion at Caltech for their hospitality. This work was supported by DARPA DSO/MOSAIC through grant N00014-02-1-0602 and NSF through grant DMR-9873573. DHS’s work is also supported by the NSF through a grant for the Institute for Theoretical Atomic, Molecular and Optical Physics at Harvard University and Smithsonian Astrophysical Observatory. DHS also acknowledges the Weitzman fund for a travel grant. Work by ACD was supported by the NSF through grant EIA-0086083 as part of the Institute for Quantum Information and the Caltech MURI Center for Quantum Networks (DAAD19-00-1-0374).

## Appendix A Bath field operators

In this appendix we describe more fully the time-local measurement bath operators $`B_t`$ introduced in section IV.1. The description in terms of finely spaced modes of the bath with a smooth density of states leads to the short memory or Markov property of the bath, that can be expressed in terms of the time-local commutation rules for $`B_t`$. In the main text we introduced the global bath operator as Eq. (26)

$$B_t=\frac{1}{\sqrt{2\pi \rho _\mathrm{d}(\omega _1)}g_\mathrm{d}\left(\omega _1\right)}\underset{n}{}g_\mathrm{d}\left(\omega _n\right)b_{\mathrm{d},n}e^{\iota \left(\omega _n\omega _1\right)t}.$$
(71)
We first derive the commutation rule Eq. (27). Substituting Eq. (71) in the commutator gives

$$[B\left(t\right),B^{}\left(t^{}\right)]=\frac{1}{2\pi \rho _\mathrm{d}\left(\omega _1\right)\left[g_\mathrm{d}\left(\omega _1\right)\right]^2}\underset{n,n^{}}{}g_\mathrm{d}\left(\omega _n\right)g_\mathrm{d}\left(\omega _n^{}\right)[b_{\mathrm{d},n},b_{\mathrm{d},n^{}}^{}]e^{i\left(\omega _n\omega _1\right)t}e^{i\left(\omega _n^{}\omega _1\right)t^{}}$$
(72)
and using $`[b_{\mathrm{d},n},b_{\mathrm{d},n^{}}^{}]=\delta _{n,n^{}}`$ we obtain

$$[B\left(t\right),B^{}\left(t^{}\right)]=\frac{1}{2\pi \rho _\mathrm{d}\left(\omega _1\right)\left[g_\mathrm{d}\left(\omega _1\right)\right]^2}\underset{n}{}\left[g_\mathrm{d}\left(\omega _n\right)\right]^2e^{i\left(\omega \omega _1\right)\left(tt^{}\right)}.$$
(73)
Changing the sum to integral form

$$\underset{n}{}_0^{\mathrm{}}𝑑\omega \rho _\mathrm{d}\left(\omega \right),$$
(74)
gives

$$[B\left(t\right),B^{}\left(t^{}\right)]=\frac{1}{2\pi \rho _\mathrm{d}\left(\omega _1\right)\left[g_d\left(\omega _1\right)\right]^2}_0^{\mathrm{}}𝑑\omega \rho _\mathrm{d}\left(\omega \right)\left[g_d\left(\omega \right)\right]^2e^{i\left(\omega \omega _1\right)\left(tt^{}\right)}.$$
(75)
Since $`\rho _\mathrm{d}\left(\omega \right)`$ and $`g_d(\omega )`$ are slowly varying functions around the ancilla oscillation frequency $`\omega =\omega _1`$ we can approximate these as$`\rho _\mathrm{d}\left(\omega \right)\rho _\mathrm{d}\left(\omega _1\right),`$ $`g_d(\omega )`$ $`g_d(\omega _1)`$ and pull them outside of the integral. Then introducing $`\epsilon =\omega \omega _1`$ and extending the lower range of the integration over $`\epsilon `$ to $`\mathrm{}`$ leads to the desired result

$$[B\left(t\right),B^{}\left(t^{}\right)]=\frac{1}{2\pi }_{\mathrm{}}^{\mathrm{}}𝑑\epsilon e^{i\epsilon \left(tt^{}\right)}=\delta \left(tt^{}\right).$$
(76)
The interaction Hamiltonian for the ancilla oscillator and the measurement bath is, from Eq. (12),

$$H_{\mathrm{int}}^\mathrm{I}=i\mathrm{}\underset{n}{}g_\mathrm{d}\left(\omega _n\right)\left[b_{\mathrm{d},n}^{}(t)a_1(t)b_{\mathrm{d},n}(t)a_1^{}(t)\right],$$
(77)
where we have moved to the interaction picture with $`a_1(t)=a_1e^{i\omega _1t}`$ and $`b_{\mathrm{d},n}(t)=b_{\mathrm{d},n}e^{i\omega _nt}`$ the ancilla and bath operators in this picture, and $`g_\mathrm{d}\left(\omega _n\right)`$ is the coupling strength. The interaction Hamiltonian can be written in terms of the bath operators $`B_t`$ as

$$H_{\mathrm{int}}^\mathrm{I}(t)=i\mathrm{}\sqrt{2\mu }\left(B_t^{}a_1B_ta_1^{}\right),$$
(78)
where the coefficient $`\mu `$ is

$$\mu \pi \varrho _\mathrm{d}\left(\omega _1\right)\left|g_\mathrm{d}\left(\omega _1\right)\right|^2$$
(79)
as before, and we have used the fact that the ancilla interacts predominantly with bath modes near frequency $`\omega _1`$ and again have assumed a smooth variation of the density of states and coupling constant so that we can make the replacements $`\rho _\mathrm{d}\left(\omega \right)\rho _\mathrm{d}\left(\omega _1\right)`$, and $`g_\mathrm{d}(\omega )`$ $`g_\mathrm{d}(\omega _1)`$.