. 21
( 27)


‚ env
ρW (t) = [HSE (t) , ρenv (t)]
i (18.6)
in the environment picture. The transformed operators,

Oenv (t) = UW 0 (t) OUW 0 (t) , (18.7)
satisfy the equations of motion
‚ env
O (t) = [Oenv (t) , HW 0 (t)] .
i (18.8)
env env env
One should keep in mind that HW 0 (t) = HS (t) + HE (t), and that the sample
Hamiltonian, HS (t), still contains all interaction terms between di¬erent degrees
of freedom in the sample. Only the sample“environment interaction is treated as a
perturbation. Thus the environment-picture sample operators obey the full Heisenberg
equations for the sample.

18.3 Averaging over the environment
In line with our general convention, we will now drop the identifying superscript ˜env™,
and replace it by the understanding that states and operators in the following discus-
sion are normally expressed in the environment picture. Exceptions to this rule will
¼ The master equation

be explicitly identi¬ed. Our immediate task is to derive an equation of motion for the
reduced density operator ρS . In pursuing this goal, we will generally follow Gardiner™s
treatment (Gardiner, 1991, Chap. 5).
The ¬rst part of the argument corresponds to the formal elimination of the reser-
voir operators in Section 14.1.2, and we begin in a similar way by incorporating the
quantum Liouville equation (18.6) and the initial density operator ρW (0) into the
equivalent integral equation,
ρW (t) = ρW (0) ’ dt1 [HSE (t1 ) , ρW (t1 )] . (18.9)

The assumption that HSE is weak”compared to HS and HE ”suggests solving
this equation by perturbation theory, but a perturbation expansion would only be
valid for a very short time. From Chapter 14, we know that typical sample correlation
functions decay exponentially:

Q1 (t + „ ) Q2 (t) ∼ e’γ„ , (18.10)

with a decay rate, γ ∼ g 2 , where g is the sample“environment coupling constant. This
exponential decay can not be recovered by an expansion of ρW (t) to any ¬nite order
in g.
As the ¬rst step toward ¬nding a better approach, we iterate the integral equation
(18.9) twice”this is suggested by the fact that γ is second order in g”to ¬nd
ρW (t) = ρW (0) ’ dt1 [HSE (t1 ) , ρW (0)]
2 t t1
+’ dt1 dt2 [HSE (t1 ) , [HSE (t2 ) , ρW (t2 )]] . (18.11)
0 0

Tracing over the environment then produces the exact equation
ρS (t) = ρS (0) ’ dt1 TrE ([HSE (t1 ) , ρW (0)])
2 t t1
+’ dt1 dt2 TrE ([HSE (t1 ) , [HSE (t2 ) , ρW (t2 )]]) (18.12)
0 0

for the reduced density operator. Since our objective is an equation of motion for
ρS (t), the next step is to di¬erentiate with respect to t to ¬nd

‚ i
ρS (t) = ’ TrE ([HSE (t) , ρW (0)])
’2 dt TrE ([HSE (t) , [HSE (t ) , ρW (t )]]) . (18.13)

This equation is exact, but it is useless as it stands, since the unknown world den-
sity operator ρW (t ) appears on the right side. Further progress depends on ¬nding
Averaging over the environment

approximations that will lead to a manageable equation for ρS (t) alone. The ¬rst sim-
plifying assumption is that the sample and the environment are initially uncorrelated:
ρW (0) = ρS (0) ρE (0). By combining the generic expression,

grν b† Qr ’ grν Q† brν ,

HSE = i (18.14)
rν r
r ν

for the sample“environment interaction with the conventional assumption, brν = 0,
and the initial factorization condition, it is straightforward to show that

TrE ([HSE (t) , ρW (0)]) = 0 . (18.15)

If one or more of the reservoirs has brν E = 0, one can still get this result by
writing HSE in terms of the ¬‚uctuation operators δbrν = brν ’ brν E , and absorbing
the extra terms by suitably rede¬ning HS and HE , as in Exercise 18.1.
Thus the initial factorization assumption always allows eqn (18.13) to be replaced
by the simpler form
‚ 1
ρS (t) = ’ 2 dt TrE {[HSE (t) , [HSE (t ) , ρW (t )]]} . (18.16)
‚t 0

Replacing ρW (t ) by ρW (0) = ρS (0) ρE (0) in eqn (18.16) would provide a perturba-
tive solution that is correct to second order, but”as we have just seen”this would
not correctly describe the asymptotic time dependence of the correlation functions.
The key to ¬nding a better approximation is to exploit the extreme asymmetry
between the sample and the environment. The environment is very much larger than
the sample; indeed, it includes the rest of the universe. It is therefore physically rea-
sonable to assume that the fractional change in the sample, caused by interaction
with the environment, is much larger than the fractional change in the environment,
caused by interaction with the sample. If this is the case, there will be no reciprocal
correlation between the sample and the environment, and the density operator ρW (t )
will be approximately factorizable at all times.
This argument suggests the ansatz

ρW (t ) ≈ ρS (t ) ρE (0) , (18.17)

and using this in eqn (18.16) produces the master equation:
‚ 1
ρS (t) = ’ 2 dt TrE {[HSE (t) , [HSE (t ) , ρS (t ) ρE (0)]]} . (18.18)
‚t 0

The double commutator in eqn (18.18) can be rewritten in a more convenient way
by exploiting the fact that typical interactions have the form

F (t) + F† (t) .
HSE (t) = (18.19)

This in turn allows the double commutator to be written as
¾ The master equation

C2 (t, t ) = [HSE (t) , [HSE (t ) , ρS (t ) ρE (0)]]

F† (t) , G (t ) + HC ,
= {[F (t) , G (t )] + HC} + (18.20)

G (t ) = [F (t ) , ρS (t ) ρE (0)] . (18.21)

18.4 Examples of the master equation
In order to go on, it is necessary to assume an explicit form for HSE . For this purpose,
we will consider the two concrete examples that were studied in Chapter 14: the single
cavity mode and the two-level atom.

18.4.1 Single cavity mode
In the environment picture, the de¬nition (14.43) of system“reservoir interaction for
a single cavity mode becomes

v („¦ν ) a† (t) bν (t) ’ HC .
HSE = i (18.22)

Since a (t) and bν (t) are evaluated in the environment picture, they satisfy the Heisen-
berg equations
‚ 1
a (t) = [a (t) , HS (t)]
‚t i
= ’iω0 a (t) + [a (t) , HS1 (t)] , (18.23)

bν (t) = ’i„¦ν bν (t) . (18.24)
By introducing the slowly-varying envelope operators a (t) = a (t) exp (iω0 t) and
bν (t) = bν (t) exp (iω0 t), we can express HSE (t) in the form (18.19), with

F (t) = ’iξ † (t) a (t) , (18.25)

where we recognize

v („¦ν ) bν (t0 ) e’i(„¦ν ’ω0 )(t’t0 )
ξ (t) = (18.26)

as the noise operator de¬ned in eqn (14.52).
The terms in [F (t) , G (t )] contain products of ξ † (t), ξ † (t ), and ρE (0) in various
orders. When the partial trace over the environment states in eqn (18.18) is performed,
the cyclic invariance of the trace can be exploited to show that all terms are propor-
tional to
ξ † (t) ξ † (t ) E = TrE ρE (0) ξ † (t) ξ † (t ) . (18.27)
Just as in Section 14.2, we will assume that ρE (0) is diagonal in the reservoir
oscillator occupation number”this amounts to assuming that ρE (0) is a stationary
Examples of the master equation

