text:
# Realistic simulations of single-spin nondemolition measurement by magnetic resonance force microscopy

## 1 Introduction

Single-spin measurement is an extremely important challenge, and necessary for the future successful development of several recent spin-based proposals for quantum information processing.  There are both direct and indirect single-spin measurement proposals. The idea behind some indirect proposals is to transform the problem of detecting a single spin into the task of measuring charge transport , since the ability to detect a single charge is now available. For direct single-spin detection, magnetic resonance force microscopy (MRFM) has been suggested  as one of the most promising techniques. To date, the MRFM technique has been demonstrated with sensitivity to a few hundred spins .

In this paper we discuss how to read out the quantum state of a single spin using the MRFM technique based on cyclic adiabatic inversion (CAI).  In this CAI MRFM technique, the frequency of the spin inversion in the rotating frame is in resonance with the mechanical vibration of an ultra thin cantilever, allowing it to amplify the otherwise extremely weak force due to the spin. These amplified vibrations can then be detected by, e.g., optical methods.

Previous studies  of the dynamics of single-spin measurement by MRFM considered only the unitary evolution of the spin and the cantilever system, without including any effects of external environments or measurement devices. Only recently, the effect of thermal noise environment on the dynamics of the spin-cantilever system in the MRFM was studied  by using the Caldeira-Leggett master equation  in the high-temperature limit.

There is, however, a macroscopic device in the MRFM setup which measures the cantilever motion and hence provides information about the spin state. To our knowledge, the back-action of the measurement device and the effect of the thermal noise on the dynamics of the cantilever-spin system for the single-spin detection problem by MRFM have not yet been investigated systematically. In this paper, we include, in our analysis, a measurement device (a fiber-optic interferometer) to monitor the position of the cantilever. We consider various relevant sources of noise and calculate the signal-to-noise ratio of the output photocurrent of the measurement device. We also develop a realistic continuous measurement model, and discuss the approximations and conditions to achieve a quantum non-demolition measurement of a single spin by MRFM. Finally, we present some simulation results of the dynamics of the single-spin measurement process.

## 2 The measurement scheme

A schematic illustration of the MRFM setup is shown in Fig 1. A uniform magnetic field, $`B_0`$, points in the positive $`z`$-direction. A single spin is placed in front of the cantilever tip which can oscillate only in the $`z`$-direction. A ferromagnetic particle (or small magnetic material) mounted on the cantilever tip produces a non-uniform magnetic field or magnetic field gradient of $`(B_z/Z)_0`$ on the single spin. As a result, a reactive force (or interaction) acts back on the magnetic cantilever tip in the $`z`$-direction from the single spin. The origin is chosen to be the equilibrium position of the cantilever tip without the presence of the spin.

In CAI, the cantilever is driven at its resonance frequency to amplify the otherwise very small vibrational amplitude. This is achieved by a modulation scheme using the frequency modulation of a rotating radio-frequency (RF) magnetic field in the $`x`$-$`y`$ plane. In this case, the rotating RF field can be represented as $`B_{1x}=B_1\mathrm{cos}\{[\omega +\mathrm{\Delta }\omega (t)]t\}`$, $`B_{1y}=B_1\mathrm{sin}\{[\omega +\mathrm{\Delta }\omega (t)]t\}`$, where the frequency modulation $`\mathrm{\Delta }\omega (t)`$ is a periodic function in time with the resonant frequency $`\omega _m`$ of the cantilever. In the reference frame rotating with the $`𝐁_1`$, the spin-cantilever Hamiltonian can be written as

$$\widehat{H}_{SZ}(t)=\widehat{H}_Z\mathrm{}\left[\omega _L\omega \mathrm{\Delta }\omega (t)\right]\widehat{S}_z\mathrm{}\omega _1\widehat{S}_xg\mu \left(\frac{B_z}{Z}\right)_0\widehat{Z}\widehat{S}_z,$$
(1)
where $`\omega _L=g\mu B_z/\mathrm{}`$ and $`\omega _1=g\mu B_1/\mathrm{}`$ are respectively the Larmor and Rabi frequencies, $`B_z`$ includes the uniform magnetic field $`B_0`$ and the magnetic field produced by the ferromagnetic particle, $`g`$ and $`\mu `$ are the $`g`$-factor and the electron or nuclear magneton, and

$$\widehat{H}_Z=\frac{1}{2m}\widehat{p}^2+\frac{m\omega _m^2}{2}\widehat{Z}^2$$
(2)
is the Hamiltonian of the cantilever in isolation (i.e., with no external magnetic field coupling it to the spin). For $`\omega =\omega _L`$, we arrive at an effective cantilever-spin Hamiltonian of the form

$$\widehat{H}_{SZ}(t)=\widehat{H}_Z2\eta \widehat{Z}\widehat{S}_z+f(t)\widehat{S}_z\epsilon \widehat{S}_x,$$
(3)
where $`f(t)=\mathrm{\Delta }\omega (t)`$, $`\eta =(g\mu /2)(B_z/Z)_0`$ and $`\epsilon =\mathrm{}\omega _1`$. We will discuss in details the rotating picture and adiabatic approximation for the spin-cantilever system in the next section.

In the following, we briefly describe the basic principle of the single-spin measurement by CAI MRFM. In the case when the adiabatic approximation is exact, the instantaneous eigenstates of the spin Hamiltonian in the rotating frame of the $`𝐁_1`$ field are the spin states parallel or antiparallel to the direction of the effective magnetic field $`𝐁^{\mathrm{eff}}(t)=(\epsilon ,0,f(t))`$, denoted as $`|v_\pm (t)`$, respectively. We define an operator $`\widehat{S}_z^{}`$ for the component of spin along this axis. Note that the initial spin state in the laboratory frame has the same expression as the initial state in the rotating frame. Starting at a general initial spin state (in the laboratory or rotating frame) of

$$\chi (0)=a|+b|$$
(4)
in the $`\widehat{S}_z`$ representation, we can rewrite this initial state in the basis of the instantaneous eigenstates of $`\widehat{S}_z^{}`$ as

$$\chi (0)=a_{\mathrm{eff}}|v_+(0)+b_{\mathrm{eff}}|v_{}(0),$$
(5)
where

$`a_{\mathrm{eff}}`$ $`=`$ $`a\mathrm{cos}(\mathrm{\Theta }_0/2)+b\mathrm{sin}(\mathrm{\Theta }_0/2),`$ (6)
$`b_{\mathrm{eff}}`$ $`=`$ $`a\mathrm{sin}(\mathrm{\Theta }_0/2)+b\mathrm{cos}(\mathrm{\Theta }_0/2),`$ (7)
and $`\mathrm{\Theta }_0\mathrm{\Theta }(0)`$ is the initial angle between the effective magnetic magnetic field and the $`z`$-axis direction. This implies $`\mathrm{tan}[\mathrm{\Theta }(t)]=B_x^{\mathrm{eff}}(t)/B_z^{\mathrm{eff}}(t)=\epsilon /f(t)`$. It then follows from the adiabatic theorem that the spin state at time $`t`$ can be written as:

$`\chi (t)`$ $`=`$ $`a_{\mathrm{eff}}|v_+(t)\mathrm{exp}({\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _0^t}\lambda _+(t^{})𝑑t^{})`$ (8)
$`+b_{\mathrm{eff}}|v_{}(t)\mathrm{exp}({\displaystyle \frac{i}{\mathrm{}}}{\displaystyle _0^t}\lambda _{}(t^{})𝑑t^{}),`$
where $`\lambda _\pm (t)`$ are instantaneous eigenvalues. So the probabilities of finding the spin to be in the instantaneous eigenstates $`|v_\pm (t)`$ are respectively $`|a_{\mathrm{eff}}|^2`$ and $`|b_{\mathrm{eff}}|^2`$. Since the coefficients $`a_{\mathrm{eff}}`$ and $`b_{\mathrm{eff}}`$ are time independent, the probabilities $`|a_{\mathrm{eff}}|^2`$ and $`|b_{\mathrm{eff}}|^2`$ remain the same at all times. This provides us with an opportunity to measure the initial spin state probabilities at later times.

How do we measure these spin state probabilities? The idea is to transfer the information of the spin state to the state of the driven cantilever. In the interaction picture in which the state is rotating with the instantaneous eigenstates of the spin Hamiltonian, the spin-cantilever interaction can be written as $`2\eta \widehat{Z}\widehat{S}_z^{}\mathrm{cos}[\mathrm{\Theta }(t)]`$. As a result, the phase of the driven cantilever vibrations depends on the orientation of the spin states. Suppose that the initial state is a product state of the cantilever and spin parts. At a later time, due to the interaction between them, the total state becomes entangled. Monitoring the phase of the cantilever vibrations will give us the information about the spin. Numerical simulations (see Fig. 4 with reasonable parameters for the CAI approximations) indicate that as the amplitude of the cantilever vibrations increases with time, the phase difference in the oscillations for the two different initial spin eigenstates of $`\widehat{S}_z^{}`$ approaches $`\pi `$. In other words, the measurement of the single-spin states can be achieved by monitoring the phases of the cantilever vibrations at some later time $`t`$. Phase-sensitive, optical homodyne measurements of the cantilever vibrations can be performed using a fiber-optic interferometer. The main purpose of this paper is to present a realistic and detailed analysis of the single-spin measurement scheme, including the effects of the measurement device and other relevant sources of noise.

