Ψe (t) Ck Ck Ψe (t)

With this notation, the contribution of |Ψe (t) to the second term in eqn (18.149) is

(“e (t) ∆t) ρe

meas (t), where

K

Pk |φek (t) φek (t)| ,

ρe e

(t) = (18.151)

meas

k=1

†

Ψe (t) Ck Ck Ψe (t)

Pk =

e

, (18.152)

“e (t)

and

K

†

“e (t) = Ψe (t) Ck Ck Ψe (t) (18.153)

k=1

is the total transition (quantum-jump) rate of |Ψe (t) into the collection of normalized

states de¬ned by eqn (18.150). Since the coe¬cients Pk satisfy 0 Pk 1 and

e e

K

Pk = 1 ,

e

(18.154)

k=1

they can be treated as probabilities.

With this interpretation, ρe

meas has the form (2.127) of the mixed state describing

the sample after a measurement has been performed, but before the particular outcome

is known. This suggests that we interpret the second term on the right side of eqn

(18.149) as a wave packet reduction resulting from a measurement-like interaction

with the reservoir.

After summing over the ensemble, eqn (18.148) becomes

†

ρS (t + ∆t) = Udis (∆t) ρS (t) Udis (∆t) + “ (t) ∆t ρmeas (t) , (18.155)

where

P e ρe

ρmeas (t) = meas (t) , (18.156)

e

Pe “e (t)

Pe = , (18.157)

“ (t)

and

Quantum jumps

Pe “e (t)

“ (t) = (18.158)

e

is the ensemble-averaged transition rate.

A The Monte Carlo wave function algorithm

In quantum theory, a system evolves smoothly by the Schr¨dinger equation until a

o

measurement event forces a discontinuous change. This feature is the basis for the

procedure described here.

It is plausible to expect that only one of the two terms in eqn (18.155)”dissipative

evolution or wave packet reduction”will operate during a su¬ciently small time step.

We will ¬rst describe the Monte Carlo wave function algorithm (MCWFA) that follows

from this assumption, and then show that the density operator calculated in this way

is an approximate solution of the master equation (18.115).

In order to simplify the presentation we assume that the initial ensemble is de¬ned

by

states {|˜1 , . . . , |˜M } ,

(18.159)

probabilities {P1 , . . . , PM } ,

so that the index e = 1, 2, . . . , M .

In each time step, a choice between dissipative evolution and wave packet reduc-

tion”i.e. a quantum jump”has to be made. For this purpose, we note that the prob-

ability of a quantum jump during the interval (t, t + ∆t) is ∆Pe (t) = “e (t) ∆t, where

“e (t) is the total transition rate de¬ned by eqn (18.153). The discrete scheme will

only be accurate if the jump probability during a time step is small, i.e. ∆Pe (t) 1.

Consequently, the time step ∆t must satisfy “e (t) ∆t 1.

With this preparation, we are now ready to state the algorithm for integrating the

master equation in the interval (0, T ).

(1) Set e = 1 and de¬ne the discrete times tn = (n ’ 1) ∆t, where 1 n N and

(N ’ 1) ∆t = T .

(2) At the initial time t = 0, set |Ψ (0) = |Ψe (0) = |˜e .

(3) For n = 2, . . . , N choose a random number r in the interval (0, 1). If ∆Pe (tn’1 ) < r

go to (a), and if ∆Pe (tn’1 ) > r go to (b). Since we have imposed ∆Pe (t) 1,

this procedure guarantees that quantum jumps are relatively rare interruptions of

continuous evolution.

(a) In this case there is no quantum jump, and the state vector is advanced from

tn’1 to tn by dissipative evolution followed by normalization:

Udis (∆t) |Ψe (tn’1 )

|Ψe (tn ) =

†

Ψe (tn’1 ) Udis (∆t) Udis (∆t) Ψe (tn’1 )

1’ Hdis |Ψe (tn’1 )

i∆t

= , (18.160)

1 ’ ∆Pe (tn’1 )

where the last line follows from the de¬nition (18.147) of Udis (∆t).

¼ The master equation

(b) In this case there is a quantum jump, and the new state vector is de¬ned

by choosing k randomly from {1, 2, . . . , K}”conditioned by the probability

distribution Pk de¬ned in eqn (18.152)”and setting

e

