PHYSICAL REVIEW A 84, 043824 (2011)

Full counting statistics of energy fluctuations in a driven quantum resonator
A. A. Clerk
Department of Physics, McGill University, Montreal, Quebec, Canada H3A 2T8(Received 3 June 2011; published 14 October 2011)
We consider the statistics of time-integrated energy fluctuations of a driven bosonic single-mode resonator,
as measured by a quantum nondemolition (QND) detector, using the standard Keldysh prescription to define
higher moments. We find that, due to an effective cascading of fluctuations, these statistics are surprisingly
nonclassical: the low-temperature, quantum probability distribution is not equivalent to the high-temperatureclassical distribution evaluated at some effective temperature. Moreover, for a sufficiently large drive detuning
and low temperatures, the Keldysh-ordered quasiprobability distribution characterizing these fluctuations fails to
be positive-definite; this is similar to the full counting statistics of charge in superconducting systems. We argue
that this indicates a kind of nonclassical behavior akin to that tested by Leggett-Garg inequalities.
DOI: 10.1103/PhysRevA.84.043824

PACS number(s): 42.50.Lc, 72.70.+m

I. INTRODUCTION

Our study also addressed the non-Gaussian nature of the driven
energy fluctuations by calculating the third central moment
(δ m̂)3 ; surprisingly, we found that, while this quantity isalways positive classically, it could become negative in the
low-temperature quantum limit. As such, the third moment is
far more sensitive to classical-quantum differences than thesecond moment.
To fully understand the significance of this result, one needs
to consider the full probability distribution characterizing the
low-frequency fluctuations of n̂ and compare its form in theclassical and quantum limits. This is the objective of thispaper: we calculate the distribution P (m) in the long-time
limit using the standard Keldysh operator ordering [10–12].
We find that the anomalous negative value of (δ m̂)3  resultsfrom a kind of cascaded fluctuation effect [13,14], which
can be heuristically attributed to a correlated fluctuation in
the resonator temperature. We also find that this negativeskewness is a precursor of something rather dramatic: inthe quantum limit, the fluctuations of m are most naturally
described by a quasiprobability distribution P (m) which is
not positive definite. Such negative counting statistics havebeen encountered before in the study of charge transfer insuperconducting systems [12,15]; their interpretation requires
some care. As we discuss in some detail, they are indicativeof nonclassical temporal correlations and are thus somewhat
similar to having violated a Leggett-Garg inequality [16].
Detecting these effects thus represents a new way of detecting
nonclassical behavior in a driven quantum resonator.
The remainder of this paper is structured as follows: In
Sec. II, we introduce our basic model, and present our main
results for the generating function of P (m); we also give acompact review of the Keldysh ordering of higher moments
for those not familiar with this topic. In Sec. III, we discuss
the form of the distribution in the classical limit. Section IV
is dedicated to the distribution in the quantum limit, while
Sec. V is devoted to interpreting the negative quasiprobabilities
which emerge. Finally, in Sec. VI we discuss issues related tothe measurement of these effects. An appendix is included
which shows how the Keldysh operator ordering emerges
naturally in the proposed experimental realization of Ref. [9],where m is measured by using homodyne interferometry todetect the frequency shift of an auxiliary cavity. Finally, we
mention work on related problems. A fermionic analog of

The statistics of photon fluctuations in various settings isby now an almost textbook phenomena. Most familiar are
the statistics that would be measured by a photodetector.
Relatively less attention has been paid to photon fluctuations
for the case where the detection is done in a nondemolition
manner, meaning that energy quanta are measured without
destroying them. The study of such fluctuations is not just atheoretical curiosity, as quantum nondemolition (QND) photon
detection is experimentally feasible both optically [1] as wellas in both cavity QED systems [2–4] and superconducting
circuit QED systems [5,6]. In the latter systems, one uses
dispersive interactions to detect the photon number inside acavity. QND detection of phonon number in a mechanicalresonator may also soon be possible in optomechanicalsystems [7,8], where the energy of a single vibrational modein a mechanical resonator is directly coupled to the frequencyof an optical cavity.
Motivated by developments in optomechanics, we recently
investigated the low-frequency energy fluctuations of a driven,
damped, single-mode harmonic resonator, focusing on thepossibility of measuring these fluctuations nondestructively
using an optomechanical cavity [9]. In the zero-temperature
quantum limit, the instantaneous state of such a resonator issimply a coherent state, yielding a Poissonian distribution ofphonon number. Our focus was instead on understanding how
the mechanical phonon number n̂ fluctuated in time. Thesefluctuations are characterized by a power spectral density
Snn [ω] or, equivalently, by the second central moment of thetime-integrated phonon number m̂: tm̂ ≡dt  n̂(t  ).
(1)
0

As could be anticipated, both Snn [ω] and (δ m̂)2  = m̂2  −m̂2 have a low-frequency “shot-noise” term proportional tothe average number of phonons induced by the drive, n̄dr ;detecting this shot-noise contribution would be direct evidence
for the quantization of the mechanical resonator’s energy.
As the quantum signature here scales as n̄dr  1, measuring
these low-frequency energy fluctuations is an easier way ofdetecting quantum behavior than attempting to resolve theinstantaneous phonon number and individual quantum jumps.
1050-2947/2011/84(4)/043824(10)043824-1©2011 American Physical Society

A. A. CLERKPHYSICAL REVIEW A 84, 043824 (2011)

the present problem, the full counting statistics of electronic
charge fluctuations in a chaotic quantum dot, were studied
in Ref. [17]. The (non-QND) photodetection statistics ofintracavity photons (as would be measured by a photodetector
in a cavity) was addressed in Refs. [18,19].
II. MODEL AND CALCULATION
A. Statement of the problem

Our damped, driven harmonic mode is described by the
HamiltonianĤ ≡ Ĥ0 + Ĥγ = h̄ωM ĉ† ĉ − h̄f (eiωD t ĉ + H.c.) + Ĥγ . (2)Here, the first term describes the resonator (frequency ωM ,number operator n̂ = ĉ† ĉ), Ĥγ describes the damping (at a rateγ ) and heating of the oscillator by a thermal bath, and f is themagnitude of the coherent oscillator driving force (frequency
ωD = ωM + δ). We take Hγ to correspond to the standard
model of a linear coupling to an Ohmic oscillator bath. We
also define the dimensionless oscillator force susceptibility asχ=1.
1 + 4(δ/γ )2(3)