## 3 The rotating picture and adiabatic approximation

We assume an effective cantilever-spin Hamiltonian of the form (3) where for the moment we let $`f(t)`$ and $`\epsilon `$ be arbitrary, and $`\widehat{H}_Z`$ is the Hamiltonian given by (2). It is useful to group this into three terms

$$\widehat{H}_{SZ}(t)=\widehat{H}_Z+\widehat{H}_I+\widehat{H}_S(t),$$
(9)
where

$`\widehat{H}_I`$ $``$ $`2\eta \widehat{Z}\widehat{S}_z,`$
$`\widehat{H}_S(t)`$ $``$ $`f(t)\widehat{S}_z\epsilon \widehat{S}_x.`$ (10)
The state of the cantilever-spin system evolves according to the Schrödinger equation

$$\frac{d|\psi (t)}{dt}=\frac{i}{\mathrm{}}\widehat{H}_{SZ}(t)|\psi (t).$$
(11)
In realistic cases, the spin part of the Hamiltonian (representing precession under the magnetic field) gives an evolution which is very rapid compared to the reaction time of the cantilever. It therefore makes sense to switch to an interaction picture in which the state is rotating along with this precession. We do this by introducing a (partial) time translation operator

$$\widehat{U}_S(t):\mathrm{exp}\frac{i}{\mathrm{}}\left[_0^t\widehat{H}_S(t^{})dt^{}\right]:,$$
(12)
where $`::`$ indicates that the integral is to be taken in a time-ordered sense; this unitary operator obeys the differential equation

$$\frac{d\widehat{U}_S(t)}{dt}=\frac{i}{\mathrm{}}\widehat{H}_S(t)\widehat{U}_S(t).$$
(13)
We then introduce the state $`|\stackrel{~}{\psi }(t)`$ in the rotating picture:

$$|\stackrel{~}{\psi }(t)\widehat{U}_S^{}(t)|\psi (t),$$
(14)
with $`|\psi (t)`$ the solution of the original Schrödinger equation (11) at time $`t`$. The evolution equation for $`|\stackrel{~}{\psi }(t)`$ is