distribution”so that the correlation functions in eqn (18.27) all vanish. This assump-
tion is convenient, but it is not strictly necessary. A more general treatment”that
includes, for example, a reservoir described by a squeezed state”is given in Walls and
Milburn (1994, Sec. 6.1).
When ρE (0) is stationary, [F (t) , G (t )] and its adjoint will not contribute to the
master equation. By contrast, the commutator F† (t) , G (t ) is a sum of terms con-
taining products of ξ (t), ξ † (t ), and ρE (0) in various orders. In this case, the cyclic
invariance of the partial trace produces two kinds of terms, proportional respectively
ξ † (t) ξ (t ) E = ncav κδ (t ’ t ) (18.28)
ξ (t) ξ † (t ) = (ncav + 1) κδ (t ’ t ) , (18.29)
where ncav is the average number of reservoir oscillators at the cavity-mode frequency
ω0 .
The explicit expressions on the right sides of these equations come from eqns (14.74)
and (14.75), which were derived by using the Markov approximation. Thus the master
equation also depends on the Markov approximation, in particular on the assump-
tion that the envelope operator a (t) is essentially constant over the memory interval
(t ’ Tmem /2, t + Tmem /2).
After evaluating the partial trace of the double commutator C2 (t, t )”see Exercise
18.2”the environment-picture form of the master equation for the ¬eld is found to be
‚ κ
ρS (t) = ’ (ncav + 1) a† (t) a (t) ρS (t) + ρS (t) a† (t) a (t) ’ 2a (t) ρS (t) a† (t)
‚t 2
’ ncav a (t) a† (t) ρS (t) + ρS (t) a (t) a† (t) ’ 2a† (t) ρS (t) a (t) .
The slowly-varying envelope operators a (t) and a† (t) are paired in every term;
consequently, they can be replaced by the original environment-picture operators a (t)
and a† (t) without changing the form of the equation. The right side of the equation
of motion is therefore entirely expressed in terms of environment-picture operators, so
we can easily transform back to the Schr¨dinger picture to ¬nd

ρS (t) = LS ρS (t) + Ldis ρS (t) . (18.31)
The Liouville operators LS ”describing the free Hamiltonian evolution of the sam-
ple”and Ldis ”describing the dissipative e¬ects arising from coupling to the environ-
ment”are respectively given by
ω0 a† a + HS1 (t) , ρS (t)
LS ρS (t) = (18.32)
Ldis ρS (t) = ’ (ncav + 1) a† , aρS (t) + ρS (t) a† , a
2 (18.33)
’ ncav a, a† ρS (t) + ρS (t) a, a† .
The master equation

The operators we are used to, such as the Hamiltonian or the creation and annihi-
lation operators, send one Hilbert-space vector to another. By contrast, the Liouville
operators send one operator to another operator. For this reason they are sometimes
called super operators.

A Thermal equilibrium again
In Exercise 14.2 it is demonstrated that the average photon number asymptotically
approaches the Planck distribution. With the aid of the master equation, we can study
this limit in more detail. In this case, HS1 (t) = 0, so we can expect the density operator
to be diagonal in photon number. The diagonal matrix elements of eqn (18.31) in the
number-state basis yield

dpn (t)
= ’κ {(ncav + 1) n + ncav (n + 1)} pn (t)
dt (18.34)
+ κ (ncav + 1) (n + 1) pn+1 (t) + κncav n pn’1 (t) ,

where pn (t) = n |ρ (t)| n . The ¬rst term on the right represents the rate of ¬‚ow
of probability from the n-photon state to all other states, while the second and third
terms represent the ¬‚ow of probability into the n-photon state from the (n + 1)-photon
state and the (n ’ 1)-photon state respectively.
In order to study the approach to equilibrium, we write the equation as

dpn (t)
= Zn+1 (t) ’ Zn (t) , (18.35)
Zn (t) = nκ {(ncav + 1) pn (t) ’ ncav pn’1 (t)} . (18.36)
The equilibrium condition is Zn+1 (∞) = Zn (∞), but this is the same as Zn (∞) = 0,
since Z0 (t) ≡ 0. Thus equilibrium imposes the recursion relations

(ncav + 1) pn (∞) = ncav pn’1 (∞) . (18.37)

This is an example of the principle of detailed balance; the rate of probability
¬‚ow from the n-photon state to the (n ’ 1)-photon state is the same as the rate of
probability ¬‚ow of the (n ’ 1)-photon state to the n-photon state. The solution of this
recursion relation, subject to the normalization condition

pn (∞) = 1 , (18.38)

is the Bose“Einstein distribution
(ncav )
pn (∞) = . (18.39)
(ncav + 1)
Examples of the master equation

18.4.2 Two-level atom
For the two-level atom, the sample“reservoir interaction Hamiltonian HSE is given by
eqns (14.131)“(14.133). In this case, the operator F (t) in eqn (18.19) is the sum of two
terms: F (t) = Fsp (t) + Fpc (t), that are respectively given by
√ √
Fsp (t) = i w21 b† (t) S 12 (t) = i w21 b† (t) σ ’ (t) (18.40)
in in

√ √
w11 c† (t) S 11 (t) + w22 c† (t) S 11 (t)
Fpc (t) = i 1,in 2,in

√ 1 ’ σ z (t) 1 + σ z (t)
w11 c† (t) w22 c† (t)
=i + .
1,in 2,in
2 2
The envelope operators σ ’ (t) and σ z (t) are related to the environment-picture forms
by σ ’ (t) = σ’ (t) exp (iω21 t) and σ z (t) = σz (t), while the operators b† (t) and c† (t)
in q,in
are the in-¬elds de¬ned by eqns (14.146) and (14.147) respectively.
We will assume that the reservoirs are uncorrelated, i.e. ρE (0) = ρsp (0) ρpc (0),
and that the individual reservoirs are stationary. These assumptions guarantee that
most of the possible terms in the double commutator will vanish when the partial trace
over the environment is carried out.
After performing the invigorating algebra suggested in Exercise 18.4.2, the surviv-
ing terms yield the Schr¨dinger-picture master equation

ρS (t) = LS ρS (t) + Ldis ρS (t) , (18.42)
where the Hamiltonian part,
1 ω21
LS ρS (t) = σz + HS1 (t) , ρS (t) , (18.43)
i 2
includes HS1 (t). The dissipative part is
Ldis ρS (t) = ’ (nsp + 1) {[σ+ , σ’ ρ] + [ρσ+ , σ’ ]}
’ nsp {[σ’ , σ+ ρ] + [ρσ’ , σ+ ]}
’ [ρσz , σz ] , (18.44)
where nsp is the average number of reservoir excitations (photons) at the transition
frequency ω21 . The phase-changing rate in the last term is
wpc = (2npc,q + 1) wqq , (18.45)
2 q=1

where npc,q is the average number of excitations in the phase-perturbing reservoir
coupled to the atomic state |µq .
The master equation

18.5 Phase space methods
In Section 5.6.3 we have seen that the density operator for a single cavity mode can be
described in the Glauber“Sudarshan P (±) representation (5.165). As we will show be-
low, this representation provides a natural way to express the operator master equation
(18.31) as a di¬erential equation for the P (±)-function. In single-mode applications”
and also in some more complex situations”this equation is mathematically identical
to the Fokker“Planck equation studied in classical statistics (Risken, 1989, Sec. 4.7).
By de¬ning an atomic version of the P -function”as the Fourier transform of a
properly chosen quantum characteristic function”it is possible to apply the same
techniques to the master equation for atoms, but we will restrict ourselves to the
simpler case of a single mode of the radiation ¬eld. The application to atomic master
equations can be found, for example, in Haken (1984, Sec. IX.2) or Walls and Milburn
(1994, Chap. 13).
For the discussion of the master equation in terms of P (±), it is better to use
the alternate convention in which P (±) is regarded as a function, P (±, ±— ), of the
independent variables ± and ±— . In this notation the P -representation is
d2 ±
|± P (±, ±— ; t) ±| .
ρS (t) = (18.46)
The function P (±, ±— ; t) is real and satis¬es the normalization condition
d2 ±
P (±, ±— ; t) = 1 , (18.47)
but it cannot always be interpreted as a probability distribution. The trouble is that,
for nonclassical states, P (±, ±— ; t) must take on negative values in some region of the
18.5.1 The Fokker“Planck equation
In order to use the P -representation in the master equation, we must translate the
products of Fock space operators, e.g. a, a† , and ρ, in the master equation into the
action of di¬erential operators on the c-number function P (±, ±— ; t). For this purpose
it is useful to write the coherent state |± as
|± = e’|±| |±; B ,
where the Bargmann state |±; B is