We will be interested in the statistics of the time-integrated
energy m̂ [cf. Eq. (1)] in the case where the oscillator has
equilibrated to both the driving force and thermal bath long
before the initial time t = 0. We will also focus exclusively on
the long-time limit (e.g., an integration time t which is longcompared to 1/γ ). The average and second moment of m̂ areeasily found by solving the Heisenberg-Langevin equations
for our system [20,21]. In the long-time limit, the averagem̂ ∼ [n̄dr + n̄th ]t, where n̄th denote the thermal number ofoscillator quanta (determined by the bath temperature) andn̄dr = (2f/γ )2 χ is the average number of quanta due to thedriving force.
For the second central moment, we find in the long-timelimit
(δ m̂)2  ∼ (δ m̂)2 dr + (δ m̂)2 th ,defined the variance as m̂2 , an expression which is naturally
symmetrized in terms of n̂ products; that is, t tm̂2  =dt1
dt2 n̂(t1 )n̂(t2 )
0
0
 t

1 t
=
dt1
dt2 {n̂(t1 ),n̂(t2 )}.
(7)
2 0
0
If we define higher moments in the same way (e.g., define
the j th moment to be m̂j ), they too would be naturally
symmetrized; one might expect that this is then a sensible
way to proceed.
Unfortunately, being sensible is not enough to guarantee
physical relevance: similar to the standard theory of photodetection [22], one must instead model the actual detection
scheme to properly understand how the measured moments
correspond to a given operator ordering. As we are interested
in nondestructive detection, the answer here will not be thenormal-ordering prescription used in photodetection. A similar
problem arises in the measurement of current fluctuationsin quantum coherent conductors; the answer emerging fromstudies of this question is the so-called Keldysh operator
ordering. This ordering appears naturally in a number ofidealized measurement setups [10–12]; it can also be given
an elegant motivation using a path-integral formulation of
the Keldysh field-theoretic technique [23,24]. For the second
moment, the ordering coincides with the simple definition in
Eq. (7); for higher moments, the ordering prescriptions haveno simple intuitive form (see, e.g., Ref. [9] for the explicit
form for the third moment).
C. Using an auxiliary qubit to obtain P(m)We give here a quick derivation of the Keldysh ordering anduse it to derive P (m) for our system. Following Ref. [11], we
consider an idealized method for measuring P (m), in which n̂
couples dispersively to the σ̂z operator of an auxiliary two-levelsystem (TLS) with a coupling strength k/2:

(4)

Ĥint =where(δ m̂)2 th =2n̄th (1 + n̄th )tγ

represents a purely thermal contribution, whereas8n̄dr χ t
1
(δ m̂)2 dr =n̄th +
γ
2

(5)

(6)

represents extra energy fluctuations due to the driving force.
The last term in Eq. (6) here survives in the limit ofzero temperature (i.e., n̄th → 0) and is a quantum effect: itcorresponds to the shot-noise fluctuations arising from thediscreteness of the resonator’s energy.
B. Higher moments and the Keldysh ordering

Before calculating higher moments and the full distributionof m̂, we must pause to consider the operator-ordering ambiguity arising from the noncommutativity of n̂(t) at different
times. In calculating the second moment, we have naivelyh̄kn̂σ̂z .
2

(8)

As there are no other terms involving the TLS in the
Hamiltonian, we see that it simply experiences a magnetic
field ∝n̂. If n̂ were just a classical, time-dependent field n(t),then during the time interval
 t between 0 and t, the TLS would
precess an angle θ = k 0 n(t  )dt  = km. If instead m was aclassically stochastic variable described by the distribution
P (m), then the average of e−iθ over this distribution (at a fixed
coupling k) would be

−iθ
(9)
e k = dmP (m)e−ikm ≡ [k].
Thus, when viewed as a function of k, the average ofthe precession phase directly yields the moment generating
function [k] of the distribution P (m).
The above correspondence now provides a means fordefining P (m) in the quantum case [11]: we simply use thefact that the average on the left-hand side (LHS) of Eq. (9)corresponds to ρ↑↓ (t)/ρ↑↓ (0), where ρ↑↓ (t) is an off-diagonal
matrix element of the TLS’s reduced density matrix. We can043824-2

FULL COUNTING STATISTICS OF ENERGY . . .

PHYSICAL REVIEW A 84, 043824 (2011)

thus define the moment generating function [k] [and hence
P (m)] in the quantum case via:ρ↑↓ (t) [k] ≡
 = Trsys [Û (t; k)ρ̂sys (Û (t; −k))† ], (10)ρ↑↓ (0) kwhere the time evolution operator Û is defined as:

 
i t
h̄k
Ĥ (t  ) + n̂(t  )Û (t; k) = T exp −
h̄ 0
2

.

(11)

Here, Ĥ is given in Eq. (2), ρ̂sys is the initial measured
system (i.e., cavity plus bath) density matrix, and the trace
is taken over all system degrees of freedom. Furthermore,
the symbol T denotes time ordering. Equation (10) uniquely
specifies the operating ordering to use for each momentof P (m); this is the Keldysh ordering. We stress that thesame ordering emerges in the analysis of other idealizedmeasurement setups [12]; we also show in Appendix A that it
applies to a realistic setup where n̂ is coupled dispersively to a
detector cavity whose frequency is monitored using homodyne
detection.
In our case, the above scheme not only motivates the
Keldysh ordering, it also gives us a convenient way to calculate
the generating function [k]. The reduced density matrix ρ̂describing the TLS and the driven resonator (i.e., only theresonator’s dissipative environment is traced out) obeys thefollowing standard master equation:iρ̂˙ = − [Ĥ0 + Ĥint ,ρ̂] + γ (n̄th +1)D[ĉ]ρ̂ + γ n̄th D[ĉ† ]ρ,h̄(12)where Ĥint is the dispersive qubit-oscillator interaction given
in Eq. (8), Ĥ0 is defined in Eq. (2), and where for any operatorÂ we define D[Â]ρ̂ = Âρ̂ Â† − (Â† Âρ̂ + ρ̂ Â† Â)/2. We areusing the “quantum optics” version of the master equation,
which is appropriate for the high-Q limit we consider.
As shown in Ref. [25], using standard phase spacetechniques, one can solve Eq. (12) and thus directly obtain
[k]. Reference [25] used this quantity to study dephasing
and coherence revivals of the TLS; the emphasis was onunderstanding ρ↑↓ (t) for a fixed value of coupling k. In
contrast, our focus here is on how ρ↑↓ (t) behaves as a functionof k in the long-time limit, because it is this behavior which
will determine P (m) in the long-time limit.
D. Generating function for P(m)Using the procedure described above, and taking the longtime limit, the final result for the moment generating function
[k] has the simple form (cf. Eq. (22) in Ref. [25]):
[k] = dr [k]th [k],