$`{\displaystyle \frac{d|\stackrel{~}{\psi }(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{d\widehat{U}_S^{}(t)}{dt}}|\psi (t)+\widehat{U}_S^{}(t){\displaystyle \frac{d|\psi (t)}{dt}}`$ (15)
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{U}_S^{}(t)\widehat{H}_S(t)|\psi (t){\displaystyle \frac{i}{\mathrm{}}}\widehat{U}_S^{}(t)\widehat{H}_{SZ}(t)|\psi (t)`$
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_Z|\stackrel{~}{\psi }(t){\displaystyle \frac{i}{\mathrm{}}}\left[\widehat{U}_S^{}(t)\widehat{H}_I\widehat{U}_S(t)\right]|\stackrel{~}{\psi }(t)`$
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_Z|\stackrel{~}{\psi }(t)+{\displaystyle \frac{2i\eta }{\mathrm{}}}\widehat{Z}\left[\widehat{U}_S^{}(t)\widehat{S}_z\widehat{U}_S(t)\right]|\stackrel{~}{\psi }(t).`$
We can define a locked spin operator $`\widehat{S}_L(t)`$

$$\widehat{S}_L(t)\left[\widehat{U}_S^{}(t)\widehat{S}_z\widehat{U}_S(t)\right];$$
(16)
in terms of this, the equation of motion for $`|\stackrel{~}{\psi }`$ becomes

$$\frac{d|\stackrel{~}{\psi }(t)}{dt}=\frac{i}{\mathrm{}}\widehat{H}_Z|\stackrel{~}{\psi }(t)+\frac{2i\eta }{\mathrm{}}\widehat{Z}\widehat{S}_L(t)|\stackrel{~}{\psi }(t).$$
(17)
Unfortunately, it is difficult to get an exact solution for $`\widehat{U}_S(t)`$ for a general function $`f(t)`$. This means that it is also difficult to derive an exact expression for $`\widehat{S}_L(t)`$, and the rotating picture (15), while formally correct, is not very helpful.

However, while we cannot easily find an exact expression for $`\widehat{U}_S(t)`$ for general $`f(t)`$, we can easily find an approximate solution for a large class of functions. Suppose that $`\epsilon `$ is large and $`f(t)`$ is slowly varying, so that $`|f(t)|,\epsilon |f^{}(t)/f(t)|`$ for typical values of $`f(t)`$ and $`f^{}(t)`$. Then $`\widehat{H}_S(t)`$ is also slowly varying, and if a spin begins in an instantaneous eigenstate of $`\widehat{H}_S(t)`$, it will remain close to an instantaneous eigenstate of $`\widehat{H}_S(t)`$ for all times by the adiabatic theorem.

The instantaneous eigenstates of $`\widehat{H}_S(t)`$ are

$$\widehat{H}_S(t)|v_\pm (t)=\lambda _\pm (t)|v_\pm (t)\pm \lambda (t)|v_\pm (t),$$
(18)
where

$`\lambda (t)`$ $`=`$ $`\sqrt{f^2(t)+\epsilon ^2},`$
$`|v_\pm (t)`$ $`=`$ $`{\displaystyle \frac{\epsilon }{\sqrt{(f(t)\lambda (t))^2+\epsilon ^2}}}|{\displaystyle \frac{f(t)\lambda (t)}{\sqrt{(f(t)\lambda (t))^2+\epsilon ^2}}}|.`$ (19)
We use these instantaneous eigenvectors and eigenvalues to define an approximation to the unitary operator $`\widehat{U}_S(t)`$:

$$\widehat{U}_S^{}(t)=\widehat{I}|v_+(t)v_+(0)|\mathrm{e}^{i\mathrm{\Phi }(t)}+\widehat{I}|v_{}(t)v_{}(0)|\mathrm{e}^{i\mathrm{\Phi }(t)},$$
(20)
with the accumulated phase

$$\mathrm{\Phi }(t)\frac{1}{\mathrm{}}_0^t\lambda (t^{})𝑑t^{}.$$
(21)
Note that $`\mathrm{\Phi }(t)`$ obeys $`d\mathrm{\Phi }(t)/dt=\lambda (t)`$. This implies that

$`{\displaystyle \frac{d\widehat{U}_S^{}(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\lambda (t)\widehat{I}|v_+(t)v_+(0)|\mathrm{e}^{i\mathrm{\Phi }(t)}+{\displaystyle \frac{i}{\mathrm{}}}\lambda (t)\widehat{I}|v_{}(t)v_{}(0)|\mathrm{e}^{i\mathrm{\Phi }(t)}`$ (22)
$`+\widehat{I}{\displaystyle \frac{d|v_+(t)}{dt}}v_+(0)|\mathrm{e}^{i\mathrm{\Phi }(t)}+\widehat{I}{\displaystyle \frac{d|v_{}(t)}{dt}}v_{}(0)|\mathrm{e}^{i\mathrm{\Phi }(t)}`$
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_S(t)\widehat{U}_S^{}(t)+\widehat{I}{\displaystyle \frac{d|v_+(t)}{dt}}v_+(0)|\mathrm{e}^{i\mathrm{\Phi }(t)}`$
$`+\widehat{I}{\displaystyle \frac{d|v_{}(t)}{dt}}v_{}(0)|\mathrm{e}^{i\mathrm{\Phi }(t)},`$
which has the form of (13) plus some additional terms. From the definition (19) of $`|v_\pm (t)`$, we see

$$\frac{d|v_\pm (t)}{dt}=\pm \frac{1}{2}\frac{\epsilon }{\lambda ^2(t)}\frac{df(t)}{dt}|v_{}(t).$$
(23)
Provided that $`f(t)`$ is slowly varying, the additional terms in (22) will be small.

Just as before, we can define a rotating picture, now using the unitary transformation $`\widehat{U}_S^{}(t)`$,

$$|\stackrel{˘}{\psi }(t)\left(\widehat{U}_S^{}(t)\right)^{}|\psi (t).$$
(24)
This gives us a new evolution equation for $`|\stackrel{˘}{\psi }`$:

$`{\displaystyle \frac{d|\stackrel{˘}{\psi }(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{d(\widehat{U}_S^{}(t))^{}}{dt}}|\psi (t)+(\widehat{U}_S^{}(t))^{}{\displaystyle \frac{d|\psi (t)}{dt}}`$ (25)
$`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_Z|\stackrel{˘}{\psi }(t)+{\displaystyle \frac{2i\eta }{\mathrm{}}}\widehat{Z}\left[(\widehat{U}_S^{}(t))^{}\widehat{S}_z\widehat{U}_S^{}(t)\right]|\stackrel{˘}{\psi }(t)`$
$`+\widehat{I}(|v_+(0){\displaystyle \frac{dv_+(t)|}{dt}}\mathrm{e}^{i\mathrm{\Phi }(t)}`$
$`+|v_{}(0){\displaystyle \frac{dv_{}(t)|}{dt}}\mathrm{e}^{i\mathrm{\Phi }(t)})\widehat{U}_S^{}(t)|\stackrel{˘}{\psi }(t).`$
At this point, it is helpful to introduce a new set of spin operators

$`\widehat{S}_x^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2}}\widehat{I}\left(|v_+(0)v_{}(0)|+|v_{}(0)v_+(0)|\right),`$
$`\widehat{S}_y^{}`$ $`=`$ $`{\displaystyle \frac{i\mathrm{}}{2}}\widehat{I}\left(|v_{}(0)v_+(0)||v_+(0)v_{}(0)|\right),`$
$`\widehat{S}_z^{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2}}\widehat{I}\left(|v_+(0)v_+(0)||v_{}(0)v_{}(0)|\right).`$ (26)
Using the definition (20) for $`\widehat{U}_S^{}(t)`$, we can solve for the various terms in (25):

$$(\widehat{U}_S^{}(t))^{}\widehat{S}_z\widehat{U}_S^{}(t)=\frac{f(t)}{\lambda (t)}\widehat{S}_z^{}\frac{\epsilon }{\lambda (t)}\left(\widehat{S}_x^{}\mathrm{cos}(2\mathrm{\Phi }(t))\widehat{S}_y^{}\mathrm{sin}(2\mathrm{\Phi }(t))\right).$$
(27)
$`\widehat{I}\left(|v_+(0){\displaystyle \frac{dv_+(t)|}{dt}}\mathrm{e}^{i\mathrm{\Phi }(t)}+|v_{}(0){\displaystyle \frac{dv_{}(t)|}{dt}}\mathrm{e}^{i\mathrm{\Phi }(t)}\right)\widehat{U}_S^{}(t)=`$
$`{\displaystyle \frac{i\epsilon }{\mathrm{}\lambda ^2(t)}}{\displaystyle \frac{df(t)}{dt}}\left(\widehat{S}_x^{}\mathrm{sin}(2\mathrm{\Phi }(t))+\widehat{S}_y^{}\mathrm{cos}(2\mathrm{\Phi }(t))\right).`$ (28)
Substituting (2628) into (25), we get

$`{\displaystyle \frac{d|\stackrel{˘}{\psi }(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_Z|\stackrel{˘}{\psi }(t)+{\displaystyle \frac{2i\eta }{\mathrm{}}}\widehat{Z}{\displaystyle \frac{f(t)}{\lambda (t)}}\widehat{S}_z^{}|\stackrel{˘}{\psi }(t)`$ (29)
$`+{\displaystyle \frac{2i\eta }{\mathrm{}}}\widehat{Z}{\displaystyle \frac{\epsilon }{\lambda (t)}}\left(\widehat{S}_x^{}\mathrm{cos}(2\mathrm{\Phi }(t))\widehat{S}_y^{}\mathrm{sin}(2\mathrm{\Phi }(t))\right)|\stackrel{˘}{\psi }(t)`$
$`+i{\displaystyle \frac{\epsilon }{\mathrm{}\lambda ^2(t)}}{\displaystyle \frac{df(t)}{dt}}\left(\widehat{S}_x^{}\mathrm{sin}(2\mathrm{\Phi }(t))+\widehat{S}_y^{}\mathrm{cos}(2\mathrm{\Phi }(t))\right)|\stackrel{˘}{\psi }(t).`$
Note that this equation is still exact—it is equivalent to the original Schrödinger equation (11). However, we can see that if $`|f(t)|,\epsilon `$ are large, then $`\mathrm{\Phi }(t)`$ will be a rapidly growing function, and the last two terms of equation (29) will oscillate very rapidly compared to the first two terms. Over a short period relative to the response time of the cantilever they will essentially average away to nothing. In this limit, therefore, we can reasonably make a rotating-wave approximation, to get the approximate evolution equation

$$\frac{d|\stackrel{˘}{\psi }(t)}{dt}\frac{i}{\mathrm{}}\left(\widehat{H}_Z2\eta (f(t)/\lambda (t))\widehat{Z}\widehat{S}_z^{}\right)|\stackrel{˘}{\psi }(t).$$
(30)
This is equivalent to making an exact adiabatic approximation, as described in section 2. We can see how this approximation compares to the complete Hamiltonian for a reasonable set of parameter values in figures 2 and 3. This set of parameters was chosen to match those of Berman et al. , as was the set of times plotted in figure 3. Comparison shows that our results match their unitary simulations to good precision. For the rest of this paper we will be using the rotating wave approximation, and representing states in the rotating frame. For simplicity, we henceforth omit the accent from the state $`|\stackrel{˘}{\psi }`$.

In this rotating-wave approximation, if the spin begins in an instantaneous eigenstate of $`\widehat{H}_S(t)`$, it will remain in an instantaneous eigenstate at all times. If it begins in a superposition of the two eigenstates, the spin and cantilever degrees of freedom will become entangled, with the two components of the wavefunction corresponding to the two spin directions remaining undisturbed for all times. Monitoring the position of the cantilever then serves as a nondemolition measurement of the spin, as we would wish.

Note that the corrections to the adiabatic approximation include terms which can flip the spin. These terms must remain small for the system to be a true nondemolition measurement. The result of the spin measurement manifests itself as a $`\pi `$ phase shift in the oscillation of the cantilever. We can see this in figure 4.

## 4 The thermal environment

Unfortunately, in practice we cannot treat the cantilever as an isolated system. It is coupled at least weakly to the vibrational modes of the bulk, and is therefore subject to dissipation and thermal noise. Since the cantilever can be treated as a single harmonic oscillator, we can model the effects of this thermal bath by the well-known Caldeira-Leggett  master equation in the high-temperature limit:

$$\dot{\rho }=\frac{i}{\mathrm{}}[\widehat{H}_{SZ}(t),\rho ]\frac{i\gamma _m}{\mathrm{}}[\widehat{Z},\{\widehat{p},\rho \}]\frac{\gamma _m}{2\mathrm{}^2}[\widehat{Z},[\widehat{Z},\rho ]],$$
(31)
where the parameters are

$`\gamma _m`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }}{2m}},`$
$`\mathrm{}`$ $`=`$ $`{\displaystyle \frac{\mathrm{}}{2\sqrt{mkT}}},`$ (32)
$`m`$ is the cantilever mass, $`T`$ is the temperature, $`k`$ is Boltzmann’s constant (or the equivalent for our system of units), and $`\mathrm{\Gamma }`$ is the strength of the coupling to the thermal bath. We can interpret $`\gamma _m`$ (with units of inverse time) as the dissipation rate and $`\mathrm{}`$ (with units of length) as the thermal de Broglie wavelength.

A feature of this equation is that it doesn’t necessarily preserve the positivity of $`\rho `$ on short time scales (though at long times it is well-behaved) . This arises because of the approximations which are made in the derivation, which become invalid at very short times. While this may be physically unimportant, it can be inconvenient; in particular, if we wish to unravel the evolution into a stochastic Schrödinger equation  (as we will in section 6), it is necessary to start with a master equation in Lindblad form,

$$\dot{\rho }=\frac{i}{\mathrm{}}[\widehat{H},\rho ]+\underset{j}{}\left[2\widehat{L}_j\rho \widehat{L}_j^{}\{\widehat{L}_j^{}\widehat{L}_j,\rho \}\right],$$
(33)
for some Hermitian $`\widehat{H}`$ and set of general Lindblad operators $`\{\widehat{L}_j\}`$. The Caldeira-Leggett equation (31) is not of this form, which is why it can violate positivity of $`\rho `$.

The exact quantum Brownian motion master equation was shown  not to have the Lindblad form, but rather requires time-dependent coefficients to ensure the positivity of the density matrix at short times. However, by keeping more terms from the high- or medium-temperature-limit expansion in a consistent way, Diósi  has shown that the Caldeira-Leggett equation can be replaced by another master equation which is of Lindblad form, and which agrees with it except at very short times when the equation’s validity is questionable in any case. This is done by adding a term to (31) of the form $`(\gamma _m\mathrm{}^2/2\mathrm{}^2)[\widehat{p},[\widehat{p},\rho ]]`$. The procedure is analogous to completing the square. If we choose the ansatz

$$\widehat{L}=A\widehat{Z}+iB\widehat{p}$$
(34)
with real $`A,B`$, plug it into the equation (33), and equate it to the Caldeira-Leggett equation (31) plus the additional term, we get

$`\dot{\rho }`$ $`=`$ $`(i/\mathrm{})[\widehat{H},\rho ]A^2[\widehat{Z},[\widehat{Z},\rho ]]B^2[\widehat{p},[\widehat{p},\rho ]]`$ (35)
$`+iAB(2\widehat{Z}\rho \widehat{p}+\widehat{Z}\widehat{p}\rho +\rho \widehat{Z}\widehat{p}+2\widehat{p}\rho \widehat{Z}\widehat{p}\widehat{Z}\rho \rho \widehat{p}\widehat{Z})`$
$`=`$ $`(i/\mathrm{})[\widehat{H}_{SZ}(t),\rho ]{\displaystyle \frac{\gamma _m}{2\mathrm{}^2}}[\widehat{Z},[\widehat{Z},\rho ]]{\displaystyle \frac{\gamma _m\mathrm{}^2}{2\mathrm{}^2}}[\widehat{p},[\widehat{p},\rho ]]`$
$`+{\displaystyle \frac{i\gamma _m}{\mathrm{}}}(\widehat{p}\rho \widehat{Z}\widehat{Z}\widehat{p}\rho \widehat{Z}\rho \widehat{p}+\rho \widehat{p}\widehat{Z}),`$
which implies that

$`A`$ $`=`$ $`\sqrt{\gamma _m/2\mathrm{}^2},`$
$`B`$ $`=`$ $`\sqrt{\gamma _m\mathrm{}^2/2\mathrm{}^2},`$
$`\widehat{H}`$ $`=`$ $`\widehat{H}_{SZ}(t)+(\gamma _m/2)(\widehat{Z}\widehat{p}+\widehat{p}\widehat{Z})\widehat{H}_{SZ}^{}(t).`$ (36)
So the Lindblad operator for this equation is

$$\widehat{L}=\sqrt{\gamma _m/2}\left((1/\mathrm{})\widehat{Z}+i(\mathrm{}/\mathrm{})\widehat{p}\right),$$
(37)
and the effective Hamiltonian, going to the rotating picture and making use of the approximation derived in section 3, is

$$\widehat{H}_{SZ}^{}(t)=\frac{1}{2m}\widehat{p}^2+\frac{m\omega _m^2}{2}\widehat{Z}^22\eta (f(t)/\lambda (t))\widehat{Z}\widehat{S}_z^{}+(\gamma _m/2)(\widehat{Z}\widehat{p}+\widehat{p}\widehat{Z}).$$
(38)
In order for the cantilever to be an effective measurement device, the loss rate must be very low: $`\omega _m\gamma _m`$.

## 5 The effects of monitoring

In order to serve as a measurement scheme, we must have some way of monitoring the motion of the cantilever. Because of the microscopic scale of the motion, this is not so easily done. One approach is to use optical interferometry to measure the cantilever position.

As shown in figure 1, the cantilever forms one side of an optical microcavity and the cleaved end of the fiber forms the other side. As the cantilever moves, the resonant frequency of the cavity changes. Because the timescale of the cantilever’s motion is very long compared to the optical timescale, we can treat the effects of this in the adiabatic limit. The cavity mode is also subject to driving by an external laser, and has a very high loss rate. The full master equation  for the cantilever-spin-cavity system in the interaction picture is

$`\dot{\rho }`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho ]+2\widehat{L}\rho \widehat{L}^{}\widehat{L}^{}\widehat{L}\rho \rho \widehat{L}^{}\widehat{L}`$ (39)
$`i[E(\widehat{a}^{}+\widehat{a})+\widehat{a}^{}\widehat{a}(\mathrm{\Delta }+\kappa \widehat{Z}),\rho ]`$
$`+(\gamma _c/2)(2\widehat{a}\rho \widehat{a}^{}\widehat{a}^{}\widehat{a}\rho \rho \widehat{a}^{}\widehat{a}),`$
where $`\widehat{H}_{SZ}^{}(t)`$ and $`\widehat{L}`$ are the Hamiltonian and Lindblad operator for the cantilever and spin given by eqns. (37) and (38), $`E`$ is the strength of the laser driving, $`\mathrm{\Delta }`$ is the detuning from the “neutral” cavity frequency, $`\kappa `$ is the coupling strength of the cantilever to the cavity mode, and $`\gamma _c`$ is the loss rate of the cavity.