√ |n .
|±; B = (18.49)

The virtue of the Bargmann states is that they are analytic functions of ±. More
precisely, for any ¬xed state |Ψ the c-number function


Ψ |±; B = Ψ |n (18.50)

is analytic in ±. In the same sense, ±; B| is an analytic function of ±— , so it is inde-
pendent of ±.
Phase space methods

Since |±; B is proportional to |± , the action of a on the Bargmann states is just

a |±; B = ± |±; B . (18.51)

The action of a† is found by using eqn (18.49) to get

±n √

a |±; B = n + 1 |n + 1

|±; B .
= (18.52)
The adjoint of this rule is

±; B| a =
±; B| . (18.53)
In the Bargmann notation, the P -representation of the density operator is

d2 ±
P (±, ±— ; t) e’|±| |±; B ±; B| .
ρS (t) = (18.54)

The rule (18.51) then gives

d2 ±
P (±, ±— ; t) e’|±| ± |±; B
aρS (t) = ±; B|
d2 ±
±P (±, ±— ; t) |± ±| .
= (18.55)

Applying the rule (18.52) yields

d2 ± ‚

P (±, ±— ; t) e’|±|
|±; B
a ρS (t) = ±; B| , (18.56)
π ‚±

but this is not expressed in terms of a di¬erential operator acting on P (±, ±— ; t).
Integrating by parts on ± leads to the desired form:

d2 ± ‚ —
a† ρS (t) = ’ P (±, ±— ; t) e’±± |±; B ±; B|
π ‚±
d± ‚
±— ’ P (±, ±— ; t) |± ±| .
= (18.57)
π ‚±

This result depends on the fact that the normalization condition requires P (±, ±— ; t)
to vanish as |±| ’ ∞.
Combining eqns (18.55) and (18.57) with their adjoints gives us the translation
aρS (t) ” ±P (±, ±— ; t) ρS (t) a ” (± ’ ‚/‚±— ) P (±, ±— ; t)
a† ρS (t) ” (±— ’ ‚/‚±) P (±, ±— ; t) ρS (t) a† ” ±— P (±, ±— ; t) .
The master equation

Applying the rules in eqn (18.58) to eqn (18.32)”for the simple case with HS1 =
0”and to eqn (18.33) yields the translations
1 ‚ ‚—
ω0 a† a, ρS (t) ” iω0 ± P (±, ±— ; t)
LS ρS (t) = ±’ (18.59)

i ‚± ‚±
“ ‚ ‚
[±P (±, ±— ; t)] + [±— P (±, ±— ; t)]
Ldis ρS (t) ” —
2 ‚± ‚±
P (±, ±— ; t) ,
+ “ncav (18.60)

for the Hamiltonian and dissipative Liouville operators respectively.
In the course of carrying out these calculations, it is easy to get confused about
the correct order of operations. The reason is that products like a† aρ”with operators
standing on the left of ρ”and products like ρa† a”with operators standing to the
right of ρ”are both translated into di¬erential operators acting from the left on the
function P (±, ±— ; t).
Studying a simple example, e.g. carrying out a direct derivation of both a† aρ and
ρa† a, shows that the order of the di¬erential operators is reversed from the order of the
Fock space operators when the Fock space operators stand to the right of ρ. Another
way of saying this is that one should work from the inside to the outside; the ¬rst
di¬erential operator acting on P corresponds to the Hilbert space operator closest to
ρ. This rule gives the correct result for Fock space operators to the left or to the right
of ρ.
The master equation for an otherwise unperturbed cavity mode is, therefore, rep-
resented by
‚ ‚ ‚
P (±, ±— ; t) = [Z (±) P (±, ±— ; t)] + [Z — (±) P (±, ±— ; t)]

‚t ‚± ‚±

P (±, ±— ; t) ,
+ “ncav (18.61)


+ iω0 ± . (18.62)
Z (±) =
We can achieve a ¬rmer grip on the meaning of this equation by changing variables
from (±, ±— ) to u = (u1 , u2 ), where u1 = Re ± and u2 = Im ±. In these variables,
P (±, ±— ; t) = P (u; t), and the ±-derivative is
‚ 1 ‚ ‚
= . (18.63)
‚± 2 ‚u1 ‚u2
In this notation, the master equation takes the form of a classical Fokker“Planck
equation in two dimensions:
‚ D0 2
P (u; t) = ’∇ · [F (u) P (u; t)] + ∇ P (u; t) , (18.64)
‚t 2
Phase space methods

D0 = “ncav /2 (18.65)
is the di¬usion constant, and we have introduced the following shorthand notation:

‚X1 ‚X2
∇·X = + ,
‚u1 ‚u2
“ “
F (u) = (’ Re Z, ’ Im Z) = ω0 u 2 ’ u1 , ’ω0 u1 ’ u2 , (18.66)
2 2
2 2
‚ ‚
+ .
‚u1 ‚u2

The ¬rst- and second-order di¬erential operators in eqn (18.64) are respectively called
the drift term and the di¬usion term.

A Classical Langevin equations
The Fokker“Planck equation (18.64) is a special case of a general family of equations
of the form
‚ 1 ‚ ‚
P (u; t) = ’∇ · [F (u, t) P (u; t)] + Dmn (u, t) P (u; t) , (18.67)
‚t 2 m=1 n=1 ‚um ‚un

where u = (u1 , . . . , uN ), F = (F1 , . . . , FN ), Dmn is the di¬usion matrix, and

X · Y = X1 Y1 + · · · + XN YN . (18.68)

For the two-component case, given by eqn (18.64), the di¬usion matrix is diagonal,
Dmn = D0 δmn , so it has a single eigenvalue D0 > 0. The corresponding condition in
the general N -component case is that all eigenvalues of the di¬usion matrix D are
positive, i.e. D is a positive-de¬nite matrix. In this case D has a square root matrix
B that satis¬es D = BB T .
When D is positive de¬nite, then eqn (18.67) is exactly equivalent to the set of
classical Langevin equations (Gardiner, 1985, Sec. 4.3.5)
dun (t)
= Cn (u, t) + Bnm (u, t) wm (t) , (18.69)
dt m=1

where the un s are stochastic variables and the wm s are independent white noise sources
of unit strength, i.e. wm (t) = 0 and

wm (t) wn (t ) = δmn δ (t ’ t ) . (18.70)

In particular, the Langevin equations corresponding to eqn (18.64) are

du (t)
= F (u) + D0 w (t) . (18.71)
¼ The master equation

These real Langevin equations are essential for numerical simulations, but for ana-
lytical work it is useful to write them in complex form. This is done by combining
± = u1 + iu2 with eqns (18.62) and (18.66) to get

d± (t)
= ’Z (± (t)) + 2D0 · (t) , (18.72)
where ± (t) is a complex stochastic variable, and

· (t) = √ [w1 (t) + iw2 (t)] (18.73)
is a complex white noise source satisfying

· — (t) · (t ) = δ (t ’ t ) .
· (t) = 0 , · (t) · (t ) = 0 , (18.74)