(13)





γt
th [k] = exp −2(14b)
In taking the long-time limit, we have simply dropped termsin [k] which decay exponentially in time as exp(−γ t/2) or
faster.
Recalling that the j th cumulant of m is given by
dj
i j dkj ln [k]|k=0 , one can easily check that Eqs. (14a)and (14b) yield the same second and third cumulants obtained
from the Heisenberg-Langevin approach. We see that thepurely thermal fluctuations described by th are independentof the additional drive-induced fluctuations described by dr ;further note that these purely thermal fluctuations vanish in
the limit of zero temperature. In the remainder of the paper we
focus on the case n̄dr  n̄th ,1 and thus focus attention on the
distribution P (m) generated by dr [k].
III. ENERGY FLUCTUATION STATISTICS IN THE
CLASSICAL LIMIT

To gain some intuition, it is useful to first consider the driven
energy counting statistics in the classical, high-temperature
limit. Formally, one transforms the distribution describedby Eq. (14a) t to a distribution describing the time-integrated
energy s = 0 dt  E(t  ) = h̄ωM m, One can then rigorously takethe h̄ → 0 limit. Transforming back to our original variable m,
one finds
−ik n̄dr t.
(15)
dr,cl [k] = exp
1 + 4iχ n̄th (k/γ )In the long-time limit of interest, the corresponding probability
distribution function can be found within a saddle-point
approximation, yielding√ √1−( m − n̄dr t)2
Pcl (m) 
,
(16)
exp
2σcl2
8π mσcl2
with
σcl2 = 2n̄th χ /γ .

(14a)

(17)

We see that, classically, m is well approximated as√being thesquare of a Gaussian random variable with mean n̄dr t andstandard deviation σcl . This of course implies that, even in theclassical limit, m is not itself a Gaussian variable.
The above behavior is easily understood. Writing thecomplexcavity amplitude a(t) in terms of its mean value√n̄dr e−iωd t and a thermally fluctuating part δa(t), we have
 t
2
√



m(t) =
dt   n̄dr + eiωd t δa(t  )0

where th [k] describes a purely thermal (drive-independent)
contribution, and dr [k] describes additional fluctuations
related to the drive. One finds−ik n̄dr t
dr [k] = exp,
1 + 4iχ n̄th + 12 (k/γ ) − χ (k/γ )2




k2
1 k
− 2 −1 .
1 + 4i n̄th +
2 γ
γ

=

 t
0



√
δX(t) + iδY (t) 2
dt  n̄dr +
√
 .
2


(18)

In the second line, we have written the fluctuation δa(t)in terms of real-valued quadratures δX(t) and δY (t). In thelarge n̄dr limit, only the low-frequency part of the “intensityquadrature” noise δX(t) will be enhanced by n̄dr . To a good043824-3

A. A. CLERKPHYSICAL REVIEW A 84, 043824 (2011)

approximation, we may thus drop the δY (t) contribution andreplace δX(t) by its low-frequency part. We thus have

 t
2
√
1
n̄dr t + √dt  δX(t  )
m(t)
2t 0
√
≡ [ n̄dr t + δX(t)]2 .
(19)
One can easily confirm that the time-averaged intensityquadrature noise δX(t) defined above is a Gaussian randomvariable with zero mean and variance 2n̄th χ /γ ; Eq. (19) isthus in agreement with Eq. (16). It is worth noting that a fullclassical calculation of P (m) (including purely thermal effects)
yields an answer in complete agreement with the classicallimit of the moment generating functions given in Eqs. (14a)and (14b).
IV. ENERGY FLUCTUATION STATISTICS IN THE
QUANTUM REGIME
A. Basic results

It is tempting to make a simple extrapolation of the classicalenergy statistics to the quantum regime. Again, in the large n̄dr
limit it is the amplification of the thermal intensity-quadraturefluctuations δX(t) [cf. Eq. (19)] which determine P (m); one
might expect that the only difference in the quantum case is thatthese quadrature fluctuations are driven by both thermal andzero-point force noise. We would thus expect the distribution
to again be given by Eq. (16), with the simple modificationthat the variance σcl in Eq. (17) should be increased to include
zero-point fluctuations via the substitution n̄th → n̄th + 1/2.
However, as already mentioned in the introduction, thisis not the case. Instead, the full quantum moment generating
function in Eq. (14a) is related to the classical one [cf. Eq. (15)]by the simple substitution
dr [k; n̄th ] = dr,cl [k; n̄th → n̄eff [k]],k1n̄eff [k] = n̄th + + i .
2
4γ

(20a)
(20b)

Thus, one shifts n̄th both by the constant 1/2 (reflecting theinclusion of zero-point force noise), as well as by an imaginary
k-dependent term.
The k dependence of n̄eff implies nontrivial quantum
corrections to the third cumulant and higher involving a kind of
feedback, whereby higher-order cumulants depend on the formof lower-order cumulants. In the case of the third cumulant
m3 , one finds

1
d
3
3
,
(21)
m  = m cl − 3
m2 cl
d n̄th4γ
where we use mj cl to denote the naive expectation forthe j th cumulant: the j th classical cumulant obtained from
Eq. (15), with the substitution n̄th → n̄th + 1/2.
Heuristically, this feedback of lower moments into higher
moments is analogous to the situation in so-called “cascaded”
Langevin approaches [13,14]. One could heuristically obtainthe feedback term in Eq. (21) using such an approach, startingwith the assumption that the effective thermal occupation n̄th
in the classical distribution fluctuates in a way that is driven by(and hence correlated with) δm. Assuming that the fluctuationsof δ n̄th are slow compared to those of n̂ (allowing a two-stepaveraging procedure), one obtains
dm3  = m3 cl + 3
m2 cl δ n̄th (t)δm(t). (22)d n̄thThis recovers Eq. (21) if we take δ n̄th (t)δm(t) = −1/4γ
While usually derived in a heuristic fashion, cascaded
Langevin approaches have been used successfully to understand higher cumulants of current fluctuations in electronic
conductors [13,14]. Here, we stress that this picture emerges
directly from a fully quantum calculation.
Turning to the explicit form of the third cumulant, evaluating Eq. (21) yields
n̄dr t6
m3  = 2 χ 2 24(1 + 2n̄th )2 −
.
(23)
γ
χ
The first term here is just the expected classical answer: it
is always positive. The second term here is the nontrivial“feedback” quantum correction; as it involves the second
moment of the classical distribution, it is lower order in thesusceptibility χ than the first term. As a result, this correction
can make the skewness negative for a sufficiently small χ ,something that is impossible classically. Furthermore, in thelimit of a strongly detuned drive (i.e., χ  1), the quantum
skewness has a much larger magnitude (by a factor 1/χ ) thanthe corresponding classical answer. The quantum “feedback”
corrections similarly enhance all higher moments over thecorresponding classical answer in the large-detuning limit.
B. Negative probabilities at large drive detuning