Suppose now that we perform homodyne measurement  on the light which escapes from the cavity. We would like to replace the equation (39) above with an equation for the conditional evolution of $`\rho `$, conditioned on the output photocurrent $`I_c(t)`$. The conditional evolution equation for our system then becomes  (in Itô calculus form)

$`d\rho `$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho ]dt+\left(2\widehat{L}\rho \widehat{L}^{}\widehat{L}^{}\widehat{L}\rho \rho \widehat{L}^{}\widehat{L}\right)dt`$ (40)
$`i[E(\widehat{a}^{}+\widehat{a})+\widehat{a}^{}\widehat{a}(\mathrm{\Delta }+\kappa \widehat{Z}),\rho ]dt`$
$`+(\gamma _c/2)\left(2\widehat{a}\rho \widehat{a}^{}\widehat{a}^{}\widehat{a}\rho \rho \widehat{a}^{}\widehat{a}\right)dt`$
$`+\sqrt{\gamma _ce_d}\left(\widehat{a}\rho +\rho \widehat{a}^{}\widehat{a}+\widehat{a}^{}\rho \right)dW_t,`$
where $`0e_d1`$ is the detector efficiency and $`dW_t`$ is a real stochastic differential variable which obeys the statistics

$$M[dW_t]=0,M[dW_tdW_s]=\delta (ts)dt.$$
(41)
This noise is related to the output photocurrent

$$I_c(t)=\beta \left[\gamma _ce_d\widehat{a}+\widehat{a}^{}_t+\sqrt{\gamma _ce_d}\frac{dW_t}{dt}\right],$$
(42)
where $`\beta `$ is a constant giving the device’s range of response.

We want to operate in the “bad cavity” limit where $`\gamma _c\omega _m`$. This means that the cavity mode will approach equilibrium on a timescale very short compared to that of the cantilever’s motion, so that the cavity mode can be adiabatically eliminated  from this equation, leaving an equation in terms of the spin and cantilever position alone.

Let the detuning vanish $`\mathrm{\Delta }0`$ and the coupling $`\kappa `$ to the cantilever be very small. If we initially neglect this coupling altogether, we can solve for the steady-state of the cavity mode in isolation from the cantilever:

$`i[E(\widehat{a}^{}+\widehat{a}),\rho ]+(\gamma _c/2)\left(2\widehat{a}\rho \widehat{a}^{}\widehat{a}^{}\widehat{a}\rho \rho \widehat{a}^{}\widehat{a}\right)`$ $`=`$ $`0,`$
$`\left(\widehat{a}\rho +\rho \widehat{a}^{}\widehat{a}+\widehat{a}^{}\rho \right)`$ $`=`$ $`0,`$ (43)
which implies that $`\rho =|\alpha _0\alpha _0|`$, where $`\widehat{a}|\alpha _0=\alpha _0|\alpha _0`$ is a coherent state with

$$\alpha _0=\frac{2iE}{\gamma _c}.$$
(44)
Now let us restore the coupling $`\kappa `$ between the cantilever and the cavity mode. If this coupling is very small, then the state of the cavity mode will remain very close to the state $`|\alpha _0`$. In this case, it is very useful to switch to a displaced basis  for the cavity mode. We switch from the operators $`\widehat{a},\widehat{a}^{}`$ to displaced operators

$`\widehat{b}\widehat{a}\alpha _0,`$
$`\widehat{b}^{}\widehat{a}^{}\alpha _0^{},`$ (45)
and displaced number states

$$\widehat{b}^{}\widehat{b}|n=n|n.$$
(46)
Obviously $`|0=|\alpha _0`$, and $`|1=\widehat{a}^{}|\alpha _0\alpha _0^{}|\alpha _0`$.

We now make the ansatz of keeping the two lowest displaced number states $`|0,1`$ of the cavity mode and neglecting the rest.  We then write the full density matrix for the spin-cantilever-cavity system as

$$\rho (t)=\rho _0(t)|00|+\rho _1(t)|10|+\rho _1^{}(t)|01|+\rho _2(t)|11|,$$
(47)
where $`\rho _{0,1,2}`$ are operators which act on the Hilbert space of the cantilever and spin, and $`\rho _{0,2}`$ are self-adjoint. The reduced density matrix of the spin-cantilever alone is obtained by tracing out the cavity mode, yielding

$$\rho _{SZ}(t)=\rho _0(t)+\rho _2(t).$$
(48)
If we substitute the definitions (45) and (47) into the stochastic master equation (40) and collect terms, we get a set of coupled equations in the operators $`\rho _{0,1,2}`$:

$`d\rho _0`$ $`=`$ $`\left({\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho _0]+2\widehat{L}\rho _0\widehat{L}^{}\widehat{L}^{}\widehat{L}\rho _0\rho _0\widehat{L}^{}\widehat{L}\right)dt`$ (49)
$`{\displaystyle \frac{4i\kappa E^2}{\gamma _c^2}}[\widehat{Z},\rho _0]dt+{\displaystyle \frac{2\kappa E}{\gamma _c}}(\widehat{Z}\rho _1+\rho _1^{}\widehat{Z})dt+\gamma _c\rho _2dt`$
$`+\sqrt{\gamma _ce_d}\left(\rho _1+\rho _1^{}\rho _0\mathrm{Tr}\{\rho _1+\rho _1^{}\}\right)dW_t,`$
$`d\rho _1`$ $`=`$ $`\left({\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho _1]+2\widehat{L}\rho _1\widehat{L}^{}\widehat{L}^{}\widehat{L}\rho _1\rho _1\widehat{L}^{}\widehat{L}\right)dt`$ (50)
$`i\kappa \widehat{Z}\rho _1dt{\displaystyle \frac{4i\kappa E^2}{\gamma _c^2}}[\widehat{Z},\rho _1]dt{\displaystyle \frac{2\kappa E}{\gamma _c}}(\widehat{Z}\rho _0\rho _2\widehat{Z})dt(\gamma _c/2)\rho _1dt`$
$`+\sqrt{\gamma _ce_d}\left(\rho _2\rho _1\mathrm{Tr}\{\rho _1+\rho _1^{}\}\right)dW_t,`$
$`d\rho _2`$ $`=`$ $`\left({\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho _2]+2\widehat{L}\rho _2\widehat{L}^{}\widehat{L}^{}\widehat{L}\rho _2\rho _2\widehat{L}^{}\widehat{L}\right)dt`$ (51)
$`\left(i\kappa +{\displaystyle \frac{4i\kappa E^2}{\gamma _c^2}}\right)[\widehat{Z},\rho _2]dt{\displaystyle \frac{2\kappa E}{\gamma _c}}(\widehat{Z}\rho _1^{}+\rho _1\widehat{Z})dt\gamma _c\rho _2dt`$
$`\sqrt{\gamma _ce_d}\rho _2\mathrm{Tr}\{\rho _1+\rho _1^{}\}dW_t.`$
Both $`\rho _1`$ and $`\rho _2`$ contain damping terms, which imply that they will remain small at all times provided $`\kappa \widehat{Z}`$ is sufficiently small compared to $`\gamma _c`$. (This also implies that our ansatz is reasonable for sufficiently small $`\kappa `$.)

By making use of the above equations, we can find the evolution equation for the reduced density matrix $`\rho _{SZ}`$:

$`d\rho _{SZ}(t)`$ $`=`$ $`d\rho _0(t)+d\rho _2(t)`$ (52)
$`=`$ $`\left({\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho _{SZ}]+2\widehat{L}\rho _{SZ}\widehat{L}^{}\widehat{L}^{}\widehat{L}\rho _{SZ}\rho _{SZ}\widehat{L}^{}\widehat{L}\right)dt`$
$`{\displaystyle \frac{4i\kappa E^2}{\gamma _c^2}}[\widehat{Z},\rho _{SZ}]dt+{\displaystyle \frac{2\kappa E}{\gamma _c}}[\widehat{Z},\rho _1\rho _1^{}]dti\kappa [\widehat{Z},\rho _2]dt`$
$`+\sqrt{\gamma _ce_d}\left(\rho _1+\rho _1^{}\rho _{SZ}\mathrm{Tr}\{\rho _1+\rho _1^{}\}\right)dW_t.`$
If we keep only terms to second order in $`\kappa \widehat{Z}`$ we can neglect the $`\rho _2`$ term. This leaves only the terms proportional to $`\rho _1\pm \rho _1^{}`$, which we need only know to leading order in $`\kappa \widehat{Z}`$. Provided (as we have already assumed) that the cantilever moves slowly compared to the timescale set by $`\gamma _c`$ and that $`\kappa \widehat{Z}`$ can be treated as small, then to leading order $`d\rho _1`$ vanishes; $`\rho _1`$ remains in an approximate equilibrium state. If we make use of this assumption we can (again to leading order) solve for $`\rho _1\pm \rho _1^{}`$:

$`\rho _1+\rho _1^{}`$ $``$ $`{\displaystyle \frac{4\kappa E}{\gamma _c^2}}\{\widehat{Z},\rho _{SZ}\},`$
$`\rho _1\rho _1^{}`$ $``$ $`{\displaystyle \frac{4\kappa E}{\gamma _c^2}}[\widehat{Z},\rho _{SZ}],`$ (53)
which when inserted into (52) gives us a closed evolution equation for $`\rho _{SZ}`$:

$`d\rho _{SZ}(t)`$ $`=`$ $`\left({\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{SZ}^{}(t),\rho _{SZ}]+2\widehat{L}\rho _{SZ}\widehat{L}^{}\widehat{L}^{}\widehat{L}\rho _{SZ}\rho _{SZ}\widehat{L}^{}\widehat{L}\right)dt`$ (54)
$`{\displaystyle \frac{4i\kappa E^2}{\gamma _c^2}}[\widehat{Z},\rho _{SZ}]dt{\displaystyle \frac{8\kappa ^2E^2}{\gamma _c^3}}[\widehat{Z},[\widehat{Z},\rho _{SZ}]]dt`$
$`+\sqrt{\gamma _ce_d}{\displaystyle \frac{4\kappa E}{\gamma _c^2}}\left(\widehat{Z}\rho _{SZ}+\rho _{SZ}\widehat{Z}2\rho _{SZ}\mathrm{Tr}\{\rho _{SZ}\}\right)dW_t.`$
(Note that we have absorbed a factor of $`1`$ into $`dW_t`$.)

Examining the terms in (54), we see that by eliminating the cavity mode we get another effective term in the Hamiltonian, and another Lindblad operator. We can therefore write this stochastic master equation in the form

$`d\rho _{SZ}(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{eff}}(t),\rho _{SZ}]dt+{\displaystyle \underset{j=1}{\overset{2}{}}}\left(2\widehat{L}_j\rho _{SZ}\widehat{L}_j^{}\widehat{L}_j^{}\widehat{L}_j\rho _{SZ}\rho _{SZ}\widehat{L}_j^{}\widehat{L}_j\right)dt`$ (55)
$`+\sqrt{e_d/2}\left((\widehat{L}_2\widehat{L}_2)\rho _{SZ}+\rho _{SZ}(\widehat{L}_2\widehat{L}_2)\right)dW_t,`$
where we now make the definitions

$`\widehat{L}_1`$ $`=`$ $`\sqrt{\gamma _m/2}\left((1/\mathrm{})\widehat{Z}+i(\mathrm{}/\mathrm{})\widehat{p}\right)`$
$`\widehat{L}_2`$ $`=`$ $`\sqrt{8\kappa ^2E^2/\gamma _c^3}\widehat{Z},`$
$`\widehat{H}_{\mathrm{eff}}(t)`$ $`=`$ $`{\displaystyle \frac{1}{2m}}\widehat{p}^2+{\displaystyle \frac{m\omega _m^2}{2}}\widehat{Z}^22\eta (f(t)/\lambda (t))\widehat{Z}\widehat{S}_z^{}`$ (56)
$`+{\displaystyle \frac{4\kappa E^2}{\gamma _c^2}}\widehat{Z}+(\gamma _m/2)(\widehat{Z}\widehat{p}+\widehat{p}\widehat{Z}).`$
Note that the term $`4\kappa E^2\widehat{Z}/\gamma _c^2`$ is a constant force, which just displaces the equilibrium position of the cantilever. It can be eliminated simply by changing the origin of $`\widehat{Z}`$, and is in any case small for reasonable values of the parameters. The output from the homodyne measurement now corresponds to a measurement of the cantilever position $`\widehat{Z}`$:

$$I_c(t)=\beta \left(\frac{8e_d\kappa E}{\gamma _c}\widehat{Z}+\sqrt{\gamma _ce_d}\frac{dW_t}{dt}\right).$$
(57)
As we shall see in the next section, we can further unravel this stochastic master equation (55) into a stochastic Schrödinger equation for pure states. This further unraveling provides considerable improvement in numerical efficiency, though it does not represent an actual measurement process.

## 6 Pure state unraveling

The stochastic master equation (55) represents the evolution of the cantilever-spin system, conditioned on the photocurrent measurement record $`I_c(t)`$. If we averaged over all possible measurement records, the $`dW_t`$ terms would average to zero, and we would be left with an ordinary deterministic master equation for the cantilever and spin. It is for this reason that the stochastic master equation is therefore often referred to as an unraveling of the average master equation.

For numerical purposes, it is often much easier to solve an equation for a pure state vector rather than a density matrix . It is therefore useful to unravel equation (55) still further to an equation which preserves pure states. We do this by introducing a second stochastic process.

First, let us idealize to perfect detector efficiency $`e_d=1`$. We then introduce the new master equation

$`d\rho _{SZ}(t)`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}[\widehat{H}_{\mathrm{eff}}(t),\rho _{SZ}]dt+{\displaystyle \underset{j=1}{\overset{2}{}}}\left(2\widehat{L}_j\rho _{SZ}\widehat{L}_j^{}\widehat{L}_j^{}\widehat{L}_j\rho _{SZ}\rho _{SZ}\widehat{L}_j^{}\widehat{L}_j\right)dt`$ (58)
$`+{\displaystyle \underset{j=1}{\overset{2}{}}}\left((\widehat{L}_j\widehat{L}_j)\rho _{SZ}+\rho _{SZ}(\widehat{L}_j\widehat{L}_j)\right)dW_{jt},`$
where the Hamiltonian and Lindblad operators are the same as in (56) and we now have two independent noise processes represented by stochastic differential variables $`dW_{1t}`$ and $`dW_{2t}`$ which satisfy