|Ψe (tn ) = |φek (tn’1 ) , (18.161)

i.e. |Ψe (tn ) jumps to one of the states permitted by the second term in eqn

(18.148).

(4) Repeat step (3) Ntraj times to get Ntraj discrete representations

{|Ψej (tn ) , 1 N } , j = 1, . . . , Ntraj

n (18.162)

of the state vector. These representations are distinct, due to the random choices

made in each time step. The density operator that evolves from the original pure

state |˜e is then given by

Ntraj

1

|Ψej (tn ) Ψej (tn )| .

ρe (tn ) = (18.163)

Ntraj j=1

(5) Replace e by e + 1. If e + 1 M go to step (2). If e + 1 > M go to step (6).

(6) The density operator ρ (t) that evolves from the initial density operator ρ (0)”

de¬ned by the ensemble (18.159)”is given by

M

Pe ρe (tn ) .

ρ (tn ) = (18.164)

e=1

The computational cost of this method scales as Ntraj N , where N is the dimen-

sionality of the sample Hilbert space HS . Consequently, the MCWFA would not be

very useful as a technique for solving the master equation, if the required number of

trials is itself of order N . Fortunately, there are applications with large N for which

one can get good statistics with Ntraj N.

B Proof that the MCWFA generates a solution

If each of the density operators ρe (t) satis¬es the master equation, then so will the

overall density operator de¬ned by eqn (18.164); therefore, it is su¬cient to give the

proof for a single ρe (t). For a su¬ciently large number of trials, the evolution of the

pure state operators,

ρej (tn ) = |Ψej (tn ) Ψej (tn )| , (18.165)

is e¬ectively given by step (2a) with probability 1 ’ ∆Pe (tn’1 ) and by step (2b) with

probability ∆Pe (tn’1 ). In other words,

ρej (tn ) = (1 ’ ∆Pe (tn’1 )) Ψdis (tn ) Ψdis (tn )

ej ej

Pk (tn’1 ) |φek (tn’1 ) φek (tn’1 )| ,

e

+ ∆Pe (tn’1 ) (18.166)

k

½

Quantum jumps

where

1’ Hdis |Ψej (tn’1 )

i∆t

Ψdis (tn ) = . (18.167)

ej

1 ’ ∆Pe (tn’1 )

The |φek (tn’1 ) s are de¬ned by substituting |Ψej (tn’1 ) for |Ψe (tn’1 ) in eqn (18.150).

Using the de¬nitions of ∆Pe , Pk , and Hdis in this equation and neglecting O ∆t2 -

e

terms leads to

ρej (tn ) ’ ρej (tn’1 ) i

= ’ [HS , ρej (tn’1 )] + Ldis ρej (tn’1 ) . (18.168)

∆t

Averaging this result over the trials, according to eqn (18.163), and taking the limit

∆t ’ 0 shows that ρe (t) satis¬es the master equation (18.115).

Laser-induced ¬‚uorescence—

18.7.4

For a concrete application of the MCWFA, we return to the trapped three-level ion

considered in Section 18.7.1. For this example, however, we replace the incoherent

source driving 3 ” 1 by a coherent laser ¬eld E L e’iωL t that is close to resonance,

i.e. |ωL ’ ω31 | ωL . In the interests of simplicity, we also drop the ¬eld driving

3 ” 2. The semiclassical approximation for the laser is applied by substituting E(+) ’

E L e’iωL t in the general results (11.36) and (11.40) of Section 11.1.4.

In the resonant wave approximation, the Schr¨dinger-picture Hamiltonian is HS =

o

HS0 + HS1 , where

HS0 = q Sqq , (18.169)

q

HS1 = „¦L S31 e’iωL t + HC , (18.170)

and „¦L = ’d31 · E L / is the Rabi frequency for the laser driving the 1 ” 3 transition.

The Sqp s are the atomic transition operators de¬ned in Section 11.1.4, and the labels

q and p range over the values 1, 2, 3.

The form of the dissipative operator Ldis for the three-level ion can be inferred from

the result (18.44) for the two-level atom, by identifying each pair of levels connected

by a decay channel with a two-level atom. For example, the lowering operator σ’ in

eqn (18.44) will be replaced by S13 for the 3 ’ 1 decay channel, and the remaining

transitions are treated in the same way.