The enhanced role of the nontrivial quantum corrections(arising from the k dependence of n̄eff [k]) in the large-detuning
limit |δ|  γ is even more apparent if one looks at the form of
the full distribution P (m). One finds that, at sufficiently low
temperature and large detuning, these quantum corrections
lead to P (m) becoming nonpositive definite (see Fig. 1). Thiscan be demonstrated analytically by just using the first four
cumulants of P (m). The solution of the Hamburger momentproblem [26] is a necessary and sufficient set of conditions for aset of moments to correspond to a positive-definite probability
distribution. Letting Cj to denote the j th cumulant scaled by
the variance [e.g., C3 = m3 /(m2 )3/2 ], the lowest-order
Hamburger positivity constraint involving the third moment is(C3 )2  C4 + 2.

(24)

In the long-time limit and for large detunings, C3 ∝ χ −1/2
while C4 is independent of χ . The above constraint is thus
violated by P (m) for a sufficiently large detuning; in the large-t
limit, the condition for violation becomes(δ/γ )2 > 89 (n̄dr tγ )(1 + 2n̄th )3 .

(25)

Heuristically, the quantum “feedback” contribution to the third
moment [second term in Eq. (21)] has made it too large for
P (m) to be positive definite.1

1

We stress that Eq. (24) is a sufficient but not necessary condition
for P (m) to exhibit negativity; higher-order Hamburger constraints
are violated at even smaller values of |δ|.

043824-4

FULL COUNTING STATISTICS OF ENERGY . . .

PHYSICAL REVIEW A 84, 043824 (2011)
V. INTERPRETATION

P m, t

The interpretation of negative quasiprobabilities of the kindfound here (negative “counting statistics”) was first given by
Nazarov and Kindermann [12]. We begin by quickly summarizing their findings, and then extend their interpretation
to argue that negative counting statistics correspond to thesame sort of nonclassical behavior tested by Leggett-Garg
inequalities; in particular, they imply that a “macrorealistic”
classical picture of the fluctuations of n(t) is not possible. At
the same time we comment on important differences between
negative counting statistics and violations of the Leggett-Garg
inequalities; the two are not necessarily equivalent.

0

0

4

2

0

2

m m

4

6

8

10

m2

FIG. 1. (Color online) Distribution P (m) of integrated energy
fluctuations m of a driven single-mode resonator, evaluated at a timet = 10/γ a drive detuning of δ = 10γ , and a driving force strength
which yields an average number of cavity quanta n̄dr = 5. The various
curves correspond to different resonator temperatures: n̄th = 3 (black
wide-dashed), n̄th = 0.25 (red, small-dashed), and n̄th = 0 (solidblue). As discussed in the text, for large drive detunings and low
temperatures, the distribution P (m) fails to be positive definite.

In Fig. 1, we plot the distribution P (m) as obtained from
Eq. (14a) for a large detuning δ = 10γ and for various bathtemperatures; as the temperature is lowered, the distribution fails to be positive for values of m > m. Figure 2plots the integrated negativity of the distribution, N[P ] =− dmP (m)θ (−P (m)) as a function of temperature. One
clearly sees that increasing the drive detuning causes thenegativity to emerge at progressively higher temperatures.
In the large-detuning limit of interest, one can show
analytically (see Appendix B) that it is indeed the anomalously
large magnitude of the third cumulant [second term in Eq. (23)]which causes the distribution to become nonpositive definite.
As shown in Appendix B, the large value of the third
cumulant causes the distribution to have the form of an Airyfunction convolved with a Gaussian; the oscillations of the
Airy function cause the P (m) to drop below 0 for values ofm  m.

0.1

NP

0.01
0.001

10 4
10 5
10 6
10 7
0.0

0.5

1.0

1.5

nth
FIG. 2. (Color online) Integrated total negativity of the distribution P (m) as a function of temperature n̄th for a time t = 10/γ anda driving force amplitude which yields n̄dr = 5. From left to right,
the three curves correspond to drive detunings of δ = 3γ , δ = 5γ ,and δ = 7γ . One sees that increasing the magnitude of the detuningcauses the negativity to persist to higher temperatures.

A. Interpretation of Keldysh-ordered quasiprobabilities as
describing “intrinsic” fluctuations

Nazarov and Kindermann considered an alternate idealized
setup for measuring m̂, where the detector is an infinitely heavy
mass, and the quantity to be measured (in our case n̂) is linearly
coupled to the position x̂ of the detector:
Hint = −Ax̂ n̂.

(26)

The detector-oscillator interaction is turned on for a time t.
Classically, the interaction would simply shift the detector
momentum an amount Ant = Am, while (due to its infinitemass) its position would be unchanged. Correspondingly, one
would expect that, in the quantum case, the final momentum
distribution of the detector mass would be a convolution of
the detector’s initial momentum distribution and the desired
probability distribution P (m). The only additional complication is that there may be a backaction effect: the distribution
P (m) may itself (via Hint ) depend on the value of the detector
mass position x. In our case, the interaction Hamiltonianin Eq. (26) implies that different values of x correspondto different values of the resonator frequency (and hence
drive detuning δ). As a result, different values of x will leadto different distributions P (m; δ = δ0 + Ax). Including thiseffect, one would then expect the following relation between
the detector mass Wigner function W (x,p; t) before and after
the interaction:W (x,p; t) = dmP (m; δ = δ0 +Ax)W (x,p − Am; 0). (27)This relation was rigorously derived in Ref. [12], with P (m,δ)
being the usual Keldysh-ordered distribution we have been
considering.
Consider the case where the detector mass is initially in a
Gaussian state with zero means, a momentum variance σp =Aσimp and a position variance σx = h̄σBA /A. We could thenuse the final momentum distribution of the mass to infer the
distribution P (m) to bePmeas (m; δ) ≡ A dxW [x,m/A; t]