The equivalence of the Fokker“Planck equation and the classical Langevin equa-
tions for a positive-de¬nite di¬usion matrix is important in practice, since the nu-
merical simulation of the Langevin equations is usually much easier than the direct
numerical solution of the Fokker“Planck equation itself.
For some problems”e.g. when the appropriate reservoir is described by a squeezed
state”the di¬usion matrix derived from the Glauber“Sudarshan P -function is not
positive de¬nite, so the Fokker“Planck equation is not equivalent to a set of classical
Langevin equations. In such cases, another representation of the density operator may
be more useful (Walls and Milburn, 1994, Sec. 6.3.1).

18.5.2 Applications of the Fokker“Planck equation
A Coherent states are robust
Let us begin with a simple example in which ncav = 0, so that the di¬usion term in eqn
(18.64) vanishes. If we interpret the reservoir oscillators as phonons in the cavity walls,
then this model describes the idealized situation of material walls at absolute zero.
Alternatively, the reservoir could be de¬ned by other modes of the electromagnetic
¬eld, into which the particular mode of interest is scattered by a gas of nonresonant
atoms. In this case, it is natural to assume that the initial reservoir state is the vacuum.
In other words, the universe is big and dark and cold.
The terms remaining after setting ncav = 0 can be rearranged to produce

P (u; t) + F (u) · ∇P (u; t) = “P (u; t) . (18.75)
Let us study the evolution of a ¬eld state initially de¬ned by P (u; 0) = P0 (u). The
general technique for solving linear, ¬rst-order, partial di¬erential equations like eqn
(18.75) is the method of characteristics (Zauderer, 1983, Sec. 2.2), but we will employ
an equivalent method that is well suited to the problem at hand.
The ¬rst step is to introduce an integrating factor, by setting

P (u; t) = P (u; t) e“t , (18.76)
Phase space methods

so that

P (u; t) + F (u) · ∇P (u; t) = 0 . (18.77)
The second step is to transform to new variables (u , t ) by

u = V (u, t) , t = t , (18.78)

where we require u = u at t = 0, and also assume that the function V (u, t) is linear
in u, i.e.
Vn (u, t) = Gnm (t) um . (18.79)

The reason for trying a linear transformation is that the coe¬cient vector,

’“/2 ω0
, (18.80)
Fj (u) = Wjl ul , where W =
’ω0 ’“/2

is itself linear in u.
The chain rule calculation in Exercise 18.5 yields expressions for the operators ‚/‚t
and ‚/‚ul in terms of the new variables, so that eqn (18.77) becomes

‚ dG (t) ‚
P (u ; t ) + G (t) W + ul P (u ; t) = 0 . (18.81)
‚t dt ‚un
n l

Choosing the matrix G (t) to satisfy

dG (t)
+ G (t) W = 0 (18.82)
ensures that the coe¬cient of ‚/‚un vanishes identically in u, and this in turn simpli¬es
the equation for P (u ; t ) to

P (u ; t ) = 0 . (18.83)
Thus P (u ; t ) = P (u ; 0), but t = 0 is the same as t = 0, so P (u ; t ) = P0 (u ).
In this way the solution to the original problem is found to be

P (u, t) = e“t P0 (V (u, t)) , (18.84)

and the only remaining problem is to evaluate V (u, t). This is most easily done by
writing G (t) as
b (t) b2 (t)
G (t) = 1 , (18.85)
c1 (t) c2 (t)
and substituting this form into eqn (18.82). This yields simple di¬erential equations
for the vectors b (t) and c (t), with initial conditions b (0) = (1, 0) and c (0) = (0, 1).
The solution of these auxiliary equations gives
¾ The master equation

cos (ω0 t) ’ sin (ω0 t)
G (t) = e“t/2 R (t) = e“t/2 , (18.86)
sin (ω0 t) cos (ω0 t)
so that
P (u; t) = P0 e“t/2 R (t) u e“t . (18.87)
Thus, in the absence of the di¬usive term, the shape of the distribution is un-
changed; the argument u is simply scaled by exp (“t/2) and rotated by the angle ω0 t.
In the complex-± description the solution is given by

P (±, ±— ; t) = P0 e(“/2+iω0 )t ±, e(“/2’iω0 )t ±— e“t . (18.88)

This result is particularly interesting if the ¬eld is initially in a coherent state |±0 ,
i.e. the initial P -function is P0 (u) = δ2 (u ’ u0 ). In this case, the standard properties
of the delta function lead to
P (u; t) = e“t δ2 e“t/2 R (t) u ’ u0 = δ2 (u ’ u (t)) , (18.89)

u (t) = e’“t/2 [R (t)] u0 . (18.90)
The conclusion is that a coherent state interacting with a zero-temperature reservoir
will remain a coherent state, with a decaying amplitude
± (t) = ±0 e’“t/2 e’iω0 t . (18.91)
Consequently, the time-dependent joint variance of a† and a vanishes at all times:
V a† , a ; t = ± (t) a† a ± (t) ’ ± (t) a† ± (t) ± (t) |a| ± (t) = 0 . (18.92)
In other words, coherent states are robust: scattering and absorption will not destroy
the coherence properties, as long as the environment is at zero temperature.
This apparently satisfactory result raises several puzzling questions. The ¬rst is that
the initially pure state remains pure, even after interaction with the environment. This
seems to contradict the general conclusion, established in Section 18.1, that interaction
with a reservoir inevitably produces a mixed state for the sample.
The resolution of this discrepancy is that the general argument is true for the exact
theory, while the master equation is derived with the aid of the approximation”see
eqn (18.17)”that back-action of the sample on the reservoir can be neglected. This
means that the robustness property of the coherent states is only as strong as the
approximations leading to the master equation. Furthermore, we will see in Section
18.6 that coherent states are the only pure states that can take advantage of this
loophole in the general argument of Section 18.1.
The second di¬culty with the robustness of coherent states is that it seems to
violate the ¬‚uctuation dissipation theorem. The ¬eld su¬ers dissipation, but there is
no added noise. Consequently, it is a relief to realize that the strength of the noise
term in the equivalent classical Langevin equations (18.71) vanishes for ncav = 0.
Further reassurance comes from the operator Langevin approach, in particular eqn
(14.74), which shows that the strength of the Langevin noise operator also vanishes
for ncav = 0.
Phase space methods

B Thermalization of an initial coherent state
The coordinates de¬ned by eqn (18.78) are also useful for solving eqn (18.64), the
Fokker“Planck equation with di¬usion. According to eqn (18.86), the transformation
from u to u is a rotation followed by scaling with exp (“t/2). The operator ∇2 on
the right side of eqn (18.64) is invariant under rotations, so ∇2 = e“t ∇ 2 and the
Fokker“Planck equation becomes

‚ D0 “t 2
e ∇ P (u ; t ) ,
P (u ; t ) = (18.93)
‚t 2
which is the di¬usion equation with a time-dependent di¬usion coe¬cient. The Fourier
d2 u P (u ; t ) e’iq ·u ,
P (q ; t ) = (18.94)

then satis¬es the ordinary di¬erential equation

d D0 “t 2
P (q ; t ) = e q P (q ; t ) , (18.95)
dt 2
which has the solution
D0 “t
e ’ 1 q 2 P 0 (q ) .
P (q ; t ) = exp (18.96)

For the initial coherent state, P0 (u) = δ2 (u ’ u0 ), one ¬nds

P 0 (q ) = exp [’iq · u0 ] , (18.97)

and the inverse transform can be explicitly evaluated to yield

(u ’ u (t))2
exp ’ , (18.98)
P (u; t) =
πw (t) w (t)

where u (t) is given by eqn (18.90) and

w (t) = ncav 1 ’ e’“t . (18.99)