There are two important simpli¬cations in the present case. The ¬rst is that the

phase-changing collision term in eqn (18.44) is absent for an isolated ion. The second

simpli¬cation is the assumption that the reservoirs coupled to the three transitions”

i.e. the modes of the radiation ¬eld near resonance”are at zero temperature. This

approximation is generally accurate at optical frequencies, since kB T ωopt for any

reasonable temperature.

One can use these features to show that Ldis is de¬ned by

“31

Ldis ρS = ’ (S31 S13 ρS + ρS S31 S13 ’ 2S13 ρS S31 )

2

“32

’ (S32 S23 ρS + ρS S32 S23 ’ 2S23 ρS S32 )

2

“2

’ (S21 S12 ρS + ρS S21 S12 ’ 2S12 ρS S21 ) . (18.171)

2

¾ The master equation

This expression for Ldis can be cast into the general Lindblad form (18.117) by setting

K = 3 and de¬ning the operators

C1 = “31 S13 , C2 = “32 S23 , C3 = “2 S12 , (18.172)

corresponding respectively to the decay channels 3 ’ 1, 3 ’ 2, and 2 ’ 1.

The Rabi frequency „¦L is small compared to the laser frequency ωL , so the

Schr¨dinger-picture master equation,

o

‚

ρS (t) = [HS , ρS (t)] + Ldis ρS (t) ,

i (18.173)

‚t

involves two very di¬erent time scales, 1/ωL 1/„¦L . Di¬erential equations with this

feature are said to be sti¬, and it is usually very di¬cult to obtain accurate numerical

solutions for them (Press et al., 1992, Sec. 16.6). In the case at hand, this di¬culty

can be avoided by transforming to the interaction picture.

The general results in Section 4.8 yield the transformed master equation

‚I †

ρS (t) = HS1 , ρI (t) + U0 (t) Ldis ρS (t) U0 (t) ,

I

i (18.174)

S

‚t

where U0 (t) = exp (’iHS0 t/ ) and the transform of any operator X is X I (t) =

†

U0 (t) XU0 (t). Applying this rule to the transition operators gives

†

Sqp (t) = U0 (t) Sqp U0 (t) = eiωqp t Sqp ,

I

(18.175)

and this in turn leads to

†

U0 (t) Ldis ρS (t) U0 (t) = Ldis ρI (t) . (18.176)

S

Thus we arrive at the useful conclusion that Ldis has the same form in both pictures.

The transformed interaction Hamiltonian is

HS1 = „¦L S31 e’iδt + HC ,

I

(18.177)

where δ = ωL ’ ω31 . The interaction-picture master equation (18.174) is not sti¬, but

it still has time-dependent coe¬cients. This annoyance can be eliminated by a further

transformation

ρS (t) = eitF ρI (t) e’itF , (18.178)

S

where

F= fq Sqq . (18.179)

q

The algebra involved here is essentially identical to the original transformation to

the interaction picture, and it is not di¬cult to show that the equation of motion

¿

Quantum jumps

for ρ (t) will have constant coe¬cients provided that the parameters fq are chosen to

satisfy

f3 ’ f1 = δ . (18.180)

The simple solution f1 = f2 = 0, and f3 = δ, leads to

‚

ρ (t) = H S1 , ρS (t) + Ldis ρS (t) ,

i (18.181)

‚t S

where the transformed interaction Hamiltonian is

⎡ ¤

0 0 „¦— L

H S1 = ⎣ 0 0 0 ¦ . (18.182)

„¦L 0 ’δ

We are now in a position to calculate all the bits and pieces that are needed for the

direct solution of the master equation (18.181), or the application of the MCWFA. We

leave the algebra as an exercise for the reader and proceed directly to the numerical

solution of the master equation. The density operator for this problem is represented

by a 3 — 3 hermitian matrix which is determined by nine real numbers. Thus the

master equation in this case consists of nine linear, ordinary di¬erential equations

with constant coe¬cients. There are many packaged programs that can be used to

solve this problem.

Of course, this means that we do not really need the MCWFA, but it is still useful

to have a solvable problem as a check on the method. In Fig. 18.7 we compare the

direct solution to the average over 48 trials of the MCWFA. The match between the

averaged results and the direct solution can be further improved by using more trials

in the average, but it should already be clear that the MCWFA is converging on a