1

= 
dm
dδ  P (m ; δ  )4π 2 σimp σBA(m − m )2(δ − δ  )2.
× exp −
exp
−
2
2
2σimp
2σBA

043824-5(28)

A. A. CLERKPHYSICAL REVIEW A 84, 043824 (2011)

The above equation provides us with an unambiguous way
to interpret P (m). It tells us that the Keldysh-ordered P (m,δ)
should be regarded as the “intrinsic” distribution describing thefluctuations of m for a given fixed value of δ. In contrast, themeasured distribution Pmeas (m,δ) is this intrinsic distribution
corrupted by the addition of measurement uncertainty. Thereis both an imprecision uncertainty σimp in m coming from themomentum uncertainty in the detector mass initial state, and abackaction uncertainty σBA in δ coming from the detector mass
position uncertainty. The Heisenberg uncertainty principle
applied to the detector mass implies that σimp σBA  1/2:
one can never eliminate both these sources of measurement
uncertainty. Nonetheless, Eq. (28) gives us a way to define theunderlying, measurement-noise-free distribution P (m).
Turning to the issue of positivity, we must have that themeasured distribution Pmeas (m,δ) is positive definite, as it isjust the final momentum distribution of the test mass. In thespecial case where P (m) is independent of x (i.e., a true QND
measurement where there is no backaction), this constraint
immediately implies that P (m) be positive definite. However,in the case relevant here, where backaction is important (i.e.,
different values of x affect the fluctuations of m), this is nolonger required: P (m) can exhibit negativity in these cases
without violating the positivity of Pmeas (m). This is precisely
what we we find in the driven number-fluctuation statistics at
low temperature and large drive detuning.
For our driven cavity, P (m,δ) will only change appreciably
when δ is varied an amount ∼γ . To avoid a sizable backaction,
we would thus want σBA  γ . The Heisenberg uncertainty
principle then implies σimp  1/γ , which sets a limit to thescale m of any negative regions in P (m,δ). These constraints
are indeed obeyed by our results.
Finally, note that, in the limit of a very weak measurement(i.e., σimp → ∞), the detector backaction σBA may tend tozero without violating the uncertainty principle, and hence themeasured distribution will just be a convolution of imprecision
noise with the Keldysh-ordered quasiprobability P [m,δ].
As first noted by Bednorz and Belzig [27], one can thus
relate Keldysh-ordered quasiprobabilities to the concepts of“weak measurement” and “weak values” as developed by
Vaidman and coworkers [28,29]. As discussed in Ref. [27],this correspondence allows one to link nonpositive definite
Keldysh quasiprobabilty distributions to the possibility ofmeasuring a nonclassical weak value.

B. Significance of negative probabilities

With Eq. (28) in hand, we can now view the failureof P (m) to be positive definite as a clear manifestation of
nonclassical behavior in our driven resonator. Classically,
we would naturally think of the “intrinsic” fluctuations ofm (i.e., free from backaction effects) in terms of random
trajectories n(t) and a corresponding distribution function.
The failure of the Keldysh-ordered P (m) to be positive means
that, even in the most highly idealized measurement setups,
we cannot interpret the measured statistics in this way: they
do not correspond to having added measurement noise to anunderlying classical stochastic process. This interpretation can
only be salvaged if one is willing to accept that the intrinsic

distribution function describing the fluctuations is not positive
definite.
Not surprisingly, negativity in P (m) and the corresponding nonclassical behavior only emerges at sufficiently low
temperatures. More subtle, however, is the dependence on
drive detuning: negativity only occurs at a sufficiently largemagnitude detuning |δ| and is always enhanced by increasing|δ|. On a basic level, this is consistent with Eq. (28), which tells
us that P (m) can only be negative if there is a backaction effectassociated with measuring m. In our case, this backaction
effect vanishes to leading order when δ = 0 [9]. Hence, it
is reasonable that obtaining negativity requires a nonzero
detuning δ.
We stress that this nonclassicality discussed here is very different than that associated with a nonpositive definite Wigner
function; here, we are not characterizing the instantaneousstate of a system, but rather the time-integrated fluctuationsof an observable. There is, rather, a much more naturalconnection to the kind of nonclassical temporal correlations
that lead to violations of Leggett-Garg inequalities (LGIs) [16].
These inequalities constrain temporal correlations of a given
observable O(t) in any classical theory which satisfies(1) Macrorealism: O(t) has a definite value at all times.
(2) Noninvasive measurability: O(t) can be measured without any backaction disturbance that would alter its subsequentevolution.
We note that two recent experiments have reported violationof an LGI [30,31].
As we have discussed, a nonpositive definite P (m) also
is indicative of a failure of macrorealism; the “intrinsic”
fluctuations of n(t) cannot be thought of as a classicalstochastic process. Furthermore, the standard violation of an
LGI involves a qubit undergoing Larmor oscillations [16].
Similarly, in our driven resonator negative probabilities only
emerge for large-magnitude drive detunings |δ|  γ , a regime
where correlation functions of n̂(t) have a strong oscillatory
behavior.
Despite these connections, there are important differencesfrom the violation of an LGI. In some sense, the nonclassicality
associated with the negativity of P (m) is stronger than thatassociated with the violation of an LGI. Violating an LGI could
simply be interpreted as indicating that the measurement is
indeed invasive (i.e., there is backaction), without necessarily
indicating a failure of macrorealism. In contrast, negativity ofP (m) tells us more than simply “backaction exists.” It tells usthat a natural way of including backaction effects classically(as additional measurement noise which smears an intrinsic
probability distribution, [cf. Eq. (28)]) is impossible. Note
that, in our system, backaction effects remain important whenmeasuring P (m) at nonzero detuning even in the more classicalregime of high temperatures; nonetheless, the distribution
exhibits no negativity here.
Yet another difference is that while standard LGIs involve
two-time correlation functions, the nonclassical behavior
found here only manifests itself when one considers higherorder correlation functions. Recall that the second moment
given in Eq. (7) has a completely classical form, where
the third moment and higher have nontrivial quantum corrections stemming from the k dependence of the effective
thermal occupancy factor n̄eff [k]. This is in keeping with the043824-6

FULL COUNTING STATISTICS OF ENERGY . . .