$$M[dW_{jt}]=0,dW_{it}dW_{js}=\delta (ts)\delta _{ij}dt.$$
(59)
If we take the mean of (58) over $`dW_{1t}`$ we recover equation (55). We can think of the additional stochastic process as representing a fictitious additional measurement, whose outcome we average over to recover the state which is conditioned on the actual measurement.

However, equation (58) has a great advantage over (55). If $`\rho _{SZ}`$ is initially a pure state $`\rho _{SZ}=|\psi _{SZ}\psi _{SZ}|`$, it will remain a pure state at all times, the state of course depending on the stochastic processes $`W_1`$ and $`W_2`$. We can recover the solution of (55) by averaging

$$\rho _{SZ}(t)=M_{W_1}[|\psi _{SZ}(t)\psi _{SZ}(t)|].$$
(60)
It would be useful to replace equation (58) with an explicit evolution equation for $`|\psi _{SZ}`$ instead of $`\rho _{SZ}`$. This equation is the quantum state diffusion equation with real noise :

$`d|\psi _{SZ}`$ $`=`$ $`{\displaystyle \frac{i}{\mathrm{}}}\widehat{H}_{\mathrm{eff}}(t)|\psi _{SZ}dt+{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \underset{j}{}}\left(2\widehat{L}_j^{}\widehat{L}_j\widehat{L}_j^{}\widehat{L}_j|\widehat{L}_j|^2\right)|\psi _{SZ}dt`$ (61)
$`+{\displaystyle \frac{1}{\sqrt{2}}}{\displaystyle \underset{j}{}}\left(\widehat{L}_j\widehat{L}_j\right)|\psi _{SZ}dW_{jt}.`$
The nonlinearity of this equation arises to preserve the norm.