1/“) u (t) ’ 0, and the P -function approaches the thermal
For long times (t
distribution given by eqn (5.176); in other words, the ¬eld comes into equilibrium with
1/“, we see that w (t) ∼ ncav “t
the cavity walls as expected. At short times, t 1
and the initial delta function is recovered.

C A driven mode in a lossy cavity
In Section 5.2 we presented a simple model for generating a coherent state of a sin-
gle mode in a lossless cavity. We can be sure that losses will be present in any real
experiment, so we turn to the Fokker“Planck equation for a more realistic treatment.
The master equation

The o¬-resonant term in the Heisenberg equation (5.38) de¬ning our model can
safely be neglected, so the situation is adequately represented by the simpli¬ed Hamil-
HS (t) = ω0 a† a ’ W e’i„¦t a† ’ W — ei„¦t a , (18.100)

that leads to the Liouville operator

1 1
ω0 a† a, ρS (t) ’ W e’i„¦t a† + W — ei„¦t a, ρS (t) .
LS ρS (t) = (18.101)
i i

After including the new terms in the master equation and applying the rules (18.58),
one ¬nds an equation of the same form as eqn (18.61), except that the Z (±) function
is replaced by

+ iω0 ± ’ iW e’i„¦t .
Z (±) = (18.102)
Instead of directly solving the Fokker“Planck equation, it is instructive to use the
equivalent set of classical Langevin equations. Substituting the new Z (±) function
into the general result (18.72) yields

d± (t) “
+ iω0 ± (t) + iW e’i„¦t +
=’ 2D0 · (t) , (18.103)
dt 2

which has the solution

± (t) = ± (0) e’(iω0 +“/2)t + ±coh (t) + 2D0 ‘ (t) , (18.104)

e’(iω0 +“/2)t e(i∆+“/2)t ’ 1
±coh (t) = (18.105)
i∆ + “/2
is a de¬nite (i.e. nonrandom) function, and
dt1 e’(iω0 +“/2)(t’t1 ) · (t1 ) .
‘ (t) = (18.106)

The initial value, ± (0), is a complex random variable, not a de¬nite complex number.
The average of any function f (± (0)) is given by

d2 ± (0) P0 (± (0) , ±— (0)) f (± (0)) ,
f (± (0)) = (18.107)

but special problems arise if the initial state is not classical. For a classical state”
i.e. P0 (±, ±— ) 0”standard methods can be used to draw ± (0) randomly from the
distribution, but these methods fail when P0 (±, ±— ) is negative. For these nonclas-
sical states, the c-number Langevin equations are of doubtful utility for numerical
Phase space methods

For the problem at hand, the initial state is the vacuum, with the positive distri-
bution P0 (± (0) , ±— (0)) = δ (± (0)). The initial value ± (0) and the noise term ‘ (t)
both have vanishing averages, so the average value of ± (t) is given by

± (t) = ±coh (t) . (18.108)

For the nondissipative case, “ = D0 = 0, the average agrees with eqn (5.41); but,
when dissipation is present, the long time (t 1/“) solution approaches
e’i„¦t .
± (t) = (18.109)
i∆ + “/2

Thus the decay of the average ¬eld due to dissipation”shown by eqn (18.91)”is
balanced by radiation from the classical current, and the average ¬eld amplitude has
a de¬nite phase determined by the phase of the classical current. This would also be
true if the sample were described by the coherent state ρcoh (t) = |±coh (t) ±coh (t)|,
so it will be necessary to evaluate second-order moments in order to see if eqn (18.104)
corresponds to a true coherent state.
We will ¬rst investigate the coherence properties of the state by using the explicit
solution (18.104) to get
’(iω0 +“/2)t
± (t) = ± (0) e + ±coh (t) + 2D0 ‘ (t)

= ±2 (t) + 2D0 ‘2 (t) . (18.110)

The simple form of the second line depends on two facts: (i) ±coh (t) is a de¬nite
function; and (ii) the distribution of ± (0) is concentrated at ± (0) = 0. A further
simpli¬cation comes from using eqn (18.105) to evaluate ‘2 (t) . The result is a double
integral with integrand proportional to · (t1 ) · (t2 ) , but eqn (18.74) shows that this
average vanishes for all values of t1 and t2 . The ¬nal result is then

±2 (t) = ±2 (t) = ± (t) 2 , (18.111)

which also agrees with the prediction for a true coherent state.
Before proclaiming that we have generated a true coherent state in a lossy cavity,
we must check the remaining second-order moment, |± (t)| , which represents the
average of the number operator a† a. Since ± (t) is concentrated at the origin, we can
simplify the calculation by setting ± (0) = 0 at the outset. This gives us
2 2 2
|± (t)| = |±coh (t)| + 2D0 |‘ (t)| . (18.112)

Combining eqns (18.106) and (18.74) leads to
t t
dt2 e’(’iω0 +“/2)(t’t1 ) e’(iω0 +“/2)(t’t2 ) · — (t1 ) · (t2 )
|‘ (t)| = dt1
0 0
1 ’ e’“t
’“(t’t1 )
dt1 e = , (18.113)

The master equation

so that
1 ’ e’“t = ncav 1 ’ e’“t ,
|± (t)|2 ’ |±coh (t)|2 = (18.114)

where we used eqn (18.65) to get the ¬nal result. The left side of this equation would
vanish for a true coherent state, so the state generated in a lossy cavity is only coherent
if ncav = 0, i.e. if the cavity walls are at zero temperature.

The Lindblad form of the master equation—
The master equations (18.31) and (18.42) share three important properties.
(a) The trace condition, Tr [ρS (t)] = 1, is conserved.
(b) The positivity of ρS is conserved, i.e. Ψ |ρS (t)| Ψ 0 for all states |Ψ and all
times t.
(c) The equations are derivable from a model of the sample interacting with a collec-
tion of reservoirs.
The most general linear, dissipative time evolution that satis¬es (a), (b), and (c)
is given by
= LS ρS + Ldis ρS , (18.115)
LS ρS = [HS (t) , ρS ] (18.116)
describes the Hamiltonian evolution of the sample, and the dissipative term has the
Lindblad form (Lindblad, 1976)
1 † † †
Ldis ρS = ’ Ck Ck ρS + ρS Ck Ck ’ 2Ck ρS Ck . (18.117)

Each of operators C1 , C2 , . . . , CK acts on the sample space HS and there can be a
¬nite or in¬nite number of them, depending on the sample under study.
One can see by inspection that there are two Lindblad operators, i.e. K = 2, for
the single-mode master equation (18.31):
“ncav a† .
C1 = “ (ncav + 1)a , C2 = (18.118)
A slightly longer calculation”see Exercise 18.6”shows that there are three operators
for the master equation (18.42) describing the two-level atom.
The Lindblad form (18.117) for the dissipative operator can be used to investigate
a variety of questions. For example, in Section 2.3.4 we introduced a quantitative
measure of the degree of mixing by de¬ning the purity of the state ρS as P (t) =
Tr ρ2 (t) 1. One can show from eqn (18.115) that the time derivative of the purity
d † †
P (t) = ’2 Tr ρS (t) Ck Ck ρS (t) ’ ρS (t) Ck ρS (t) Ck . (18.119)
At ¬rst glance, it may seem natural to assume that interaction with the environment
can only cause further mixing of the sample state, so one might expect that the time
Quantum jumps

derivative of the purity is always negative. If ρS (0) is a mixed state this need not be
true. For example, the purity of a thermal state would be increased by interaction
with a colder reservoir, as seen in Exercise 18.7.
On the other hand, for a pure state there is no way to go but down; therefore, the
intuitive expectation of declining purity should be satis¬ed. In order to check this, we
evaluate eqn (18.119) for an initially pure state ρS (0) = |Ψ Ψ|, to ¬nd
d † †
P (t) = ’2 Ψ Ck Ck Ψ ’ Ψ Ck Ψ Ψ |Ck | Ψ
dt t=0 k=1

= ’2 Ψ δCk δCk Ψ 0, (18.120)

δCk = Ck ’ Ψ |Ck | Ψ . (18.121)
Thus the Lindblad form guarantees the physically essential result that initially pure
states cannot increase in purity (Gallis, 1996).
The appearance of an inequality like eqn (18.120) prompts the following question:
Are there any physical samples possessing states that saturate the inequality? We can
answer this question in one instance by studying the master equation (18.31) with a

zero-temperature reservoir. In this case eqn (18.118) gives us C2 = 0 and C1 = “a,
so that eqn (18.120) becomes
d †
P (t) = ’2“ Ψ (a ’ ±) (a ’ ±) Ψ 0, (18.122)
dt t=0

where ± = Ψ |a| Ψ . The inequality can only be saturated if |Ψ satis¬es
a |Ψ = ± |Ψ , (18.123)
i.e. when |Ψ is a coherent state. For all other pure states, interaction with a zero-
temperature reservoir will decrease the purity, i.e. the state becomes mixed.

18.7 Quantum jumps
18.7.1 An elementary description of quantum jumps
The notion of quantum jumps was a fundamental part of the earliest versions of
the quantum theory, but for most of the twentieth century it was assumed that the
phenomenon itself would always be unobservable, since there were no experimental
methods available for isolating and observing individual atoms, ions, or photons.
This situation began to change in the 1980s with Dehmelt™s proposal (Dehmelt,
1982) for an improvement in frequency standards based on observations of a single
ion, and the subsequent development of electromagnetic traps (Paul, 1990) that made
such observations possible.
The following years have seen a considerable improvement in both experimental
and theoretical techniques. The improved experimental methods have made possible
the direct observation of the quantum jumps postulated by the founders of quantum
The master equation

Fig. 18.1 A three-level ion with dipole-al-
lowed transitions 3 ” 2 and 3 ” 1, indicated
by wavy arrows, and a dipole-forbidden tran-
sition 2 ’ 1, indicated by the light dashed
arrow. The heavy double arrows denote strong
incoherent couplings on the 3 ” 1 and 3 ’ 2

A A three-level model
It is always good to have a simple, concrete example in mind, so we will study a single,
trapped, three-level ion, with the level structure shown in Fig. 18.1. The dipole-allowed
transitions, 3 ’ 1 and 3 ’ 2, have Einstein-A coe¬cients “31 and “32 respectively, so
the total decay rate of level 3 is “3 = “31 + “32 . Since the dipole-forbidden transition,
2 ’ 1, has a unique ¬nal state, it is described by a single decay rate “2 , which is small
compared to both “31 and “32 .
In addition to the spontaneous emission processes, we assume that an incoher-
ent radiation source, at the frequency ω31 , drives the ion between levels 1 and 3 by
absorption and spontaneous emission. As explained in Section 1.2.2, both of these
processes occur with the rate W31 = B31 ρ (ω31 ), where ρ (ω31 ) is the energy density of
the external ¬eld and B31 is the Einstein-B coe¬cient.
When level 3 is occupied, the ion can isotropically emit ¬‚uorescent radiation, i.e.
radiation at frequency ω31 . Another way of saying this is that the ion scatters the pump
light in all directions. Consequently, observing the ¬‚uorescent intensity”say at right
angles to the direction of the pump radiation”e¬ectively measures the population of
level 3.
We will further assume that the levels 2 and 3 are closely spaced in energy, com-
pared to their separation from level 1, so that ω32 ω31 . From eqn (4.162) we know
that the Einstein-A coe¬cient is proportional to the cube of the transition frequency;
therefore, the transition rate for 3 ’ 2 will be small compared to the transition rate
for 3 ’ 1, i.e. “32 “31 . In some cases, the small size of “32 may cause an excessive
delay in the transition from 3 to 2, so we also allow for an incoherent driving ¬eld on
the 3 ” 2 transition such that

“32 W32 = B32 ρ (ω32 ) “31 . (18.124)

Under these conditions, the ion will spend most of its time shuttling between levels
1 and 3, with infrequent transitions from 3 to the intermediate level 2. The forbidden
transition 2 ’ 1 occurs very slowly compared to 3 ’ 1 and 3 ’ 2, so level 2 e¬ectively
traps the occupation probability for a relatively long time. When this happens the
¬‚uorescent signal will turn o¬, and it will not turn on again until the ion decays back
to level 1. We will refer to these transitions as quantum jumps.1

1 Itwould be equally correct”but not nearly as exciting”to refer to this phenomenon as ˜inter-
rupted ¬‚uorescence™.
Quantum jumps

During the dark interval the ion is said to be shelved and |µ2 is called a shelving
state. The shelving e¬ect is emphasized when the 1 ” 3 transition is strongly saturated
and the state |µ2 is long-lived compared to |µ3 , i.e. when
W31 “3 “2 . (18.125)
During the bright periods when ¬‚uorescence is observed, the state vector |Ψion will
be a linear combination of |µ1 and |µ3 ; in other words, |Ψion is in the subspace H13 .

B A possible experimental realization
As a possible experimental realization of the three-level model, consider the intermit-
tent resonance ¬‚uorescence of the strong Lyman-alpha line, emitted by a singly-ionized
helium ion (He+ ) in a Paul trap. One advantage of this choice is that the spectrum is
hydrogenic, so that it can be calculated exactly.
The complementary relation between theory and experiment guarantees the pres-
ence of several real-world features that complicate the situation. The level diagram
in part (a) of Fig. 18.2 shows not one, but two intermediate states, 2S1/2 and 2P1/2 ,
that are separated in energy by the celebrated Lamb shift, ∆EL / = 14.043 GHz.
The 2S1/2 -level is a candidate for a shelving state, since there is no dipole-allowed
transition to the 1S1/2 ground state, but the 2P1/2 -level does have a dipole-allowed
transition, 2P1/2 ’ 1S1/2 . This adds unwanted complexity.
An additional theoretical di¬culty is caused by the fact that the dominant mech-
anism for the transition 2S1/2 ’ 1S1/2 is a two-photon decay. This is a problem,
because the reservoir model introduced in Section 14.1.1 is built on the emission or
absorption of single reservoir quanta; consequently, the standard reservoir model would
not apply directly to this case.
Fortunately, these complications can be exploited to achieve a closer match to our
simple model. The ¬rst step is to apply a weak DC electric ¬eld E 0 to the ion. In this



(a) (b)

Fig. 18.2 (a) Level diagram for the He+ ion. The spacing between the 2S1/2 and 2P1/2 levels
¬ ¬
« «
is exaggerated for clarity. (b) The unperturbed ¬2S1/2 and ¬2P1/2 states are replaced by the
Stark-mixed states |µ2 and |µ2 . The wavy arrows indicate dipole-allowed decays, while the
solid arrows indicate incoherent driving ¬elds. A dipole-allowed decay 2 ’ 1 is not shown,
since 2 is e¬ectively isolated by the method explained in the text.
¼ The master equation

application ˜weak™ means that the energy-level shift caused by the static ¬eld is small
compared to the Lamb shift, i.e.
2S1/2 |E 0 · d| 2P1/2 ∆EL . (18.126)
In this case there will be no ¬rst-order Stark shift, and the second-order Stark e¬ect
(Bethe and Salpeter, 1977) mixes the 2S1/2 and 2P1/2 states to produce two new
|µ2 = CS 2S1/2 + CP 2P1/2 , (18.127)
|µ2 = CS 2S1/2 + CP 2P1/2 , (18.128)
as illustrated in part (b) of Fig. 18.2.
A second-order perturbation calculation”using the Stark interaction HStark =
’d · E 0 ”shows that |CS | |CP |, i.e. |µ2 is dominantly like 2S1/2 , while |µ2 is
mainly like 2P1/2 . The states |µ1 and |µ3 of the simple model pictured in Fig. 18.1
are identi¬ed with 1S1/2 and 2P3/2 respectively.
Since neither |µ2 nor |µ2 has de¬nite parity, the dipole selection rules now allow
single-photon transitions 2 ’ 1 and 2 ’ 1. The rate for 2 ’ 1 is proportional
to |CP | , so a proper choice of |E 0 | will guarantee that the single-photon process
dominates the two-photon process, while still being slow compared to the rate “31 for
the Lyman-alpha transition.
By the same token, there are dipole-allowed transitions from 3 to both 2 and 2 .
The unwanted level 2 can be e¬ectively eliminated by applying a microwave ¬eld
resonant with the 3 ’ 2 transition, but not with the 3 ’ 2 transition. The strength
of this ¬eld can be adjusted so that the stimulated emission rate for 3 ’ 2 is large
compared to the spontaneous rates for 3 ’ 2 and 3 ’ 2, but small compared to
the stimulated and spontaneous rates for the 3 ’ 1 transition. These settings ensure
that the population of |µ2 will remain small at all times and that |µ2 is an e¬ective
shelving state.
The main practical di¬culty for this experiment is that the pump would have to
operate at the vacuum-UV wavelength, 30.38 nm, of the Lyman-alpha line of the He+
ion. One possible way around this di¬culty is to use the radiation from a synchrotron
light source.
The transition 3 ’ 2 is primarily due to the microwave-frequency transition
2P3/2 ’ 2S1/2 , which occurs at 44 GHz. Our assumption that the spontaneous emis-
sion rate for this transition is small compared to the transition rate for the Lyman-
alpha transition is justi¬ed by the rough estimate
A (Lyman alpha) ν (Lyman alpha)
∼ ∼ 1016 (18.129)
A 2P3/2 ’ 2S1/2 ν 2P3/2 ’ 2S1/2

(Bethe and Salpeter, 1977), which uses the values 2.47 — 1015 Hz and 11 GHz for the
Lyman-alpha transition and the 2P3/2 ’ 2S1/2 microwave transition in hydrogen
The combination of the low rate for the 3 ’ 2 transition and the long lifetime of
the shelving level 2 will permit easy observation of interrupted resonance ¬‚uorescence
at the helium ion Lyman-alpha line, i.e. quantum jumps.
Quantum jumps

For hydrogen, the lifetime of the 2P3/2 state is 1.595 ns, so the estimate (18.129)
tells us that the lifetime for the microwave transition is approximately 2 — 107 s, i.e.
of the order of a year. The lifetime of the same transition for a hydrogenic ion scales
as Z ’4 , so for Z = 2 the microwave transition lifetime is 1.5 — 106 seconds, which is
about a month.
This is still a rather long time to wait for a quantum jump. The solution is to
adjust the strength of the resonant microwave ¬eld driving the 3 ” 2 transition to
bring this lifetime within the limits of the experimentalist™s patience.

C Rate equation analysis
The assumption that the driving ¬eld is incoherent allows us to extend the rate equa-
tion approximation (11.190) for two-level atoms to our simple model to get
= ’ (“3 + W31 + W32 ) P3 + W31 P1 + W32 P2 , (18.130)
= ’ (“2 + W32 ) P2 + (“32 + W32 ) P3 , (18.131)
= ’W31 P1 + (W31 + “31 ) P3 + “2 P2 . (18.132)
Adding the equations shows that the sum of the three probabilities is constant:
P1 + P2 + P3 = 1 . (18.133)
The inequalities (18.125) suggest that the adiabatic elimination rule (11.187) can
be applied to the rate equations (18.130)“(18.132). To see how the rule works in this
case, it is useful to express the rate equations in terms of the probability P31 = P3 + P1
that the ionic state is in H13 , and the inversion D31 = P3 ’ P1 . The new form of the
rate equations is
d 1 1
D31 = ’ 2W31 + “31 + “32 + W32 D31
dt 2 2
1 1
’ “31 + “32 + W32 P31 + (W32 ’ “2 ) P2 , (18.134)
2 2
d 1 1
P31 = ’ (“32 + W32 ) P31 ’ (“32 + W32 ) D31 + (“2 + W32 ) P2 , (18.135)
dt 2 2
d 1 1
P2 = ’ (“2 + W32 ) P2 + (“2 + W32 ) P31 + (“2 + W32 ) D31 . (18.136)
dt 2 2
The rate multiplying D31 on the right side of eqn (18.134) is much larger than any
other rate in the equations; therefore, D31 (t) will rapidly decay to the steady-state
solution of eqn (18.134), i.e.

W32 ’ “2
“31 + 1 “32 + 1 W32
=’ 2 2
D31 P31 + P2 . (18.137)
2W31 + “31 + 1 “32 + 1 W32 2W31 + “31 + 1 “32 + 1 W32
2 2 2 2

The coe¬cients of the probabilities P31 and P2 are very small, so we can set D31 0
in the rest of the calculation.
¾ The master equation

With this approximation, the remaining rate equations are
= ’Ron P2 + Ro¬ P31 (18.138)
= Ron P2 ’ Ro¬ P31 , (18.139)
where Ron = “2 + W32 is the rate at which the ¬‚uorescence turns on, and Ro¬ =
(“32 + W32 ) /2 is the rate at which ¬‚uorescence turns o¬. Solving eqns (18.138) and
(18.139) for P31 (t) yields
P31 (t) = P31 (0) e’(Ron +Roff )t + 1 ’ e’(Ron +Roff )t . (18.140)
Ron + Ro¬
The ¬‚uorescent intensity IF (t) is proportional to P31 (t), so IF (t) evolves smoothly
from its initial value IF (0) to the steady-state value
IF ∝ . (18.141)
Ron + Ro¬
This result is completely at odds with the ¬‚ickering on-and-o¬ behavior predicted
above. The source of this discrepancy is the fact that the quantities P1 , P2 , and P3 in
the rate equations (18.138) and (18.139) are unconditional probabilities. This means
that P1 , for example, is the probability that the ion is in level 1 without regard to its
past history or any other conditions. Another way of saying this is that P1 refers to
an ensemble of ions which have reached level 1 in all possible ways.
Before the development of single-ion traps, resonance ¬‚uorescence experiments
dealt with dilute atomic gases, and the total ¬‚uorescence signal would be correctly
described by eqn (18.140). In this case, the on-and-o¬ behavior of the individual atoms
would be washed out by averaging over the random ¬‚uorescence of the atoms in the
gas. For a single trapped ion, the smooth behavior in eqn (18.140) can only be recov-
ered by averaging over many observations, all starting with the ion in the same state,
e.g. the ground state.
In addition to the inability of the rate equations to predict quantum jumps, it
is also the case that statistical properties”such as the distribution of waiting times
between jumps”are beyond their reach. Thus any improvement must involve putting
in some additional information; that is, reducing the size of the ensemble.
The ¬rst step in this direction was taken by Cook and Kimble (1985) who intro-
duced the conditional probability P31,n (t, t + T ) that the ion is in H13 after making
n transitions between H13 and |µ2 during the interval (t, t + T ). The number of tran-
sitions de¬nes a subensemble of ions with this history. The complementary object
P2,n (t, t + T ) is the probability that the ion is in level |µ2 after n transitions between
H13 and |µ2 during the interval (t, t + T ).
By using the approximations leading to eqns (18.138) and (18.139), it is possi-
ble to derive an in¬nite set of coupled rate-like equations for P31,n (t, t + T ) and
P2,n (t, t + T ), with n = 0, 1, . . .. This approach permits the calculation of various
statistical features of the quantum jumps, but it is not easy to connect it with the
more re¬ned quantum-jump theories to be developed later on.
Quantum jumps

D A stochastic model
We will now consider a simple on-and-o¬ model which is qualitatively similar to the
more sophisticated quantum-jump theories. In this approach, the analytical treatment
based on conditional probabilities is replaced by an equivalent stochastic simulation.
We ¬rst assume that the ¬‚uorescent intensity can only have the values I = 0 (o¬)
or I = IF (on). If the signal is on at time t, then the probability that it will turn o¬
in the interval (t, t + ∆t) is ∆po¬ = Ro¬ ∆t. Conversely, if the signal is o¬ at time t,
then the probability that it will turn on in the interval (t, t + ∆t) is ∆pon = Ron ∆t.
For su¬ciently small ∆t, we can assume that only one of these events occurs.
The ¬‚uorescent intensities In at the discrete times tn = (n ’ 1) ∆t can then be
calculated by the following algorithm.

For In = 0 choose a random number r in (0, 1) ;
then set In+1 = 0 if ∆pon < r or In+1 = IF if ∆pon > r .
For In = IF choose a random number r in (0, 1) ;
then set In+1 = IF if ∆po¬ < r or In+1 = 0 if ∆po¬ > r .

The random choices in this algorithm are a special case of the rejection method (Press
et al., 1992, Sec. 7.3) for choosing random variables from a known distribution.
From a physical point of view, the algorithm is an approximate embodiment of the
collapse postulate for measurements in quantum theory. The value In is the outcome of
a measurement of the ¬‚uorescent intensity at t = tn , so it corresponds to a collapse of
the state vector of the ion into the state with the value In . If In+1 = In the subsequent
collapse at t = tn+1 is into the same state as at t = tn . For In+1 = In the collapse at
tn+1 is into the other state, so we see a quantum jump.
A typical2 sequence of quantum jumps is shown in Fig. 18.3. Random sequences of
binary choices (dots and dashes) of this kind are called random telegraph signals.
This plot exhibits the expected on-and-o¬ behavior for a single ion, but the smooth
¬‚uorescence curve predicted by the rate equations is nowhere to be seen.
In order to recover an approximation to eqn (18.140), we consider M experiments,
all starting with I1 = IF , and de¬ne the average ¬‚uorescent intensity at time tn by
In,av = In,j , (18.143)
M j=1

where In,j is the ¬‚uorescent intensity at time tn for the jth run. A comparison of In,av
with the values predicted by eqn (18.140) is shown in Fig. 18.4, for M = 100.

E Experimental evidence
We have demonstrated a simple model displaying quantum jumps and a plausible
experimental realization for it, but the question remains if any such phenomena have
2 Thestochastic algorithm (18.142) gives a di¬erent plot for each run with the same input para-
meters. The ˜typical™ plot shown here was chosen to illustrate the e¬ect most convincingly. This kind
of data selection is not unknown in experimental practice.
The master equation







5 10 15 20

Fig. 18.3 Normalized ¬‚uorescent intensity I/IF versus time (in units of the radiative lifetime
1/“b of the shelving state). In these units, Ron = 1.6, Ro¬ = 0.3, and ∆t = 0.1. The initial
intensity is I (0) = IF .






2.5 5 7.5 10 12.5 15 17.5 20

Fig. 18.4 Fluorescent intensity (normalized to IF and averaged over 100 runs) versus time
(measured in units of the radiative lifetime of the shelving state). The initial intensity in each
run is I (0) = IF , and the parameter values are those used in Fig. 18.3.

been seen in reality. For this evidence we turn to an experiment in which intermittent
¬‚uorescence was observed from a single, laser-cooled Ba+ ion in a radio frequency trap
(Nagourney et al., 1986).
The complementary relation between theory and experiment is in full play in this
case, as seen by comparing the level diagram for this experiment”shown in Fig. 18.5”
with Fig. 18.1. Fortunately, the complications involved in the real experiment do not
change the essential nature of the e¬ect, which is seen in Fig. 18.6.
Quantum jumps

62 23/2
62 21/2
614.2 nm

493.4 nm
52 ,5/2
649.9 nm
455.4 nm
52 ,3/2

62 51/2

Fig. 18.5 Level structure of Ba+ . The states in the simple three-level model discussed in
¬ ¬ ¬
« « «
the text are |µ1 = ¬62 S1/2 , |µ2 = ¬62 P3/2 , and |µ3 = ¬52 D5/2 , which is the shelf state.
The remaining states are only involved in the laser cooling process indicated by the heavy
solid lines. The 1 ” 2 transition is driven by an incoherent source (a lamp) indicated by the
light solid line. (Reproduced from Nagourney et al. (1986).)
Fluorescence photon counts / second

Lamp on
62 51/2


52 ,5/2
0 100 200
Time (s)

Fig. 18.6 A typical trace of the 493 nm ¬‚uorescence from the 62 P1/2 -level showing the
quantum jumps after the hollow cathode lamp is turned on. The atom is de¬nitely known to
be in the shelf level during the low ¬‚uorescence periods. (Reproduced from Nagourney et al.

Quantum jumps and the master equation—

Many features of quantum-jump experiments are well described”at least semi-quantit-
atively”by the rate equation approximation for the conditional probabilities, e.g.
P31,n , or by the equivalent stochastic simulation; but the rate equation model has
de¬nite limitations. The most important of these is the restriction to incoherent ex-
citation of the atomic states. Many experiments employ laser excitation, which is
inherently coherent in character.
The master equation

The e¬ort required to incorporate coherence e¬ects eventually led to the creation of
several closely related approaches to the problem of quantum jumps. These techniques
are known by names like the Monte Carlo wave function method, quantum trajectories,
and quantum state di¬usion. Sorting out the relations between them is a complicated
story, which we will not attempt to tell in detail. For an authoritative account, we
recommend the excellent review article of Plenio and Knight (1998) which carries the
history up to 1999.
We will present a brief account of the Monte Carlo wave function technique for the
solution of the master equation. The other approaches mentioned above are technically
similar; but they di¬er in the original motivations leading to them, in their physical
interpretations, and in the kinds of experimental situations they can address.
There are two complementary views of these theoretical approaches. One may re-
gard them simply as algorithms for the solution of the master equation, or as concep-
tually distinct views of quantum theory. The discussion therefore involves both com-
putational and fundamental physics issues. We will ¬rst consider the computational
aspects of the Monte Carlo wave function technique, and then turn to the conceptual
relations between this method and the approaches based on quantum trajectories or
quantum state di¬usion.
The master equation (18.115) is a di¬erential equation describing the time evolu-
tion of the sample density operator. Except in highly idealized situations”for which
analytical solutions are known”the solution of the master equation requires numeri-
cal methods. Even for the apparently simple case of a single cavity mode, the sample
Hilbert space HS is in¬nite dimensional, so the annihilation operator a is represented
by an in¬nite matrix. A direct numerical attack would therefore require replacing
HS by a ¬nite-dimensional space, e.g. the subspace spanned by the number states
|0 , . . . , |M ’ 1 . This would entail representing the creation and annihilation opera-
tors and the density operator by M — M matrices.
In some situations, such as those discussed in Section 18.5.2, an alternative ap-
proach is to replace the in¬nite-dimensional space HS by the two-dimensional quan-
tum phase space, and to use”for a restricted class of problems”the Fokker“Planck
equation (18.61) or the equivalent classical Langevin equation (18.72). In general, this
method will fail if the di¬usion matrix D is not positive de¬nite.
The master equation for an atom can also be represented by a Fokker“Planck
equation on a ¬nite-dimensional phase space, but the collection of problems amenable
to this treatment is restricted by the same kind of considerations, e.g. a positive-de¬nite
di¬usion kernel, that apply to the radiation ¬eld. In many cases the center-of-mass
motion of the atom can be neglected”or at least treated classically”so the sample
Hilbert space is ¬nite dimensional. In this situation the master equation for a two-level
atom is simply a di¬erential equation for a 2 — 2 hermitian matrix. This is equivalent
to a set of four coupled ordinary di¬erential equations, so it is not computationally
Unfortunately, in the real world of experimental physics, atoms often have more


. 21
( 27)