PHYSICAL REVIEW A 84, 043824 (2011)conclusions of Ref. [32], which argues that, for general observables (i.e., continuous, unbounded), having access only to
two-point correlation functions is not sufficient to demonstratenonclassical behavior.
In closing our discussion here, we note that one could
attempt to derive a generalized LGI involving correlations
between n(t) evaluated at a finite number of measurement
times. As discussed above, such an inequality would necessarily involve more that two-point correlation functions. A similar
generalized Bell inequality (involving fourth moments) was
derived in Ref. [32]. Evaluating such an inequality would
necessarily involve a knowledge of time-resolved phononcounting statistics (similar to time-resolved electron counting
statistics [33]), and would involve going beyond the long-time
limit considered here. While these issues can be addressedusing a similar calculation technique as in Sec. II C, it
introduces considerable additional complexity. We will thus
address this topic in a separate future work.
VI. MEASUREMENT ISSUES

We end with a discussion of how one might experimentally
detect evidence of the nonclassical photon and phonon number
fluctuations described in this paper.

To resolve this nonlinear contribution to the ac Stark shift,
one requires sufficient qubit coherence. Because we must allow
the qubit phase to evolve long enough both to be in the longtime limit of the photon-number fluctuations and to resolve thek 3 stark shift, we need the total dephasing rate (including thecontribution from T1 processes) to satisfy


k3
n̄
<
min
γ
,
χ
.
(30)
ϕ
dr
γ2
One unavoidable contribution to ϕ will come from thedispersive coupling and the even moments of P (m); it iseasy to see that this contribution satisfies Eq. (30) as long as
one has a strong dispersive coupling k  γ . This regime has
been achieved in several recent circuit QED experiments [5,6].
One also needs the intrinsic, noncavity dephasing of the qubit
to be sufficiently small. Given recent advances in extending
the coherence of superconducting qubits [34], this also would
appear to be feasible.
An alternate approach would be to measure the order k 3
term in the ac Stark shift via simple spectroscopy, where one
directly drives the qubit and measures its state as a function ofthis drive frequency. In order to avoid complications arisingfrom the spectroscopy drive modifying the cavity state [andhence P (m)], one would need to use, for example, a second
cavity for the spectroscopy [5].

A. Qubit plus resonator measurement

One approach would be to experimentally implement themodel of Sec. II C, where a qubit is coupled dispersively tothe photon number operator n̂ of a driven cavity [cf. Eq. (8)].
The evolution of the qubit phase (i.e., its off-diagonal density
matrix element ρ↑↓ ) at various values of the dispersive coupling k directly yields the cumulant generating function of thedesired distribution P (m) [cf. Eq. (10)]. Such a measurement
could be contemplated in a cavity QED or circuit QED
systems, where a two-level system (atom or superconducting
qubit) is coupled to an electromagnetic cavity. The required
phase information could be extracted by using a Ramseyinterference technique, similar to the seminal experiments ofRefs. [3,4]. Unlike those experiments, the focus here is very
different: the goal is to learn about the way the cavity photon
number fluctuates over a time t  1/γ , as opposed to probing
its instantaneous value.
In order to detect evidence of the nonclassical photon
number fluctuations discussed here, it would be sufficient
to see that the third moment m3  is negative. Given thedispersive qubit-cavity interaction, the odd moments of P (m)will contribute in the long-time limit t  1/γ to the ac Starkshift
qb of the qubit frequency. One hasm k 3 m3 + O(k 5 )−qb = lim k
t→∞
t
6
t
= k n̄dr +k3n̄dr χ + O(k 5 χ 3 ).
γ2

(29)

In the second line, we have assumed that the driven cavity is in
the interesting regime of zero temperature and strong detuning(i.e., χ  1), where we expect strong quantum effects. We seethat anomalous negative skewness of P (m) manifests itself
in the sign of the k 3 term in the ac Stark shift of the qubit
frequency.

B. Measurement with a generic linear-response detector

Perhaps a more general way to measure the fluctuations ofm would be to weakly couple the photon (or phonon) number
operator n̂ we wish to measure to the input port of a genericlinear detector [21], as discussed in Ref. [9]. We would thus
have a coupling of the form
Hint,lin = h̄n̂ · F̂(31)where F̂ is a detector operator and a generalized force.
We would then monitor some other detector observable,
say Iˆ, whose value depends linearly on n̂. The dispersively
coupled optomechanical setup for detecting phonon shot noise
analyzed in Ref. [9] falls into this general description. In thatcase, F̂ is the photon number operator of the optical cavity used
to detect mechanical quanta, and Iˆ is the homodyne current.
In this sort of generic setup, the statistics of the detector
output I can be used to to extract the statistics of m. Of
course (similar to the idealized detector of Sec. V A), thiscorrespondence will be hindered by the presence of bothmeasurement imprecision noise (i.e., the intrinsic fluctuationsin Iˆ) as well as backaction noise (i.e., the effective fluctuations
in detuning resulting from fluctuations in F̂ ). The simplest
evidence for nonclassical fluctuations of m comes fromthe anomalous sign of the third moment; we thus need to
ask whether measurement imprecision and backaction would
obscure the intrinsic negativity of the skewness.
To that end, we first note that measurement imprecision
here can be treated as an additive Gaussian noise process, andhence will not change the value of the third moment. As forthe backaction fluctuations, they are equivalent to having phasefluctuations on the mechanical drive. Treating these backaction
phase fluctuations along the same lines as Ref. [35], we findthat they yield an additional additive contribution to m3 

043824-7

A. A. CLERKPHYSICAL REVIEW A 84, 043824 (2011)

which is always positive, and which in the large n̄dr limit takes
the form
SF F n̄dr 2 96n̄dr t [1 − 12(δ/γ )2 ]2. (32)
m3 BA =γγ 2 (1 + 4δ 2 /γ 2 )4
By choosing a sufficiently small measurement strength (i.e.,
SF F ) and sufficiently large detuning, one can still have thetotal third moment be negative. In the large detuning limit,
the backaction-induced skewness scales as 1/δ 4 , whereas theintrinsic, negative skewness scales as 1/δ 2 . Note also thatthe backaction contribution to the second moment m2  forthis generic linear-detector setup was discussed in Ref. [9], inwhich one finds that the minimum possible total added noiseis achieved for SF F ∝ 1/n̄dr .
VII. CONCLUSIONSIn this paper, we have shown that the full-counting statistics
of energy fluctuations in a driven quantum resonator canbecome negative for sufficiently low temperature and large
drive detuning. This negativity arises from the same kind ofquantum correction that leads to a negative third moment [9],something that is impossible classically. We have argued
that the failure of the quasiprobability distribution describing
P (m) to be positive definite is similar to having violated a
Leggett-Garg inequality and implies that a “macrorealistic,”
classical picture for the fluctuations of phonon or photon
number is not possible. We have thus shown in a relatively
simple setting that higher moments of such counting statistics
can be used to detect nonclassical behavior. It would be
extremely interesting to investigate whether similar effects
manifest themselves in other system.
ACKNOWLEDGMENTS