## 7 Numerical simulation

We have simulated this system using the C++ quantum state diffusion library  to numerically solve both the unitary evolution with Hamiltonian (30) and the stochastic equation (61). All of the figures in this paper were generated using this software.

We chose our parameters based on those used by Berman et al., . These values are (in arbitrary units):

$`\mathrm{}`$ $`=`$ $`\omega _m=m=1,`$
$`\eta `$ $`=`$ $`0.3,`$
$`\epsilon `$ $`=`$ $`400.0,`$
$`\gamma _m`$ $`=`$ $`\omega _m/Q=0.00001,`$
$`k_BT`$ $`=`$ $`10000.0,`$ (62)
where $`Q`$ is the quality factor of the cantilever. The driving force $`f(t)`$ takes the form

$$f(t)=\{\begin{array}{cc}6000+300t\hfill & \text{if }0t20\text{;}\hfill \\ 1000\mathrm{sin}(t20)\hfill & \text{if }t>20\text{.}\hfill \end{array}$$
(63)
If we make contact with physical values for actual cantilevers used in experiments, we have $`\omega _m10^5\mathrm{s}^1`$ and $`m10^{12}\mathrm{kg}`$. The value of $`k_BT`$ above then corresponds to a temperature of around $`0.1\mathrm{K}`$, which is within the bounds of experimental feasibility, though rather lower than the temperatures used in current experiments (around 3K) . Since $`\eta =(g\mu /2)(B_z/Z)_0`$, the value of $`\eta `$ corresponds to a field gradient of about $`1.5\times 10^7\mathrm{T}/\mathrm{m}`$, which is higher than current experiments by roughly two orders of magnitude , but hopefully this too will improve with time. The cantilever would undergo displacements of about a nanometer.

Alternatively, rather than increasing the field gradient we could achieve similar numbers by lowering the spring constant of the cantilever, for instance by shrinking the mass of the cantilever. Lowering the mass by a factor of 100 has the same relative effect on $`\eta `$ as increasing the field gradient by a factor of ten.

We then might ask about realistic parameters for the monitoring. A typical cavity size $`L`$ is about a micrometer, with a laser frequency of $`\omega _c1.4\times 10^{15}\mathrm{s}^1`$. This cavity is generally quite lossy; reasonable quality factors might be in the range $`Q_c`$10–100. The parameter $`E`$ is a function of the laser power, $`E=\sqrt{P\gamma _c/\mathrm{}\omega _c}=\sqrt{P/\mathrm{}Q_c}`$. For $`P1\mu \mathrm{W}`$ and $`Q_c100`$ we have $`E10^{13}\mathrm{s}^1`$. The coupling between the cantilever and the cavity is given by a geometric factor $`\kappa =\omega _c/L1.4\times 10^{21}(\mathrm{m}\mathrm{s})^1`$. In arbitrary units, this gives coefficients

$`{\displaystyle \frac{8\kappa E}{\gamma _c}}`$ $`=`$ $`1.9\times 10^3,`$
$`{\displaystyle \frac{4\kappa E^2}{\gamma _c^2}}`$ $`=`$ $`7\times 10^2,`$
$`\sqrt{{\displaystyle \frac{8\kappa ^2E^2}{\gamma _c^3}}}`$ $`=`$ $`0.07.`$ (64)
The first value is the multiplier in (57); the second gives the equilibrium displacement of the cantilever; the third is the coefficient of the Lindblad operator $`\widehat{L}_2`$.

One question we can now easily address is how quickly the state of the spin collapses onto eigenstates of $`\widehat{S}_z^{}`$. In figure 5 we plot $`\widehat{S}_z^{}`$ for ten different trajectories. We see that in all ten cases the spin converged to $`\pm 1/2`$ quite quickly, before $`t=80`$.

If we compare this to the results of figure 4, we see that the spin state collapses rather more quickly than the cantilever oscillations can respond. We only get a clear output signal when the two phases are well separated, which does not occur until nearly $`t=150`$. Generically, the difficulty of collapsing the spin state is much less than the difficulty of obtaining an unequivocal readout.

The curves depicted in figure 4 are idealized, without the measurement noise which will always be present in the output current (42) or (57). In figure 6 we show what actual output would look like for the set of parameters we are discussing. Note that even with the noise, the two phases (representing spin up and spin down) are clearly distinguishable. In the next section, we derive an expression for the signal-to-noise ratio in more general situations.

## 8 Signal-to-noise ratio

Since we have to detect the effect of very weak force on the cantilever by the single spin, we need very high resolution for the cantilever position measurements and a good control of the various noise sources in the MRFM device. As described in section 2, the tiny displacement of the cantilever is measured by a fiber-optic interferometer as a phase shift of the interference fringes. We shall analyze the quantum and thermal noise in this homodyne measurement scheme.

The Hamiltonian for the combined system of the spin, cantilever and cavity mode, excluding coupling to the environments, in the spin rotating frame is

$`\widehat{H}`$ $`=`$ $`\widehat{H}_Z2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}\widehat{Z}\widehat{S}_z^{}+\mathrm{}\omega _c\widehat{a}^{}\widehat{a}`$ (65)
$`+\mathrm{}E(\widehat{a}^{}e^{i\omega _0t}+\widehat{a}e^{i\omega _0t})+\mathrm{}\kappa \widehat{a}^{}\widehat{a}\widehat{Z}.`$
Here $`\omega _c`$ is the optical frequency of the cavity mode, $`\omega _0\omega _c`$ is the driving frequency of the external laser and other terms and parameters have been described in section 5. The master equation approach in section 4 is valid in high or medium temperature case. Here we analyze the noise in the Heisenberg picture, using the quantum Langevin equation approach that is valid at any temperature.

Using standard techniques,  the reservoir (environmental) variables may be eliminated, in the interaction picture with respect to $`\mathrm{}\omega _0\widehat{a}^{}\widehat{a}`$, to give the quantum Langevin equations describing the dynamics of the whole system:

$`{\displaystyle \frac{d\widehat{Z}(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{m}}\widehat{p}(t),`$ (66)
$`{\displaystyle \frac{d\widehat{p}(t)}{dt}}`$ $`=`$ $`m\omega _m^2\widehat{Z}(t){\displaystyle \frac{\mathrm{\Gamma }}{m}}\widehat{p}(t)\mathrm{}\kappa \widehat{a}^{}(t)\widehat{a}(t)+\widehat{𝒲}(t)+2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}\widehat{S}_z^{}(t),`$ (67)
$`{\displaystyle \frac{d\widehat{a}(t)}{dt}}`$ $`=`$ $`(i\omega _ci\omega _0+{\displaystyle \frac{\gamma _c}{2}})\widehat{a}(t)i\kappa \widehat{Z}(t)\widehat{a}(t)iE+\sqrt{\gamma _c}\widehat{a}_{\mathrm{in}}(t),`$ (68)
$`{\displaystyle \frac{d\widehat{S}_z^{}(t)}{dt}}`$ $`=`$ $`0,`$ (69)
$`{\displaystyle \frac{d\widehat{S}_x^{}(t)}{dt}}`$ $`=`$ $`2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}\widehat{Z}(t)\widehat{S}_y^{}(t),`$ (70)
$`{\displaystyle \frac{d\widehat{S}_y^{}(t)}{dt}}`$ $`=`$ $`2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}\widehat{Z}(t)\widehat{S}_x^{}(t).`$ (71)
In the equations, the usual optical input noise operator $`\widehat{a}_{\mathrm{in}}(t)`$ is associated with the vacuum fluctuations of the continuum of electromagnetic modes outside the cavity and its correlation function is given by

$$\widehat{a}_{\mathrm{in}}(t)\widehat{a}_{\mathrm{in}}^{}(t^{})=\delta (tt^{}).$$
(72)
The random force $`\widehat{𝒲}(t)`$ describes the thermal noise motion (quantum Brownian motion) of the cantilever at temperature $`T`$. For the case of an Ohmic environment, the thermal random force correlation is given by

$$\widehat{𝒲}(t)\widehat{𝒲}(t^{})=\frac{\mathrm{}\mathrm{\Gamma }}{\pi }[_r(tt^{})+i_i(tt^{})],$$
(73)
where

$`_r(t)`$ $`=`$ $`{\displaystyle _0^\mathrm{\Omega }}𝑑\omega \omega \mathrm{cos}(\omega t)\mathrm{coth}\left({\displaystyle \frac{\mathrm{}\omega }{2k_BT}}\right),`$ (74)
$`_i(t)`$ $`=`$ $`{\displaystyle _0^\mathrm{\Omega }}𝑑\omega \omega \mathrm{sin}(\omega t),`$ (75)
with $`\mathrm{\Omega }`$ the frequency cutoff of the reservoir spectrum. Without the presence of the external driving force from the spin, the cantilever-cavity system can be characterized by a semi-classical steady state with a new equilibrium position for the cantilever, displaced by $`Z_{st}=\kappa |\alpha _{st}|^2/(m\omega _m^2)`$ with respect to that with no external driving laser field, and the cavity mode in a coherent state $`|\alpha _{st}`$ with the amplitude given by

$$\alpha _{st}=\frac{iE}{\gamma _c/2+i\mathrm{\Delta }},$$
(76)
where $`\mathrm{\Delta }=\omega _c\omega _0\kappa ^2|\alpha _{st}|^2/(m\omega _m^2)`$ is the cavity mode detuning. By adjusting either $`\omega _0`$ or $`\omega _c`$, the detuning can be set to zero $`\mathrm{\Delta }=0`$. As a result, $`\alpha _{st}=\alpha _0=2iE/\gamma _c`$. Linearizing the quantum Langevin equations about the steady-state values and renaming with $`\widehat{Z}(t)`$, $`\widehat{a}(t)`$ the operators describing the quantum fluctuations around the classical steady state, we obtain

$`{\displaystyle \frac{d\widehat{Z}(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{1}{m}}\widehat{p}(t),`$ (77)
$`{\displaystyle \frac{d\widehat{p}(t)}{dt}}`$ $`=`$ $`m\omega _m^2\widehat{Z}(t){\displaystyle \frac{\mathrm{\Gamma }}{m}}\widehat{p}(t)\mathrm{}\kappa [\alpha _0\widehat{a}^{}(t)+\alpha _0^{}\widehat{a}(t)]`$ (78)
$`+\widehat{𝒲}(t)+2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}\widehat{S}_z^{}(t),`$
$`{\displaystyle \frac{d\widehat{a}(t)}{dt}}`$ $`=`$ $`{\displaystyle \frac{\gamma _c}{2}}\widehat{a}(t)i\kappa \alpha _0\widehat{Z}(t)+\sqrt{\gamma _c}\widehat{a}_{\mathrm{in}}(t),`$ (79)
$`{\displaystyle \frac{d\widehat{S}_z^{}(t)}{dt}}`$ $`=`$ $`0,`$ (80)
$`{\displaystyle \frac{d\widehat{S}_x^{}(t)}{dt}}`$ $`=`$ $`2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}[Z_{st}+\widehat{Z}(t)]\widehat{S}_y^{}(t),`$ (81)
$`{\displaystyle \frac{d\widehat{S}_y^{}(t)}{dt}}`$ $`=`$ $`2\eta {\displaystyle \frac{f(t)}{\lambda (t)}}[Z_{st}+\widehat{Z}(t)]\widehat{S}_x^{}(t).`$ (82)
In the bad cavity limit where $`\gamma _c\omega _m,(\mathrm{\Gamma }/m),\kappa \widehat{Z}`$ \[i.e., set $`(d\widehat{a}(t)/dt)=0`$ in (79)\], the dynamics of the field quadrature, $`\widehat{a}^{}(t)+\widehat{a}(t)`$, adiabatically follows that of the cantilever position:

$$\widehat{a}^{}(t)+\widehat{a}(t)=i\frac{4\kappa \alpha _0}{\gamma _c}\widehat{Z}(t)+\frac{2}{\sqrt{\gamma _c}}[\widehat{a}_{\mathrm{in}}(t)+\widehat{a}_{\mathrm{in}}^{}(t)].$$
(83)
Thus monitoring this field quadrature of the cavity mode via homodyne measurement corresponds to a measurement of the cantilever position and hence the state of the spin. Using the usual input-output relation,

$$\widehat{a}_{\mathrm{out}}(t)=\sqrt{\gamma _c}\widehat{a}(t)\widehat{a}_{\mathrm{in}},$$
(84)
we may define an operator corresponding to the output current

$`\widehat{I}_{\mathrm{out}}(t)`$ $`=`$ $`\sqrt{\gamma _c}\beta [\widehat{a}_{\mathrm{out}}(t)+\widehat{a}_{\mathrm{out}}^{}(t)]`$ (85)
$`=`$ $`\beta \{\gamma _c[\widehat{a}(t)+\widehat{a}^{}(t)]\sqrt{\gamma _c}[\widehat{a}_{\mathrm{in}}(t)+\widehat{a}_{\mathrm{in}}^{}(t)]\}.`$
Equation (85) is similar to (42) with $`e_d=1`$. By substituting (83) into (85), the resultant output current in the bad cavity limit is given by

$$\widehat{I}_{\mathrm{out}}(t)=\beta \left(\frac{8\kappa E}{\gamma _c}\widehat{Z}(t)+\sqrt{\gamma _c}[\widehat{a}_{\mathrm{in}}(t)+\widehat{a}_{\mathrm{in}}^{}(t)]\right).$$
(86)
This equation is also similar to (57) with $`e_d=1`$, obtained from master equation approach.

The Langevin equations for $`\widehat{S}_x^{}`$ and $`\widehat{S}_y^{}`$ effectively decouple from the other equations, since they do not appear on the right-hand-side of the equations for the other variables. Because of this, they have no effect in our estimate of the signal-to-noise ratio, and we shall drop them henceforth. Taking a Fourier transform of the linearized Langevin equations, we find, from (85), the Fourier component of the output current as

$`\widehat{I}_{\mathrm{out}}(\omega )`$ $`=`$ $`{\displaystyle \frac{\beta \gamma _c}{(i\omega \gamma _c/2)}}\{(i\omega +{\displaystyle \frac{\gamma _c}{2}})[\widehat{a}_{\mathrm{in}}(\omega )+\widehat{a}_{\mathrm{in}}^{}(\omega )]`$ (87)
$`{\displaystyle \frac{2i\kappa \alpha _0\sqrt{\gamma _c}}{m(\omega _m^2\omega ^2i\mathrm{\Gamma }\omega /m)}}[{\displaystyle \frac{\mathrm{}\kappa \alpha _0\sqrt{\gamma _c}}{(i\omega \gamma _c/2)}}[\widehat{a}_{\mathrm{in}}^{}(\omega )\widehat{a}_{\mathrm{in}}(\omega )]`$
$`+\widehat{𝒲}(\omega )+G(\omega )\widehat{S}_z^{}]]\},`$
where $`G(\omega )`$ is the Fourier transform of $`G(t)=2\eta f(t)/\lambda (t)`$. The Fourier component of the mean output current signal is then given by

$$|\widehat{I}_{\mathrm{out}}(\omega )|=\beta \gamma _c\left(\frac{2\kappa \sqrt{\gamma _c}|\alpha _0|}{m}\right)\frac{|G(\omega )|}{|D(\omega )|}\widehat{S}_z^{},$$
(88)
where

$$D(\omega )=\left(i\omega \frac{\gamma _c}{2}\right)\left(\omega _m^2\omega ^2i\omega \frac{\mathrm{\Gamma }}{m}\right).$$
(89)
The output current noise power density spectrum is defined as

$`S_{\mathrm{out}}(\omega )`$ $`=`$ $`\left\{{\displaystyle 𝑑\tau e^{i\omega \tau }\widehat{I}_{\mathrm{out}}(t)\widehat{I}_{\mathrm{out}}(t+\tau )_{G(t)=0}}\right\}_t`$ (90)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\left\{{\displaystyle 𝑑\omega ^{}e^{i(\omega +\omega ^{})t}\widehat{I}_{\mathrm{out}}(\omega ^{})I_{\mathrm{out}}(\omega )_{G(\omega )=0}}\right\}_t,`$
where the subscript $`G(t)=0`$ means evaluation in the absence of the external driving force from the spin and $`\{\mathrm{}\}_t`$ denotes the time average over $`t`$. To calculate this noise spectrum, the Fourier transform of the noise correlation functions (72)–(75) is needed and given by

$`\widehat{a}_{\mathrm{in}}(\omega )\widehat{a}_{\mathrm{in}}^{}(\omega ^{})`$ $`=`$ $`2\pi \delta (\omega +\omega ^{}),`$ (91)
$`\widehat{𝒲}(\omega )\widehat{𝒲}(\omega ^{})`$ $`=`$ $`2\pi \mathrm{}\mathrm{\Gamma }\omega \left[1+\mathrm{coth}\left({\displaystyle \frac{\mathrm{}\omega }{2k_BT}}\right)\right]\delta (\omega +\omega ^{}),`$ (92)
where in obtaining (92) the infinite frequency cutoff limit of the Ohmic thermal reservoir spectrum, $`\mathrm{\Omega }\mathrm{}`$, has been assumed. After some calculations, one can then obtain the output noise spectrum as

$`S_{\mathrm{out}}(\omega )`$ $`=`$ $`\beta ^2\gamma _c^2\{1+4\left({\displaystyle \frac{\mathrm{}\kappa ^2\gamma _c|\alpha _0|^2}{m}}\right)^2{\displaystyle \frac{1}{[(\gamma _c/2)^2+\omega ^2]|D(\omega )|^2}}`$ (93)
$`+4\left({\displaystyle \frac{\kappa ^2\gamma _c|\alpha _0|^2\mathrm{\Gamma }}{m^2}}\right){\displaystyle \frac{\mathrm{}\omega }{|D(\omega )|^2}}\mathrm{coth}\left({\displaystyle \frac{\mathrm{}\omega }{2k_BT}}\right)\}.`$
The first term in (93), independent of frequency, is the contribution from the shot noise of the photons. The next term is the “back-action” noise on the position of the cantilever by the radiation (photons). This back action is due to the random way in which photons bounce off the cantilever. The final term is the thermal noise, due to the thermal Brownian-motion fluctuation of the cantilever. Equation (93) is valid at all temperatures. The assumptions made in its derivation are the linearization around the semi-classical steady state and the infinite frequency cutoff $`\mathrm{\Omega }\mathrm{}`$. The high (or medium) temperature limit $`\mathrm{}\omega _mk_BT`$ can be obtained by approximating

$$\mathrm{}\omega \mathrm{coth}\left(\frac{\mathrm{}\omega }{2k_BT}\right)2k_BT+\frac{\mathrm{}^2\omega ^2}{6k_BT}.$$
(94)
We plot these three contributions to the noise for the simulation parameters given in section 7. We see that at the oscillator resonance $`\omega _m=1`$, thermal noise dominates.

Let us define the signal-to-noise ratio per root Hertz as

$$\mathrm{SNR}(\omega )=\frac{|\widehat{I}_{\mathrm{out}}(\omega )|}{\sqrt{S_{\mathrm{out}}(\omega )}}.$$
(95)
We are interested in evaluating SNR$`(\omega )`$ at the frequency equal to the cantilever vibration frequency, $`\omega =\omega _m`$. Note that

$$\frac{1}{|D(\omega _m)|}=\frac{1}{[(\gamma _c/2)^2+\omega _m^2]^{1/2}}\left(\frac{Q}{\omega _m^2}\right),$$
(96)
where the quality factor $`Q=m\omega _m/\mathrm{\Gamma }`$. As a result, the mean output current signal (88) at $`\omega =\omega _m`$ is enhanced by a factor of $`Q\gamma _c/[(\gamma _c/2)^2+\omega _m^2]^{1/2}`$ as compared with the $`\omega =0`$ case. However, a similar enhancement occurs in the back-action noise and the thermal noise terms. In other words, driving the cantilever at $`\omega =\omega _m`$ amplifies not only its vibration amplitudes due to the the driving force, but also the noise amplitude due to the back-action radiation pressure and thermal Brownian motion (see Fig. 7). We find $`\mathrm{SNR}(\omega =\omega _m)`$ can be written as

$$\mathrm{SNR}(\omega _m)=\frac{|G(\omega _m)|\widehat{S}_z^{}}{\sqrt{N(\omega _m)}}$$
(97)
where

$$N(\omega _m)=\frac{[(\gamma _c/2)^2+\omega _m^2]}{4\kappa ^2\gamma _c|\alpha _0|^2}\left(\frac{m\omega _m^2}{Q}\right)^2+\frac{\mathrm{}^2\kappa ^2\gamma _c|\alpha _0|^2}{[(\gamma _c/2)^2+\omega _m^2]}+\mathrm{\Gamma }\mathrm{}\omega _m\mathrm{coth}\left(\frac{\mathrm{}\omega _m}{2k_BT}\right).$$
(98)
We may set $`\widehat{S}_z^{}=\pm (1/2)`$ to estimate the signal-to-noise ratio per root Hz, corresponding respectively to the spin in the two different states in the rotating frame.

Because the driving force $`f(t)`$ is periodic, $`G(\omega )`$ is equal to a sum of delta functions at $`\omega =\omega _m,3\omega _m,5\omega _m,\mathrm{}`$. Averaging over a small interval about $`\omega _m`$, we can integrate over the delta function to get a value (for our simulation parameters) of $`\mathrm{SNR}(\omega _m)220\mathrm{s}^{1/2}`$. Thus, given a bandwidth of about 1Hz, this should be easily detectable by our measurement scheme. As mentioned in section 7, we have assumed a magnetic field gradient roughly two orders of magnitude greater than current experiments, and a much lower temperature. A single spin, therefore, would be below the edge of detectability by current experimental techniques. Steady improvement in the field strength, temperature and spring constant of these experiments, however, should soon make single-spin measurement possible.

If the dominant noise source in MRFM comes from the thermal Brownian motion of the cantilever, we can estimate the minimum detectable force (when the signal-to-noise ratio is one) by keeping only the last term of (98). In this case, with a measurement bandwith $`\mathrm{\Delta }\nu `$, we obtain from (97), (98) and (94) the usual expression of the minimum detectable force at the high-temperature limit ($`\mathrm{}\omega _mk_BT`$) as

$$F_{\mathrm{min}}=\sqrt{N(\omega _m)\mathrm{\Delta }\nu }=\sqrt{\frac{2kk_BT\mathrm{\Delta }\nu }{Q\omega _m}},$$
(99)
where $`k=m\omega _m^2`$ is the spring constant of the cantilever. We see, then, that improvement can come either from raising the force (by increasing the field gradient), lowering the temperature, or lowering the spring constant.

## 9 Conclusion

We have derived an approximate description of single-spin measurement by magnetic resonance force microscopy, including both thermal noise and measurement back-action, and used it to produce numerical simulations of a single-spin measurement. These simulations use the quantum trajectory method for open quantum systems. The parameters we assumed for this simulation were somewhat optimistic; but given the steady improvement in experimental technique, we believe that measurements of this type will be possible in the near future.

Single-spin measurements would be very useful in the construction of solid-state quantum computers, in which the spin of an electron represents a single qubit of information. Given the great interest in solid-state implementations as a possibly scalable realization of quantum computers, finding practical ways to measure single spins would be very useful. The results of our simulations suggest that magnetic resonance force microscopy is a very promising approach to this difficult problem.

## Acknowledgment

HSG would like to thank G. P. Berman, G. J. Milburn, D. V. Pelekhov and P. C. Hammel for useful discussions. HSG would also like to acknowledge support from the Hewlett-Packard Fellowship. TAB was supported by the Martin A. and Helen Chooljian Membership in Natural Sciences, and DOE Grant No. DE-FG02-90ER40542.