solution of the master equation.

Following the general practice in physics, we assume”on the basis of this special

case”that the MCWFA can be con¬dently applied in all cases. In particular, this

includes those applications for which the dimension of the relevant Hilbert space is

large compared to the number of trials needed.

Quantum trajectories—

18.7.5

The results displayed in Fig. 18.7 show that the full-blown master equation”whether

solved directly or by averaging over repeated trials of the MCWFA”does no better

than the rate equations of Section 18.7.1 in describing the phenomenon of interrupted

¬‚uorescence. This should not be a surprise, since the master equation describes the

evolution of the entire ensemble of state vectors for the ion.

What is needed for the description of quantum jumps (interrupted ¬‚uorescence)

is an improved version of the simple on-and-o¬ model used to derive the random

telegraph signal in Fig. 18.3. This is where single trials of the MCWFA come into

play. Each trial yields a sequence of state vectors

|Ψ (t1 ) , |Ψ (t2 ) , . . . , |Ψ (tN ) , (18.183)

which is a discrete sampling of a continuous function |Ψ (t) . This has led to the use

of the name discrete quantum trajectory for each individual trial of the MCWFA.

The master equation

2!

J

Fig. 18.7 The population of |µ3 as a function of time. The smooth curve represents the

direct solution of eqn (18.181) and the jagged curve is the result of averaging over 48 trials of

the Monte Carlo wave function algorithm. Time is measured in units of the decay time 1/“31

for the 3 ’ 1 transition. In these units „¦L = 0.5, δ = 0, “32 = 0.01, and “21 = 0.001.

An example of the upper-level population P3 obtained from a single quantum trajec-

tory is shown in Fig. 18.8. Once again, a judicious choice from the results for several

trajectories nicely exhibits the random telegraph signal characterizing interrupted ¬‚u-

orescence.

According to the standard rules of quantum theory, the information from a com-

pleted measurement”in particular, the collapse of the state vector”should be taken

into account immediately. In the algorithm presented in Section 18.7.3 the new infor-

2!

J

Fig. 18.8 The population of |µ3 as a function of time for a single quantum trajectory. The

parameter values are the same as in Fig. 18.7.

Quantum jumps

mation is not used until the next time step at tn + ∆t, so single trials of the Monte

Carlo wave function method are approximations to the true quantum trajectory.

A more re¬ned treatment involves allowing for the projection or collapse event to

occur one or more times during the interval ∆t, and using the dissipative Hamiltonian

to propagate the state vector in the subintervals between collapses. With this kind

of analysis, it can be shown that the Monte Carlo method is accurate to order ∆t.

Increasing the accuracy to order ∆t2 requires the inclusion of jumps at both ends of

the interval and also the possibility that two jumps can occur in succession (Plenio

and Knight, 1998).

Results like that shown in Fig. 18.8 might tempt one to believe that the Monte Carlo

technique”or the more re¬ned quantum trajectory method”provides a description of

single quantum events in isolated microscopic samples. Any such conclusion would be

completely false. A large sample of trials for the Monte Carlo technique will resemble

a corresponding set of experimental runs, but the relation between the two sets is

purely statistical. Both will yield the same expectation values, correlation functions,

etc. In other words, the Monte Carlo or quantum trajectory methods are still based

on ensembles. The di¬erence between these methods and the full master equation is

that the ensembles are conditioned, i.e. reduced, by taking experimental results into

account.

Quantum state di¬usion—

18.7.6

As explained above, the standard formulations of quantum theory do not apply to

individual microscopic samples, but rather to ensembles of identically prepared sam-

ples. Several of the founders of the quantum theory, including Einstein (Einstein et al.,

1935) and Schr¨dinger (Schr¨dinger, 1935b), were not at all satis¬ed with this feature,

o o

and there have been many subsequent e¬orts to reformulate the theory so that it ap-

plies to individual microscopic objects. One approach, which has attracted a great deal

of attention, is to replace the Schr¨dinger equation for an ensemble by a stochastic

o

equation”e.g. a di¬usion equation in the Hilbert space of quantum states”for an

individual system.

The universal empirical success of conventional quantum theory evidently requires

that the new stochastic equation should agree with the Schr¨dinger equation when ap-

o

plied to ensembles. Many such equations are possible, but symmetry considerations”