We thank T. Bhattacharya for a useful discussion. This
work was supported by the Canadian Institute for Advanced
Research and by the DARPA ORCHID program through agrant from AFOSR.
APPENDIX A : KELDYSH ORDERING FROM HOMODYNE
MEASUREMENT THEORY

Consider the measurement setup of recent optomechanicalexperiments [7,8], where the number operator n̂(t) of amechanical mode is coupled dispersively (strength g) to adriven optical mode which acts as a detector cavity. The
coupling takes the form
Hint = g n̂â † â,(A1)

where â is the annihilation operator for the measurement
cavity. By virtue of this interaction, the frequency of thedetector cavity will depend on the value of n(t). One can thus
measure the time variation of n(t) by detecting the resulting
variation in the detector cavity frequency via homodyne
detection of its output field (see, e.g., [21,36]). This involves
first mixing the detector cavity output field with a large,
classical reference beam. To leading order in g, the output
field b̂ from the mixer will have the formb̂(t) = β + B n̂(t),where β parametrizes the large magnitude of the classicalreference beam used, B ∝ g, and we have omitted vacuumnoise terms responsible for shot noise (i.e., the imprecision
noise in this measurement scheme). Without loss of generality,
we take both β and B to be real. The intensity Iˆ = b̂† b̂
of the mixer output is then measured using a photodetector.
Assuming that the constant term β 2 intensity can be subtracted
from Iˆ (by, e.g., using balanced homodyne detection), theoutput of the detector to leading order in β is
δ Iˆ = Iˆ − β 22Bβ n̂.

(A3)

We see that, in the large β limit, the output is just linearly
proportional to n̂.
Given this simple linear correspondence, one can directly
infer the value of the j th moment of m from the measured
j th-order intensity correlation function; that is,

1
β→∞ (2βB)jmj meas ≡ lim

 t
j
0 l=1

dtl T̃ [δ Iˆ(t1 ) · · · Iˆ(tj )]. (A4)On the right-hand side (RHS) of this expression, the symbolT̃ denotes that the measured intensity correlation functions
correspond to expectation values which are both normal
ordered and time ordered with respect to the b(t) and b† (t)operators; this ordering prescription is a direct consequence ofmeasuring intensity via photodetection [22,37]. For example,
for the third moment we haveT̃ I (t1 )I (t2 )I (t3 ) = b̂† (ta )b̂† (tb )b̂† (tc )b̂(tc )b̂(tb )b̂(ta ), (A5)where ta < tb < tc denotes the time-ordered listing of t1 ,t2 ,and t3 .
Equations (A2) and (A4) completely determine the correspondence between the measured moments mj meas andappropriately ordered expectation values of products of m̂(t).
Although tedious, one can now explicitly confirm that, for each
moment mj meas , the resulting ordering of m̂(t) operators isexactly the Keldysh ordering defined by Eq. (10). For the third
moment, this was done in Ref. [9].
A more compact way to see that one obtains the Keldysh
ordering at each order is to use a functional-integral formulation of the Keldysh technique; a pedagogical introduction tothis approach is given in Refs. [23,24]. In this formulation,
each bosonic operators is replaced by two time-dependent
fields, for example, b̂(t) → bσ (t),n̂(t) → nσ (t), where theindex σ = + (−) denotes the forward (backwards) Keldysh
contour. Within this approach, different operator orderings
correspond to different combinations of + and − fields. A
special role is played by the so-called “classical” field, whichis the average of + and − fields; that is,(A2)
043824-8

ncl (t) ≡n+ (t) + n− (t)
.
2

(A6)

FULL COUNTING STATISTICS OF ENERGY . . .

PHYSICAL REVIEW A 84, 043824 (2011)

At the level of a saddle-point approximation, the dynamics
of the classical field correspond to an effective classicalLangevin equation. Correlation functions of this classical field
are obtained in the usual way using the Keldysh action S
describing the system: ncl (t1 ) · · · ncl (tj ) ≡[Dφj σ ]ncl (t1 ) · · · ncl (tj ) exp[iS].
j,σ =±

˜ dr [k] of the transformed distribution takes the form
CGF 
⎞
⎛
2
1 + i 4√χ1n̄dr tγ kk
˜ dr [k] = − ⎝
⎠.

(B1)

χ
1
2 1+i
2
k
−
k
4n̄dr tγn̄dr tγ
Consider now the strong-detuning, long-time limit, suchthat n̄dr γ t → ∞ but χ n̄dr γ t is finite. In this limit
˜ dr [k] → − 1 k 2 − i √ 1k3.
2
8 χ n̄dr tγ(A7)
The φj σ (t) here denote the various fields that describe thesystem; the action is a function of these fields. By construction,
the j th-order correlation functions defined above is identicalto the j th-order, Keldysh-ordered operator expectation value
defined by Eq. (10) [23].
Turning to our homodyne measurement, we first notethat the normal-ordered, time-ordered intensity correlation
functions that are measured via photodetection can be obtainedby adding an auxiliary source term to the Keldysh action S of
the form
 t∗
Ssrc = kdt  [b−(t)b+ (t) − β 2 ].
(A8)0

Given the correspondence between the Keldysh ± fields andoperator orderings [24], one finds that derivatives of the full
Keldysh partition function (action S + Ssrc ) with respect to
k at k = 0 generate the desired normal and time-orderedcorrelation functions in the usual way.
Next, for homodyne detection, we can make the replacementbσ (t)bσ∗ (t)which results in
Ssrc

β + Bnσ (t),
 t

2kβB

(A9)

(B2)

In the quantum case, both the second and third moments
are nonvanishing in this long-time, strong-detuning limit. In
contrast, in the same limit the classical distribution would becompletely Gaussian. Thus, the third-moment term in Eq. (B2)is entirely due to the effective k dependence of the thermal
factor n̄eff in the quantum distribution. Furthermore, note that
increasing the drive detuning (and hence reducing χ ) enhances
the non-Gaussian nature of the distribution described by
Eq. (B2); this is the opposite of what happens classically, wherea large detuning suppresses fluctuations and non-Gaussian
effects, as the magnitude of thermal fluctuations at the drive
frequency are reduced.
Fourier transforming the approximate CGF in Eq. (B2)reveals that, in the large detuning limit, the distribution P (m)is a convolution of a Gaussian and an Airy function.This canbe explicitly evaluated. Defining m̃ = (m − n̄dr t)/ m2 ),we have





1
m̃
1
1
−1
m̃ − 3Ai − + 4 ,P (m̃,t)
exp
λ
2λ3
6λ
λ
4λ
(B3)
where

dt  ncl (t),λ3 = −(A10)

0

m3 
3
.
∼ √
3/2
2
8
χ
n̄dr tγ2m 

(B4)

We first shift and rescale P (m) in the full quantum case sothat it has zero mean and unit variance. Setting n̄th = 0, theIt is the oscillation of the Airy function factor in Eq. (B3)above which gives rise to the negative probabilities at m̃ > 0.
The exponential prefactor ensures that the resulting negativity
is exponentially suppressed in the long-time limit when λ  1.
However, for intermediate times (still much longer than 1/γ ),one has λ  1, and the negativity can be appreciable. Note thatthe most prominent domain of negativity in this large-λ limit
has an extent in m̃ ∼ λ; in terms of m/t, this corresponds toa range < 1. Thus, while P (m) exhibits negativity even in aseemingly classical regime where n̄dr γ t  1, it only occurs
on a scale which corresponds to less than one quantum in theresonator. This is consistent with the discussion of negative
quasiprobabilities given in Sec. V A.

[1] P. Grangier, J. A. Levenson, and J.-P. Poizat, Nature (London)
396, 537 (1998).
[2] G. Nogues, A. Rauschenbeutel, S. Osnaghi, M. Brune, J. M.
Raimond, and S. Haroche, Nature (London) 400, 239 (1999).
[3] C. Guerlin, J. Bernu, S. Deleglise, C. Sayrin, S. Gleyzes, S. Kuhr,
M. Brune, J.-M. Raimond, and S. Haroche, Nature (London)
448, 889 (2007).

[4] S. Gleyzes, S. Kuhr, C. Guerlin, J. Bernu, S. Deleglise,
U. B. Hoff, M. Brune, J.-M. Raimond, and S. Haroche, Nature(London) 446, 297 (2007).
[5] D. Schuster et al., Nature (London) 445, 515 (2007).
[6] B. R. Johnson et al., Nature Phys. 6, 663 (2010).
[7] J. D. Thompson, B. M. Zwickl, A. M. Jayich, F. Marquardt, S. M.
Girvin, and J. G. E. Harris, Nature (London) 452, 72 (2008).

where we have only retained the leading-order-in-β term in
the action. We thus see that the source field k couples tothe classical field ncl (t); it thus follows that the j th-orderintensity correlation functions (as generated by differentiationof the Keldysh partition function with respect to k) will bedirectly proportional to j th-order Keldysh-ordered correlation
functions of m.
APPENDIX B : P(m) IN THE LARGE TIME, LARGE
DETUNING LIMIT

043824-9

A. A. CLERKPHYSICAL REVIEW A 84, 043824 (2011)[8] J. C. Sankey, C. Yang, B. M. Zwickl, A. M. Jayich, and J. G. E.
Harris, Nature Phys. 6, 702 (2010).
[9] A. A. Clerk, F. Marquardt, and J. G. E. Harris, Phys. Rev. Lett.
104, 213603 (2010).
[10] L. Levitov and G. Lesovik, JETP Lett. 58, 230 (1993).
[11] L. Levitov, H. Lee, and G. Lesovik, J. Math. Phys. 37, 4845
(1996).
[12] Y. V. Nazarov and M. Kindermann, Eur. Phys. J. B 35, 413(2003).
[13] K. Nagaev, Phys. Rev. B 66, 104431 (2002).
[14] C. W. J. Beenakker, M. Kindermann, and Y. V. Nazarov, Phys.
Rev. Lett. 90, 176802 (2003).
[15] W. Belzig and Y. V. Nazarov, Phys. Rev. Lett. 87, 197006 (2001).
[16] A. J. Leggett and A. Garg, Phys. Rev. Lett. 54, 857 (1985).
[17] S. Pilgram and M. Büttiker, Phys. Rev. B 67, 235308 (2003).
[18] L. Mandel, Opt. Acta 28, 1447 (1981).
[19] M. D. Srinivas and E. B. Davies, Opt. Acta 29, 235 (1982).
[20] C. W. Gardiner and P. Zoller, Quantum Noise (Springer, Berlin,
2000).
[21] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and
R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010).
[22] R. J. Glauber, Phys. Rev. 130, 2529 (1963).
[23] A. Kamenev, in Nanophysics: Coherence and Transport, editedby H. Bouchiat et al. (Elsevier, Amsterdam, 2005), p. 177.

[24] A. Kamenev and A. Levchenko, Adv. Phys. 58, 197 (2009).
[25] A. A. Clerk and D. W. Utami, Phys. Rev. A 75, 042302(2007).
[26] D. V. Widder, The Laplace Transform (Princeton University
Press, Princeton, 1946).
[27] A. Bednorz and W. Belzig, Phys. Rev. Lett. 105, 106803(2010).
[28] Y. Aharonov, D. Z. Albert, and L. Vaidman, Phys. Rev. Lett. 60,
1351 (1988).
[29] Y. Aharonov and L. Vaidman, Phys. Rev. A 41, 11 (1990).
[30] A. Palacios-Laloy, F. Mallet, F. Nguyen, P. Bertet, D. Vion,
D. Esteve, and A. N. Korotkov, Nat. Phys. 6, 442 (2010).
[31] M. E. Goggin, M. P. Almeida, M. Barbieri, B. P. Lanyon, J. L.
O’Brien, A. G. White, and G. J. Pryde, Proc. Acad. Nat. Sci.
108, 1256 (2011).
[32] A. Bednorz and W. Belzig, Phys. Rev. B 83, 125304 (2011).
[33] A. Bednorz and W. Belzig, Phys. Rev. Lett. 101, 206803(2008).
[34] H. Paik et al., e-print arXiv:1105.4652v1 (to be published).
[35] P. Rabl, C. Genes, K. Hammerer, and M. Aspelmeyer, Phys.
Rev. A 80, 063819 (2009).
[36] D. F. Walls and G. J. Milburn, Quantum Optics, 2nd ed.
(Springer, Berlin, 2008).
[37] H. Carmichael, J. Opt. Soc. Am. B 4, 1588 (1987).

043824-10