6

Chapter 6
Stochastic Methods













Abstract In all physical processes there is an associated loss mechanism. In this chapter we shall consider how losses may be included in the quantum mechanical equations of motion. There are several ways in which a quantum theory of damping may be developed. We shall adopt the following approach: We consider the system of interest coupled to a heat bath or reservoir. We first derive an operator master equation for the density operator of the system in the Schro¨dinger or interaction picture. Equations of motion for the expectation values of system operators may directly be derived from the operator master equation. Using the quasi-probability representations for the density operator discussed in Chap. 4, the operator master equation may be converted to a c-number Fokker–Planck equation. For linear prob- lems a time-dependent solution to the Fokker–Planck equation may be found. In cer- tain nonlinear problems with an appropriate choice of representation the steady-state solution for the quasi-probability distribution may be found from which moments may be calculated.
Using methods familiar in stochastic processes the Fokker–Planck equation
may be converted into an equivalent set of stochastic differential equations. These stochastic differential equations of which the Langevin equations are one example are convenient when linearization is necessary. We begin then with a derivation of the master equation. We follow the method of Haake [1].



Master Equation


We consider a system described by the Hamiltonian HS coupled to a reservoir de- scribed by the Hamiltonian HR. The reservoir may be considered to be a large num- ber of harmonic oscillators as, for example, the modes of the free electromagnetic field or phonon modes in a solid. In some cases the reservoir may be more appro- priately modelled as a set of atomic energy levels. The derivation of the master equation is not dependent on the specific reservoir model. There is a weak inter- action between the system and the reservoir given by the Hamiltonian V . Thus the total Hamiltonian is


93


H = HS + HR + V (6.1)
Let w(t) be the total density operator of the system plus reservoir in the interaction picture. The equation of motion in the interaction picture is


dw(t) dt


i
= h¯ [V (t), w(t)] (6.2)

The reduced density operator for the system is defined by
? (t)= TrR{w(t)} (6.3)
where TrR indicates a trace over reservoir variables. We assume that initially the system and reservoir are uncorrelated so that
w(0)= ? (0) ? ?R (6.4)

where ?R is the density operator for the reservoir.
Integrating (6.2) we obtain


t
i ¸
w(t)= w(0) ? h¯
0



dt1[V (t1), w(t1 )] . (6.5)

Iterating this solution we find



? . i .n ¸t t1

w(t)= w(0)+ ?
n=1

? h¯

dt1
0 0

dt2 •••

tn?1
¸
× dtn[V (t1), [V (t2),... [V (tn), w(0)]]] . (6.6)
0

Performing the trace over reservoir variables

? . i .n ¸t t1


tn?1
¸

? (t)= ? (0)+ ?
n=1

? h¯

dt1
0 0

dt2 •••
0

dtn






where

× TrR{[V (t1), [V (t2),... [V (tn), ?R ? ? (0)]]]}
? (1 + U1(t)+ U2(t)+ ••• )? (0) (6.7)
? U (t)? (0)

. i .n

t t1
¸ ¸

Un(t)=

? h¯

tn?1
¸

TrR
0

dt1
0

dt2 •••

× dtn[V (t1), [V (t2),... [V (tn), ?R ? (•)]] .. .] . (6.8)
0



Thus



= [U?1(t)+ U?2(t)+ ••• ]U (t)?1? (t)
? l(t)? (t) (6.9)


where the generator of time development is
l(t)= [U?1(t)+ U?2(t)+ ••• ]U (t)?1 . (6.10)

We now assume that V (t) is such that

TrR(V (t)?R)= 0 . (6.11)

This ensures that U1(t)= 0. If the perturbation is weak we may drop terms from l(t)
of order higher than two. Thus

l(t)= U?2(t)

t
1 ¸
= ? h¯ 2
0


dt1 TrR[V (t), [V (t1), ?R ? (•)]] . (6.12)


Thus to second order in the perturbation



d? =
dt

t
1 ¸
? h¯ 2
0



dt1



TrR


[V (t), [V (t1), ?R ?



? (t)]] . (6.13)


The next-order correction is at least quartic in the coupling and thus we expect (6.13) to be a good approximation.
Let us now consider the case of a damped simple harmonic oscillator. In this case

V (t)= h¯ (a†?(t)ei?0t + a?†(t)e?i?0t ) (6.14)


where



and


?(t)= ?g j b je?i? jt , (6.15)
j


[b j , b†]= ? jk . (6.16)


Substituting (6.14) into (6.13) we find that the following integrals are required


t
¸
I1 =
0
t
¸
I2 =
0


dt1(?(t)?(t1))ei?0 (t+t1 ) , (6.17)


dt1(?†(t)?†(t1))e?i?0(t+t1 ) , (6.18)










which we now evaluate.


t
¸
I3 =
0
t
¸
I4 =
0


dt1(?(t)?†(t1))ei?0 (t?t1 ) , (6.19)


dt1(?†(t)?(t1))e?i?0 (t?t1 ) (6.20)

Using the definition of ?(t) we have

t

¸
I1 =
0

dt1 ?gig j(bib j )Re?i(?it+? jt1 )ei?0(t+t1 ) . (6.21)
i, j


Converting the sum over modes to a frequency-space integral


t
¸
I1 =
0


?
¸
dt1
0


?
d?1 ¸
2? ? (?1)
0


d?2
2? ? (?2)g(?1)g(?2)(b(?1)b(?2))R

× e?i(?1t+?2t1 )+i?0(t+t1 ) (6.22)

where ? (?) is the density of states function.
For a thermal bath the phase dependent correlation function (b(?1)b(?2))R = 0.
However, certain specially prepared reservoirs such as squeezed reservoirs may have
phase-dependent correlations.
In order to include these we now assume that
(b(?1)b(?2)) = 2?M(?1 )? (2?0 ? ?1 ? ?2) (6.23)

which corresponds to a multimode squeezed vacuum state with the carrier frequency equal to the cavity resonance frequency. Thus


t
¸
I1 =
0


?
¸
dt1
0


d?
2? ? (?)? (2?0 ? ?)g(?)g(2?0 ? ?)M(?)e


i(?0??)(t?t1 )



. (6.24)


Note that the time integral depends only on t ? t1. This suggests the change of vari- able ? = t ? t1 and thus


t
¸
I1 =
0


?
¸ d?
d? 2? ? (?)? (2?0 ? ?)g(?)g(2?0 ? ?)M(?)e
0


i(?0??)?



. (6.25)


We now make the first Markov approximation by assuming ? (?), g(?) and M(?)
are slowly varying functions around ? = ?0, where ?0 is very large. Thus it is convenient to make the change of variable ? = ? ? ?0 and write


t ?
¸ ¸ d? 2 2

?i??

I1 ?
0

d? 2? ? (? + ?0)g (? + ?0)M(? + ?0)e
??

, (6.26)

assuming a symmetry around ?0. Since the integral over frequency is assumed to be a rapidly decaying function of time we have extended the upper limit of the time integration to infinity. Interchanging the order of the time and frequency integral this term becomes

?
¸
I1 ?
??


d? 2? ?


2(? + ?0)g2(? + ?0)M(? + ?0)

. . 1 ..
?? (?) ? i PV ?


(6.27)

where we have used



?
d? e±i?? = ?? (?) ± i PV
0


. 1 .
?




(6.28)

with PV being the Cauchy principal value part defined by


b
PV ¸ f (?)



= lim

? ??
¸


f (?)


b
¸
d? +

?
f (?) d?



. (6.29)

? ??0 ? ? ? ?
?a ?a ?
If we now define the damping rate ? by
? ? ? 2(?0)g2(?0) (6.30)


and a term



?
?¯ = PV ¸



d? 1
? (? + ? )g (? + ? )M(? + ? ) (6.31)



then

2? ? 0 0 0
??

?

I1 = 2 M(?0 )+ i?¯ . (6.32)
Proceeding in a similar way we find

I = ? M? (? ) ? i?¯ , (6.33)
I = ? (N(? )+ 1) ? i? , (6.34)
I = ? N(? ) i?? , (6.35)
2

where the function N(?) is defined by
(b†(?)b(??)) = 2?N(?)? (? ? ??)

and is thus proportional to the intensity spectrum of the reservoir. The term ? is defined by


?
¸
? = PV


d? 1 2? ?



? 2(?0 + ?)g2(?0 + ?)(N(?0 + ?)+ 1) . (6.36)

??


This term represents a small shift in the frequency of the oscillator. In the case where the system is a two-level atom, this term contributes to the Lamb shift which we discuss in Chap. 10. In what follows we ignore the effects of ? and ?¯ .
Substituting the above results into (6.13) we find that the evolution of the sys- tem’s density operator in the interaction picture is described by the master equation


d? = ? (N + 1)(2a? a†


a†a?


? a†a)+ ? N(2a†? a


aa†?


? aa†)

dt 2

? ? 2 ? ?

+ ? M(2a†? a† ? a†a†? ? ? a†a†)+ ? M?(2a? a ? aa? ? ? aa) . (6.37)
2 2

(For convenience, we have suppressed the functional dependence of N(?0) and
M(?0 )).
If the bath is in thermal equilibrium at temperature T, M = 0 and
N(?0 )= (eh¯ ?0/kT ? 1)?1 (6.38)

which is just the mean number of bath quanta at frequency ?0. In this case the master equation considerably simplifies. If the bath temperature is zero, N = 0, and the master equation simplifies further. In general the positivity of the density operator
2

requires, |M|

? N(N + 1).

Equations of motion for the expectation values of system operators may be di-
rectly derived from the master equation, (6.37). For example, the mean amplitude of the simple harmonic oscillator is given by
d(a) = Tr .a d?ˆ . = ?


dt dt

? 2 (a) (6.39)



which has the solution


(a(t)) = (a(0)) e??t/2 . (6.40)

(Note in the Schro¨dinger picture the mean amplitude evolves as (a(t)) = (a(0))e?i?t
e??t/2 ). Thus the mean amplitude decays at a rate ?/2. The mean number of quanta
(n) = (a†a) obeys the equation
d(a†a)


dt The solution to this equation is

= ??(a†a) + ?N . (6.41)

(n(t)) = (n(0))e??t + N(1 ? e??t ) . (6.42)
In the steady state (n(t))? N and the mean number of quanta in the oscillator is equal to the mean number of quanta in the reservoir at that temperature. The role of
the terms multiplied by M can be seen by evaluating the equation of motion for (a2),
d 2 2
dt (a ) = ??(a ) + ?M . (6.43)

Thus these terms lead to a driving force on the second-order phase dependent moments.
The master equation (6.37) applies to a free damped harmonic oscillator in the interaction picture. If the harmonic oscillator is perturbed by an additional interac- tion HI, the master equation in the interaction picture becomes

d? = l ?

† † †

dt ? h¯ [HI , ? ]+ 2 (N + 1)(2a? a
+ ? N(2a†? a ? aa†? ? ? aa†)

? a a? ? ? a a)

(6.44)

+ ? M(2a†? a† ? (a†)2? ? ? (a†)2)+ ? M? (2a? a ? a2? ? ? a2) .
2 2

In general, the equations of motion for the mean amplitude, mean quantum number etc. are not as easily obtained from (6.44), as they were for the free damped har- monic oscillator. To proceed in such situations, it is desirable to convert the master equation to an equivalent c-number partial differential equation. We now discuss various ways this may be done.



Equivalent c-Number Equations


Photon Number Representation


The operator master equation may be converted into an equation for the matrix elements of ? in the number state basis:


? ?mn = ?N[(nm)1/2?
? t

1
m?1, n?1 ? 2 (m + n + 2)?mn].
1

+ ?(N + 1){[(m + 1)(n + 1)]1/2?m+1, n+1 ?

2 (m + n)?mn}

?
? 2 M

,
2[m(n + 1)]

1/2

?m?1, n+1 ? ,(n + 1)(n + 2)?m,n+2

?,m(m ? 1)?m?2, n

?
? 2 M?

,2[n(m + 1)]

1/2


?m+1,n?1
,

?,(m + 1)(m + 2)?m+2, n ? ,n(n ? 1)?m,n?2

(6.45)


where ?m, n ? (m|? |n). This gives an infinite hierarchy of coupled equations for the
off-diagonal matrix elements. When M = 0 the diagonal elements ?m, m are coupled
only amongst themselves and not coupled to the off-diagonal elements. In this case
the diagonal elements satisfy


dP(n) dt

= t+(n ? 1)P(n ? 1)+ t?(n + 1)P(n + 1) ? [t+(n)+ t?(n)]P(n) (6.46)


where we have set P(n)= (n|? |n) and defined the transition probabilities
t+(n) ? ?N(n + 1) , (6.47)
t?(n) ? ?(N + 1)n . (6.48)

In the steady state the detailed balance condition holds:
t?(n)P(n)= t+(n ? 1)P(n ? 1) (6.49)

and the steady state solution is found by iteration
n t+ (k ? 1)

Pss(n)= P(0) ?
k=1

t (k) . (6.50)


Thus the steady state solution for the ordinary damped harmonic oscillator
(M = 0) is

1
Pss(n)= 1 + N

. N .n
1 + N


. (6.51)

An optical cavity damped into a reservoir with phase-independent correlation func- tions has a power law photon number distribution of thermal light.
In the more general case M ?= 0, or when there are additional terms in the master
equation such as linear driving with the Hamiltonian H0 = h¯ [?(t)a† + ??(t)a], the
coupling of diagonal and off-diagonal matrix elements makes the photon number
representation less convenient for determining ? (t).



P Representation


An operator master equation may be transformed to a c-number equation using the Glauber–Sudarshan representation for ? . It is necessary to first establish the rules for converting operators to an equivalent c-number form. We know the relations
a|?) = ?|?) , (6.52)
(?|a† = ?? (?| . (6.53)
To derive other relations it is convenient to use the Bargmann state " ?) defined by




so that

" ?) = e1/2|?| |?) (6.54)



n

a† " ?) = ? ?

?n + 1|n + 1)

n n!
= ?
?? " ?) . (6.55)



Similarly


?
(? " a = ??? (? " . (6.56)

Hence, given the P representation
¸ 2



we find

? = d2? " ?)(? " e?|?|

P(?) (6.57)

a†? = ¸ d2? ? ( ?

? )e???? P(?) (6.58)



and integrating by parts

?? "

)( "

2 . .

a†? = ¸ d2?

" ?)(?

" e?|?|


? ?? ? ??


P(?) . (6.59)


We can thus make an operator correspondence between a† and ?? ? ? /??. A sim- ilar formula holds for a. Summarizing we have the following operator correspon-
dences:

a? ? ?P(?),
. ? .

a†? ?

?? ? ??

P(?),

. ? .

? a ?

? ? ???

P(?),

? a† ? ??P(?) . (6.60)

Consider the correspondences for operator products


a†a? ?

. ? .
?? ? ??



?P , (6.61a)

.
? a†a ? ? ?

?
???

.
??P . (6.61b)


Notice that the order of the operators in (6.61b) reverses, since acting on ? , they operate from the right, whereas on P, they operate from the left.
Note that ? and ?? are not independent variables. In terms of real variables we
may write

? = x + iy,
?? = x ? iy,
? 1 . ? ? .
= ,

?? 2

? x ? i ? y

? 1 . ?
=

+ i ? . . (6.62)

???

2 ? x ? y


To obtain a c-number equation we substitute the P representation for ? into the master equation and use the operator correspondences. This leads to the equation

2 2 2

? P(?) = . 1

. ?

? .
?

? . ?
?

? . ? .
+ ?N P(?) .

? t 2 ?

? + ?
?? ???

+ 2 M ???2 + M ??2

?????


(6.63)

When M > N this equation has “non positive-definite” diffusion, hence the P repre- sentation is unable to describe the system in terms of a classical stochastic process.
Alternative representations will be discussed later in this chapter. When M ? N
(6.63) has the form of a Fokker–Planck equation. We shall discuss some useful
properties of Fokker–Planck equations below.




Properties of Fokker–Planck Equations


A general Fokker–Planck equation in n variables may be written in the form


? P(x)= .


? 1 ?
A (x)+

? .
D (x)


P(x) . (6.64)

? t ? ? x j j

2 ? xi ? x j ij


The first derivative term determines the mean or deterministic motion and is called the drift term, while the second derivative term, provided its coefficient is positive definite, will cause a broadening or diffusion of P(x, t) and is called the diffusion term. A = (A j ) is the drift vector and D = (Dij ) is the diffusion matrix. The different
role of the two terms may be seen in the equations of motion for (xk ) and (xk xl ).


d(xk) =
dt

(Ak)


, (6.65)

d(xkxl ) =
dt

(xkAl ) +

(xl Ak)

1
+ 2 Dkl

+ Dlk)


. (6.66)


We see that Ak determines the motion of the mean amplitude whereas Dlk enters into the equation for correlations.
Thus, for the damped harmonic oscillator described by (6.63)


d(?)P =
dt

?
? 2 (?)P , (6.67)

d(???)P =
dt

??(??

?)P + ?N , (6.68)


which are equivalent to (6.39 and 6.41) derived directly from the master equa- tion (6.37). Note that the expectation values ()P are defined by integrals over
P(?, t).


Steady State Solutions – Potential Conditions


For many problems in nonlinear optics it is sufficient to know the steady state solu- tion. That is the solution after all transients have died out. We shall therefore seek a steady state solution to (6.64).
In the steady state we set the time derivatives to zero which gives


? . 1 ?
Ai(x)P(x)+

.
Dij (x)P(x)


= 0 . (6.69)

? xi ?

2 ? x j


As a first attempt consider

1 ?

Ai(x)P(x)= 2 ? x

[Dij (x)P(x)] (6.70)



which implies


D ? ln P
ij ? x j


= 2A (x) ? Dij
i ? ? x j



(x) . (6.71)

Denoting P(x)= exp[?? (x)] we wish to solve


?? (x) = 2(D?1)

.
A (x)

1 ? D jk .


F (x) . (6.72)

? ? xi

ij j

? 2 ? xk ? i


If we consider Fj(x) as a generalized force, ? (x) corresponds to a potential. The system of equations (6.72) can be solved by integration if the so called potential conditions are satisfied

?? 2?
? xi? x j


= ? Fj
? xi


= ? Fi
? x j


? 2?
=
? x j ? xi


. (6.73)


These conditions say that the function ? is well behaved and that the multivariable integral is independent of the path of integration. The potential conditions are always satisfied in the one dimensional case.
Provided the potential conditions are satisfied a steady state solution of the form
P(x)= N exp[?? (x)] (6.74)


exists where



x
¸
? (x)=



2[D(x)?1]ij


.
?A j(x)+


1 ? D jk .




dxi .

2 ? xk
0
The turning points of the potential ? correspond to the values x such that for each
j = 1, ... , n
. 1 ? D jk (x) .

A j(x) ? 2


? xk

= 0 . (6.75)


In systems where the diffusion matrix is diagonal and constant (Dij = W ?ij ) (6.72) become

?? (x)
? ? xi

= 2Ai(x)
W

. (6.76)

Hence the turning points of the potential ? correspond exactly to the determinis- tic steady state solutions, that is the steady state solutions to the first-order moment equations.



Time Dependent Solution


In the case where the drift term is linear in the variable (x) and the diffusion coef- ficient is a constant, a solution to the Fokker–Planck equation may be found using the method of Wang and Uhlenbeck [2]. We consider the Fokker–Planck equation


? P =


? 1
a (x P)+ d


? 2P



. (6.77)

? t ? ?

i ? xi i

2 ij ? xi? x j


The Greens function solution to this equation given by the initial condition
P(xi, 0)= ? (xi ? x0)


is


P(xi, x0,t)=





exp


.
? ??ij (t)?1{[xi ? x0 exp(ait)]

i ?n/2[det ?ij (t)]1/2 i
×[x j ? x0 exp(a jt)]}. (6.78)

where

?ij (t)= ?

2dij

{1 ? exp[(ai + a j)t]} .

ai + a j
The solution for a damped harmonic oscillator initially in a coherent state with
P (?, 0)= ? 2(? ? ?0) is



P(?,t)=


1
?N(1 ? e??t )



exp

?|? ? ?0e??t/2|2
N(1 ? e??t )



. (6.79)


This represents an initial coherent state undergoing relaxation with a heat bath. Its coherent amplitude decays away and fluctuations from the heat bath cause its P function to assume a Gaussian form characteristic of thermal noise. The width of the distribution grows with time until the oscillator reaches equilibrium with the heat bath.
From the above solution we may construct solutions for all initial conditions
which have a non-singular P representation. It is not, however, possible to construct


the solution for the oscillator initially in a squeezed state since no non-singular P function exists for such states. Nor can we find the solution for an oscillator damped into a squeezed bath.
We now consider alternative methods of converting the operator master equa-
tion to a c-number equation, which can be used for initial squeezed states or a squeezed bath.



Q Representation


A c-number representation for the Q function is obtained by first normally ordering all operator products. We shall make use of the following theorem [3].
If f (a, a†) is a function which may be expanded in a power series in a and a†,
then


[a, f (a, a†)] = ? f
? a†
? f
[a , f (a, a )] =
? a


, (6.80a)

. (6.80b)


The proof of these relations is as follows: We assume that we may expand f in antinormal order f (a)
[a†, f ]= ? f (a)[a†, ar (a†)s] . (6.81)
r,s

Using the following result for the commutators

[A, BC]= [A, B]C + B[A, C] (6.82)

where A, B and C are noncommuting operators, we may write

[a†, f ]= ? f (a){[a†, ar ]a†s + ar[a†, a†s]}
r,s
r,s rar?1a†s (6.83)
r,s
= ? f . (6.84)
? a

The proof of (6.80a) follows in a similar way.
We consider as an example the term
? a†a = a†? a ? [a†, ? ]a . (6.85)


Using the result above



? a†a = a†? a + ?? a . (6.86)
? a


Taking matrix elements in coherent states yields

. . . .

(?|? a†a|?) = ?

? .a†? + ?? . ?

. ? a .
. .

= ? 2 + ? ?
??


Q (?) . (6.87)


Following this procedure we may convert the master equation for the damped har- monic oscillator into a c-number equation for the Q function


? Q = ? . ?


? + ?


?? Q + ? M? ?


+ M ?

+ 2(N + 1) ? . Q .

? t 2 ? ?

? ??

2 ? ??2

? ?2

? ?? ??


(6.88)

This differs from the corresponding equation of motion for the P function only through the phase independent diffusion coefficient which is N + 1 rather than N. This is sufficient to give a positive definite diffusion matrix when the bath is in an ideal squeezed state.
To illustrate the use of the Q function consider a damped oscillator which is initially in the squeezed state |?0, r). Using the Wang and Uhlenbeck solution for
an initial ? -function and convoluting this with the Gaussian Q function for an initial squeezed state we arrive at the result

1 1 T ?1

Q(?, t)=

2? ,det ? (t)

exp[? 2 u(t) ?

(t)u(t)] (6.89)



where


. ? ? e??t/2 .
u(t)=

?? ? ??e??t/2
. ?sinh 2r cosh 2r + 1. e??t


. M N + 1.



??t

? (t)=

cosh 2r + 1 ?sinh 2r

2 + N + 1 M?

(1 ? e ).


The variances for the quadrature phase operators Xi for the oscillator are then easily found to be


1 2r
V (X1)= 4 [(e?

? 1)e


??t

+ 2(N + Re{M})(1 ? e


??t


)+ 1] , (6.90a)

1 2r
V (X2)= 4 [(e


+ 1)e

??t

+ 2(N + Re{M})(1 ? e

??t


)+ 1] . (6.90b)


In Fig. 6.1 we depict the evolution of an initial squeezed state coupled to a zero temperature reservoir (N = M = 0). The amplitude of the squeezed state damps to zero and the variances in X1 and X2 become equal at the value one corresponding to the vacuum.


Fig. 6.1 Evolution of the error contour of an initial squeezed state of a simple harmonic oscillator damped into a zero temperature heat bath














Wigner Function

Alternatively one may convert the operator master equation into a c-number equa- tion via the Wigner function. This is best accomplished by deriving an equation for the characteristic function



where

? (? )= Tr{D? } (6.91)

D = e? a†?? ?a .


An equation of motion for ? (? ) may be derived as follows

??(? ) = Tr .
? t

?? .
? t


. (6.92)


To illustrate the technique we shall derive the equation of motion for the Wigner function of a damped harmonic oscillator.
We require some operator rules to convert to differential operators. Writing D in
normal order
D = e??? ?/2e? a† e?? ?a , (6.93)

? ? ?
D = ?


D + a†D , (6.94)

?? 2

or
. ?
a†D = +


? ? .




D . (6.95)



Similarly we may show





Da =

??

. ??
2 ?

2



?
?? ?



.
D . (6.96)



Writing D in antinormal order



D = e?? ?/2e?? ?ae? a† . (6.97)



Thus


or


? D = ? ? D + Da† (6.98)
?? 2

and similarly

. ?
Da† =
??

? ? .
? 2


D (6.99)

. ?
aD =
2

?
? ?? ?

.
D . (6.100)

Then using these rules the master equation (6.37) yields the following equation for the characteristic function

??(? ) = ? . ? ?

? . 2

? t 2

?|? | ? ?

?? ? ? ? ??

? (? ) ? ?N|? | ? (? )

?M 2

?M? 2

(? ?) ? (? )
2

2 ? ? (? ) . (6.101)

The equation for the Wigner function is obtained by taking the Fourier transform of this equation since

Thus

¸
W (?)=

e? ?????? ? (? )d2? . (6.102)

¸ ¸ ? ? ? ?

e? ?????? ? ??? (? )d2? = ?

?? ???
? 2W (?)

(e? ??? ? )? (? )d2?

and

= ? ?? ??? , (6.103)

¸ e? ?????? ? ? ?

2 ? ¸

? ?????? ? 2

?? ? ? (? )d ? = ?? (e

) ? (? )d ?
?? ?

= ? ¸ ? (? ) ? (e
?? ?? ?
?

? ??????


)d2? (6.104)

= ? ?? [?W (?)] . (6.105)

Using these results we may write the equation for the Wigner function as


?W (?)

? . ?
= ? +

? ??. W

? t 2 ??

???
2 2 2

+ ? .M? ?

? . 1 . ?

.
W . (6.106)

2 ???2 + M ??2 + 2 N + 2

?? ???


A comparison of the three equations (6.63, 6.88 and 6.106) for the P, Q and Wigner functions show that they differ only in the coefficient of the diffusion term being ?N, ?(N + 1) and ?(N + 1/2), respectively. However, the additional +? and +?/2 in the equations for the Q and Wigner function are sufficient to ensure that these equations have positive definite diffusion.



Generalized P Representation


In our study of nonlinear problems we shall find systems which either do not give Fokker–Planck equations in the Q and Wigner representations or no steady state solution may readily be found. For some systems a steady state solution in terms of a Glauber–Sudarshan P representation does not exist. For such systems the complex P representation is sometimes useful in deriving a steady state solution to Fokker– Planck equations. The positive P representation is useful when it is desirable to have a Fokker–Planck equation with a positive definite diffusion term, as is necessary in order to deduce the corresponding stochastic differential equations.
Master equations may be converted to a c-number representation using the com-
plex P representation by an analogous set of operator rules used for the diagonal P
representation.
The nondiagonal coherent state projection operator is defined as

?(?)= |?)(? ?|
(? ?|?)
where (?) denotes (?, ? ). The following identities hold
. .


(6.107)

a?(?)= ??(?), a†?(?)=

? + ?
??


?(?),

. ?
?(?)a† = ?(?)? , ?(?)a =
??

.
+ ? ?(?) . (6.108)


By substituting the above identities into (4.65) defining the generalized P repre- sentation, and using partial integration (providing the boundary terms vanish) these identities can be used to generate operations on the P function depending on the representation.

a) Complex P representation



a? ? ?P(?), a†? ?

. ? .
? ? ??


P(?),

. ? .
? a† ? ? P(?), ? a ? ? ?


P(?) . (6.109)


This procedure yields a very similar equation to that for the Glauber–Sudarshan P
function. We assume that, by appropriate reordering of the differential operators,


we can reduce an operator master equation to the form [where (?, ? ) = (?) = (?(1), ?(2))]


?? = ¸ ¸
? t
c c?
¸ ¸
=

?(?) ? P(?) d? d?
? t
.
d?(1) d?(2) P(?)





A? (?) ?





+ D?? (?)



? . ?(?) .


c c?

?? ? 2

?? ? ???



(6.110)



We now integrate by parts and if we can neglect boundary terms, which may be made possible by an appropriate choice of contours, c and c?, at least one solution is
obtained by equating the coefficients of ?(?)

? P(?) . ?
= ?


A? (?)+ 1 ?

? .
D?? (?)



P(?) . (6.111)

? t ?? ?

2 ?? ? ???


This equation is sufficient to imply (6.110) but is not a unique equation because the
?(?) are not linearly independent. The Fokker–Planck equation has the same form as that derived using the diagonal P representation with ?? replaced by ? .
It should be noted that for the complex P representation, A? (?) and D?? (?) are always analytic in (?), hence if P(?) is initially analytic (6.111) preserves this analyticity as time develops.

b) Positive P Representation

The operator identities for the positive P representation are the same as (6.109) for the complex P representation. In addition, using the analyticity of ?(?, ? ) and noting that if



then

? = ?x + i?y, ? = ?x + i?y,




and

? ?(?)=
??


?

?
? ?x


?

?(?)= ?i

?
? ?y


?


?(?)

?(?)=
??


? ?x

?(?)= ?i


? ?y

?(?) . (6.112)

Thus in addition to (6.109) we also have


a†? ?

. ? .
? ? ? ?


P(?) ?

? + i ? .
? ?y



P(?),

. ? .

. ? .

? a ?

? ? ??

P(?) ?

? + i
? ?y

P(?) . (6.113)


The positive P representation may be used to give a Fokker–Planck equation with a positive definite diffusion matrix. We shall demonstrate this in the following.


We assume that the same equation (6.64) is being considered but with a positive P
representation. The symmetric diffusion matrix can always be factorized in the form

D(?)= B(?)BT(?).

We now write

A (?)= Ax(?)+ iAy(?) , (6.114)
B (?)= Bx(?)+ iBy(?) , (6.115)

where Ax, Ay, Bx, By are real. We then find that the master equation yields
?? = ¸¸ d2? d2? ?(?)(? P(?)/? t)
? t
= ¸¸ P(?)[A? (?)? x + A? (?)? y + 1 (B?? B?? ? x ? x + B?? B?? ? y ? y

x ?
+ 2B?? B?? x y

y ? 2 x
2 2

x ? ?

y y ? ?

x y ?? ?? )]?(?)d ? d ? . (6.116)

Here we have, for notational simplicity, written ? /?? ? = ? x etc., and have used the

x ?
analyticity of ?(?) to make either of the replacements
? /?? ? ? ? x ? ?i? y




(6.117)

? ?

in such a way as to yield (6.116). Now, provided partial integration is permissible, we deduce the Fokker–Planck equation:


? P(?)/? t = [?? x A? (?) ? ? y A? (?)+


[? x ? xB?? (?)B?? (?)

? x ? y

2 ? ? x x

+ 2? x ? yB?? (?)B?? (?)+ ? y ? yB?? (?)B?? (?)]} P(?) (6.118)

? ? x y

? ? y y

Again, this is not a unique time-development equation but (6.116) is a consequence of (6.118).
However, the Fokker–Planck equation (6.118) now possesses a positive semidef- inite diffusion matrix in a four-dimensional space whose vectors are

(?(1)

(2)

(1)

(2)

x , ?x , ?y , ?y ) ? (?x, ?x, ?y, ?y . (6.119)
We find the drift vector is

(1)

(2)

(1)

(2)

A (?) ? (Ax (?), Ax (?), Ay (?), Ay (?)) , (6.120)
and the diffusion matrix is

? BxBT
D (?)= ?
ByBT

BxBT ?
? (?) ? B(?)BT(?) (6.121)
ByBT



where

.Bx 0.

B(?)=

By 0

(?) (6.122)


and D is thus explicitly positive semidefinite.



6.3 Stochastic Differential Equations


A Fokker–Planck equation of the form


? P =

? 1
A (x, t)P +

? ? [B(x, t)BT(x, t)]


P (6.123)

? t ? ? ? xi i

2 ? ? xi ? x j ij


may be written in a completely equivalent form as

dx
dt = A(x, t)+ B(x)E(t) (6.124)

where E(t) are fluctuating forces with zero mean and ? correlated in time
(Ei(t)E j (t?)) = ?ij ? (t ? t?) . (6.125) We have written (6.124) in the form of a Langevin equation. The relationship be-
tween (6.123 and 6.124) may be derived more rigorously in terms of stochastic dif- ferential equations where Ito’s rules are used. However, the relation quoted in (6.123 and 6.124) is sufficient for our use. The reader is referred to the texts C.W. Gardiner for a complete discourse on stochastic differential equations and their applications to quantum noise problems.
We shall illustrate the use of the stochastic differential equation for a particle undergoing damping and diffusion in one dimension. This motion is described by the Fokker–Planck equation


? P(x) = ? ? [xP(x)] + D ?



P(x) (6.126)

? t ? x

2 ? x2


where ? is the damping coefficient and D is the diffusion coefficient. This equa- tion describes an Ornstein–Uhlenbeck process. It may, for example, describe the Brownian motion of a particle under the random influence of collisions from many particles in thermal motion where the variable x represents the particle’s velocity.
The Langevin equation equivalent to the Fokker–Planck Equation (6.126) is

x? = ?? x + ?


(t) (6.127)


where


(E(t)E(t?)) = ? (t ? t?) .



The solution to this equation is

t
x(t)= x(0)e??t + ? ¸

0




e??(t?t?)E(t?)dt? . (6.128)


If the initial condition is deterministic or Gaussian distributed, then x(t) is clearly Gaussian, with mean and variance
(x(t)) = (x(0))e??t , (6.129)
2
? ¸t ? .

Var[x(t)] =

?[x(0) ? (x(0))]e??t + ?

e??(t?t? )E(t?)dt??
?

. (6.130)

? 0

Assuming the initial condition is independent of E(t), we may write


t
Var{x(t)} = Var{x(0)}e?2?t + D
0


e?2?(t?t?)dt?

. D .
= Var{x(0)}? 2?

e?2?t

D
+ . (6.131)
2?



In the steady state


D
Var{x(t)} = 2? . (6.132)

The two time correlation function may be calculated directly, as follows
(x(t), x(s)) = (x(t)x(s))? (x(t))(x(s))


= Var{x(0)}e??(t+s) + D

.¸t


0

s
e??(t?t?)E(t?)dt? ¸

0

.
e??(s?s?)E(s?)ds?


= Var{x(0)}e??(t+s) + D

min(t,s)
e??(t+s?2t?)dt?

0

. D .


?(t+s)

D ??|t?s|

= Var{x(0)}? 2? e? + 2? e

. (6.133)



In the stationary state


D
(x(t), x(s)) = 2? e?


?|t?s|



. (6.134)

We shall now consider the equivalent Langevin equation for the Fokker–Planck equation for the damped harmonic oscillator. The Fokker–Planck equation for the P representation is

? P ? . ?
= ? +

? ??. P + ?N ?



P . (6.135)

? t 2 ??

???

?? ???


Note that we have set M = 0, as the P representation cannot be used with squeezed baths since the diffusion matrix is nonpositive definite. In such cases the Q-function could be used.
Equation (6.135) is an example of an Ornstein–Uhlenbeck process. The diffusion
matrix is




which may be factored as


D = ?N

. 0 1 .
1 0


(6.136)



where

D = BBT
. ?N .1/2 . i 1 .
B =




. (6.137)

2 ?i 1
Thus the stochastic differential equations become

d . ? . = . ?

.. . . . .. .
+

dt ??

2 0 ?
?? ?
2

?N i 1
2 ?i 1

?1(t)
?2(t)

(6.138)

where ?1(t) and ?2(t) are independent stochastic forces which satisfy
(?i(t)? j (t?)) = ?ij ? (t ? t?) . (6.139) Equation (6.138) may be written


d? =
dt

? 2 ? + ,?N?(t),




where

d?? =
dt

? 2 ? + ,?N??(t) (6.140)

1

?(t) ? ?2 [?2(t)+ i?1(t)]
is a complex stochastic force term which satisfies
(?(t)??(t?)) = ? (t ? t?) .

An alternative factorisation is
. ?N .1/2 . ei?/4 e?i?/4 .
B =










. (6.141)

2 e?i?/4 ei?/4


In this case


where


d? =
dt

? 2 ? + ,?N??



(6.142)

?? = 1

(? ei?/4 + ? e?i?/4 . (6.143)

?2 2 1
One easily verifies that (?? ?? ?) = ? (t ? t?).


The solutions derived from (6.140) are

?
(?(t)) = (?(0)) exp . t. , (6.144)
(??(t)?(t)) = (??(0)?(0))e??t + N(1 ? e??t ),
(?2)ss = (??2)ss = 0,
(???)ss = (???)ss = N , (6.145)

where ss denotes steady state.



6.3.1 Use of the Positive P Representation


The relationship (6.123 and 6.124) between the Fokker–Planck equation and the stochastic differential equation only holds if the diffusion matrix

D (x, t)= B(x, t)BT(x, t)

is positive semidefinite. In some cases use of the Glauber–Sudarshan P representa- tion will result in Fokker–Planck equations with a non-positive semi-definite diffu- sion matrix, for example, if the bath is squeezed. In such cases use of the positive P representation will give a Fokker–Planck equation with a positive semi-definite diffusion matrix



where

D (?)= B(?)BT (?) (6.146)
. Bx 0 .



and ? = (?, ? ).

B(?)=

By 0

The corresponding stochastic differential equations may be written
d . ?x . = . Ax(?) . + . Bx(?)E(t) . (6.147)

dt ?y

Ay(?)

By(?)E(t)


on recombining real and imaginary parts

d? = A(?)+ B(?)E(t) . (6.148)
dt
Apart from the substitution ?? ? ? , (6.148) is just the stochastic differential equa- tion which would be obtained by using the Glauber–Sudarshan P representation, and naively converting the Fokker–Planck equation with a non-positive definite diffusion matrix into a stochastic differential equation. In the above derivation the two formal
variables (?, ??) have been replaced by variables in the complex plane (?, ? )
that are allowed to fluctuate independently. The use of the positive P representation justifies this procedure.


6.4 Linear Processes with Constant Diffusion


For linear processes with constant diffusion coefficients a number of useful results may be proven. These may be derived from the Fokker–Planck equation using the solution (6.78) or the equivalent Langevin equation. We shall quote the results here. The Langevin equations are a useful starting point since for nonlinear processes approximate results may be obtained by a linearization procedure. We consider a process described by the Langevin equation


dx(t) dt

= ?Ax(t)+ BE(t) (6.149)


where A and B are constant matrices. This describes a multivariate Ornstein– Uhlenbeck process. Suppose AAT = ATA, then we can find an orthogonal matrix S such that

SST = 1
SAST = SATST = Diag{?1, ?2 ... ?n} . (6.150)

The two time correlation function is given by
(x(t), xT(s)) = STG(t, s)S


where



(SBBTST)ij



? t s



? t ? s

[G(t, s)]ij =


?i + ? j

(e?

i| ? | ? e? i ?

j ) . (6.151)


In the stationary state the second term in the parentheses is zero and the correlation is only a function of the difference, ? ? t ? s.
Let us define the stationary covariance matrix ? by
? = (xss(t), xT (t)) . (6.152)

Then by setting dx(t)/dt = 0 in (6.149) we find that ? obeys the equation

A? + ? AT = BBT . (6.153)

In the case of a two dimensional problem it may be shown that
(Det A)BBT + [A ? (Tr A)I]BBT[A ? (Tr A)I]T
? = 2(Tr A)(Det A) . (6.154)

The two time correlation function in the steady state may be shown to obey the same time development equation as the mean. That is

d
d? [Gss(?)] = ?AGss(?) . (6.155)

6.5 Two Time Correlation Functions in Quantum Markov Processes 117
The computation of Gss(?), therefore requires the knowledge of Gss(0)= ? and the time development equation of the mean.
It is often of more interest to view the noise in the frequency domain. We are thus lead to define the noise spectrum by,


?
1 ¸
S(?)= 2?
??

Using (6.155 and 6.153) we find


e?i??



Gss(?)d? . (6.156)


1 1 T T ?1

S(?)= 2? (A + i?)? BB (A

? i?)

. (6.157)



6.5 Two Time Correlation Functions in Quantum Markov Processes


We shall now demonstrate how two time correlation functions for operators may be derived from the master equation or the equivalent Fokker–Planck equation.
Consider a system coupled to a reservoir. W (t) is the total density operator in the
Schro¨dinger picture and H is the Hamiltonian, A and B are operators for variables to be measured, then



and

(A(t)) = Tr{AW (t)} (6.158)

(A(t + ?)B(t)) = Tr eiH ?/h¯

Ae?iH ?/h¯

BW (t .

(6.159)


while this is exact it is not particularly useful. For a system interacting with a heat bath in the Markov approximation we wish to express everything in terms of the Liouvillian for the reduced system in which heat bath variables have been traced out.
Supposing A and B are operators in the system space, then
(A(t + ?)B(t)) = Trs{A TrR[e?iH ?/h¯ BW (t)eiH ?/h¯ ]} . (6.160) The equation of motion for the term

X (?,t)= e?iH ?/h¯ BW (t)eiH ?/h¯


(6.161)



in terms of ? is


ih¯ ?
??



X (?,t)= [H , X (?,t)] . (6.162)

Proceeding in exactly the same way as for the derivation of the Markovian master equation (6.37) which may be written in the form

?
? t ? (t)= L? (t) (6.163)


where L is a Liouvillian operator, where ? = TrR{W } is the reduced density operator
for the system, we may derive the equation


? [Tr
??


R{X (?,t)}]= L{TrR


X (?,t)} (6.164)

so that the two time correlation function may be expressed as
(A(t + ?)B(t)) = Trs{AeL? B? (t)} . (6.165)



6.5.1 Quantum Regression Theorem


In cases where the master equation gives linear equations for the mean, we can de- velop a quantum regression theorem, similar to that for ordinary Markov processes. This result was first derived by Lax [4].
Suppose for a certain set of operators Yi, the master equation can be shown to yield, for any initial ?




Then we assert that
?

?
? t (Yi(t)) = ?Gij (t)(Yj (t)) . (6.166)



For

? t (Yi(t + ?)Yl (t)) = ?Gij (?)(Yj (t + ?)Yl (t)) . (6.167)

(Yi(t + ?)Yl (t)) = Trs {YieL? Yl ? (t)} (6.168)

the right-hand side is an average of Yi at time t + ?, with the choice of initial density matrix
?init = Yl ? (t) . (6.169)

Since by hypothesis, any initial ? is permitted and the equation is linear, we may generate any initial condition whatsoever. Hence, choosing ?init as defined in (6.169) the hypothesis (6.166) yields the result (6.167) which is the quantum regression theorem.



Application to Systems with a P Representation


For systems where a P representation exists the following results for normally or- dered time correlation functions may be proved
G(1)(t, ?)= (a†(t + ?)a(t)) = (??(t + ?)?(t)) , (6.170)
G(2)(t, ?)= (a†(t)a†(t + ?)a(t + ?)a(t)),
= (|?(t + ?)|2|?(t)|2 ) . (6.171)


In these cases the measured correlation functions correspond to the same correla- tion function for the variables in the P representation. For non-normally ordered correlation functions the result is not as simple.



Stochastic Unravellings


The master equation describes the dynamics of a subsystem by averaging over (tracing out) the properties of the larger “bath” to which it is coupled. Solving the master equation typically results in a mixed state. Any mixed state admits infinitely many decompositions into convex combinations of (non-orthogonal) pure states. In a stochastic unravelling of a master equation we represent the solution at any time as a convex combination of pure states each evolving under a stochastic Schro¨ dinger equation such that if we average over the noise we obtain the solution to the orig- inal master equation. This approach leads to a powerful numerical simulation tool as much less memory is required to store a pure quantum state at each time step. In Chap. 15, we give an alternative interpretation of an unravelling in terms of con- ditional states conditioned on a continuously recorded sequence of measurement results.
Consider a simple harmonic oscillator coupled to a zero temperature heat bath. The dynamics, in the interaction picture, is given by the master equation, (6.37) with N = M = 0. Solving this equation over a small time interval dt we can write

?

? (t + dt)= .

(t) ?

(a†a? (t)+ ? (t)a†a)dt. + ?a? (t)a†dt (6.172)
2


We can think of this as describing photons leaking from a single mode cavity at Poisson distributed times. Suppose there were exactly n photons in the cavity at
time t so that ? (t)= |n)(n|. Then (6.172) would become
? (t + dt)= (1 ? ?ndt)|n)(n| + ?ndt|n ? 1)(n ? 1| (6.173) We can think of this as follows. In a small increment of time dt, two events are
possible: either a single photon is lost or no photon is lost. If a photon is lost, the state of the field has one less photon so that it changes from |n) to |n ? 1) and
this event will occur with probability ?ndt. This form results from the last term of (6.172). On the other hand, if no photon is lost the state is unchanged, and this
will occur with probability 1 ? ?ndt, which arises from the first term in (6.172).
Thus (6.172) describes a statistical mixture of the two events that can occur in a
small time step dt: the first term in square brackets describes the change in the state of the cavity field given that no photon is lost in time interval dt, while the second term describes what happens to the state of the field if one photon is lost in a time interval dt.
If this interpretation is correct it suggests an answer to conditional questions such as: if no photon is lost from time t to t + dt, what is the conditional state of the field?


In this case we have no contribution from the last term in (6.172) so the conditional state is the solution to


? (t + dt)c =

[? (t)c ? ? (a†a? (t)c + ? (t)ca†a)dt]
Tr[? (t)c ? ? (a†a? (t)c + ? (t)ca†a)dt]

. 1 †

† † .

? ?c(t) ? ?dt

.a a?c(t)+ ?c(t)a a. ? (a a)c(t)?c(t)


to linear order in dt, where the subscript c is to remind us that we are dealing with a particular conditional state conditioned on a rather special history of null events and
(a†a)c(t) is the conditional average of the photon number in the state ?c(t).
We can now introduce a classical stochastic process, a conditional Poisson pro-
cess, dN(t) which is the number of photons lost in time dt. Clearly

dN(t)2 = dN(t) (6.174)
E [dN(t)] = ?(a†a)c(t) (6.175)

where E is an average over the classical stochastic variable. In terms of dN(t) we can now define a stochastic master equation
d?c(t)= dN(t)G [a]?c(t) ? ?dtH [a†a]?c(t) (6.176) where we have defined two new super-operators (that map density operators to den-
sity operators),

A? A†
G [A]? = Tr[A? c†] ? ? (6.177)
H [A]? = A? + ? A† ? Tr[A? + ? A†] (6.178)

for any operator A. Note that if we take the classical ensemble average over the noise process dN(t) we recover the original unconditional master equation in (6.172). The solution to (6.176) is the conditional state at time t conditioned on an entire fine- grained history of jump events (that is to say, the total number of jumps and the time of each jump event). Denote such a history as the sequence of jump times on the
interval [0,t) as h[t]= {t1,t2 ,...,tm}. The unconditional state is a sum over all such
histories
? (t)= ?Pr(h[t])?c(h[t]) (6.179)
h[t]

where we have explicitly indicated that the conditional state ?c is conditioned on the history of jump events, h[t] in the time interval of interest and Pr(h[t]) is the probability for each history. We have unravelled the solution to the master equation in terms of conditional stochastic events. For a point process as considered here the sum over histories has an explicit form in terms of time ordered integrals [5]

? ¸ t
? (t)= ?
m=0 0



dtm


tm t1
¸ ¸
•••
0 0



S (t ? tm)JS (tm ? tm?1) ... JS (t1)? (0) (6.180)


where the super operators are defined by

? ?
S (t)? = e? 2 ta† a? e? 2 ta† a (6.181)
J = ?a? a† (6.182)

Clearly the probability of a specific jump history is given by
Pr(h[t]) = Tr[S (t ? tm)JS (tm ? tm?1) ... JS (t1)? (0)] (6.183)

The form of (6.180) indicates that if we start in a pure state, and have access to the entire history of photon loss events h[t], the conditional state ?c(h[t]) must still be a pure state. This implies that we can write a stochastic Schro¨dinger for the damped harmonic oscillator:


d|?c(t)) =

.
dNc(t)

. a
(a†a)c(t)

.
? 1 + ?dt

. (a†a)c(t)
2

a†a .
? 2

.
? iHdt


|?c(t))

(6.184)

where we have now included the possibility of a hamiltonian part to the dynam- ics. We can show the equivalence between this equation and the stochastic master equation by considering the Ito-like expansion
d(|?c(t))(?c(t)|)= (d|?c(t)))(?c(t)| + |?c(t))(d(?c(t)|)+ (d(|?c(t))|)(d(?c(t)|)
(6.185)

and retaining all terms to first order in dt, noting that dN2 = dN.
A point process with a large rate parameter ? can be well approximated on a time scale long compared to ??1 by a white noise process. This suggests that it must be
possible to unravel the master equation in terms of white noise processes as well as the point process dN(t). In Chap. 15 we will see that such master equations give the conditional dynamics conditioned on homodyne and heterodyne measurements on the field leaving the cavity. Here we simply quote the result and show that averaging over the classical noise returns us to the unconditioned master equation.
In the case of the real valued Weiner process dW (t) we can write the homodyne
stochastic master equation for a dampled simple as

d?c(t)= ?i[H, ?c(t)]dt + D [a]?c(t)dt + dW (t)H [a]?c(t) (6.186)


where



D [A]? = A?A† + 1 (A†A? + ?A†A) (6.187)
2

and H [a] is given in (6.178)


In terms of a complex valued white noise process, dW (t) = dW1(t)+ idW2(t) where dWi(t) are independent Winer processes, we can write the heterodyne stochas- tic master equation
1
d?c(t) = ?i[H, ?c(t)]dt + D [a]?c(t)dt + ?2 (dW1(t)H [a]?c(t)
+dW2(t)H [?ia])?c(t) (6.188)

There is a connection between the quantum jump process and the two-time cor- relation function discussed in Sect. 6.5. We first define a new stochastic process, the rate or the current, as

dN i(t)=
dt

(6.189)

This is a rather singular stochastic process, comprised of a series of delta functions concentrated at the actual jump times. In physical terms this is intended to model the output of an ideal photon counting detector, with infinite response bandwidth, that detects every photon that is lost from the cavity. Define the classical current two-time correlation function

G(?,t)= E (i(t + ?)i(t)) (6.190)

Given the nature of a Poisson jump process, dN(t) can only take the values 0 or 1, so it is easy to see that we can write the two-time correlation function in terms of the conditional probability to get dN(t + ?)= 1 given a jump at time t,
G(?,t)dt2 = Pr(dN(t + ?)= 1|dN(t)= 1) (6.191)

This conditional probability is given by
Pr(dN(t + ?)= 1|dN(t)= 1)= ?2Tr[a†aeL ? a? (t)a†]dt2 (6.192) where eL t is the formal solution to the unconditional master equation evolution
written in terms of the abstract generator L as ?? = L ? . The two time correlation function is then given by

G(?,t)= ?2Tr[a†aeL ? a? (t)a†] (6.193) Note that the so-called regression theorem follows directly from the definition of G,
dG(?,t) = L G(?,t) (6.194)
d?

We usually deal with driven damped harmonic oscillators for which the system settles into a steady state, emitting photons according to the conditional Poisson process derived from the steady state solution ?? = limt?? ? (t), so we define the
stationary two-time correlation function for the current as


G(?)= ?2Tr a†


aeL ?

a??a†.


(6.195)


6.7.1 Simulating Quantum Trajectories


A mixed state for system with a Hilbert space of dimension N requires that we specify N2 complex matrix elements. On the other hand a pure state requires that we specify only N complex numbers. For this reason numerically solving the master equation is more computationally difficult than solving the Scho¨dinger equation. We can use the unravelling of a master equation in terms of a stochastic Scho¨dinger equation to make the numerical solution of master equations more tractable. In this numerical setting, the method of quantum trajectories was independently developed as the Monte-Carlo wavefunction method [6]. We will illustrate the method using the jump process.
Suppose the state at time t is |?(t)). Then in a time interval ? t, sufficiently short
compared to ??1, the system will evolve to the (unnormalised) state conditioned on
no-jump having occurred,

|?? (t + ? t)) = e?iH? t??a a? t/2|?(t)) (6.196) To compute this we implement a routine to solve the Schro¨dinger equation with the
effective non-hermitian Hamiltonian
? †
K = H ? i 2 a a (6.197)

The norm of this state is the probability that no-jump has occurred in the time inter- val ?t,






where it is easy to see that

p0 = (?? (t + ? t)|?? (t + ? t)) (6.198)
= 1 ? p (6.199)


p = ?? t(?(t)|a†a|?(t)) (6.200)


which we understand to be the probability that a jump takes place in this time inter- val. We need to ensure that p << 1.
Let us now chose a random number r uniformly distributed on the unit interval. At the end of the time interval, we compare p and r. If p < r (usually the case) we
normalise the state
?? (t + ? t ))

|?(t + ? t)) = | ??
0

(6.201)

and continue the non-hermitian evolution for a further time step. If however p > r, we implement a quantum jump via


|?? (t + ? t))? |?(t + ? t)) =

??a|?? (t + ? t)) (6.202)
p/? t


Based on our previous discussions we see that p/? t = ?(?(t)|a†a|?(t)), and the jump operation is as described by the first term in (6.184). As the simulation pro-
ceeds we accumulate the record of times at which particular jump events occur. That is to say, we have access to a sample fine grained history of the jump process, h[t]. However we are primarily interested in solving the master equation. We thus run K trials up to time t, starting from an identical initial state each time, and then form the mixed state

k
?¯ (t)= K?1 ? |?k(t)) (6.203)
k=1
as a uniform average over the K trials. Strictly speaking the probability of each of the K trials is not uniform, however one can show that for K sufficiently large
?¯ (t) ? ? (t) with a error that scales as K?1/2. In M?lmer et al. [6] more general cases
are discussed including how to simulate non zero temperature master equations or
master equations with multiple jump processes.



Exercises


The photon number distribution for a laser may be shown to obey the master equation


d
P(n)=
dt


An
1 + n/ns


P(n ? 1) ?


A(n + 1)
1 + (n + 1)/ns


P(n) ? ?nP(n)+ ?(n + 1)P(n + 1),

where A is related to the gain, ns is the saturation photon number and ? is the cavity loss rate.
Use detailed balance to show that the steady state solution is
. Ans .n
?
Pss(n)= N (n + n )!

where N is a normalisation constant.
The interaction picture master equation for a damped harmonic oscillator driven by a resonant linear force is


d? = i?[a + a†, ? ]+ ? (2a? a†


a†a?


? a†a) .

dt 2 ? ?
Show that the steady state solution is the coherent state |2i?/?).
A model for phase diffusion of a simple harmonic oscillator is provided by the
master equation

d? =
dt

??[a†a, [a†a, ? ]] .

Show that the Q function obeys the Fokker–Planck equation.

Further Reading 125

? Q ? . ?
=

?Q +

? ?? Q + 2 ?

|?|2Q ?

? 2
?2Q ?

? ??2Q. .

? t 2 ??

???

?????

??2

???2


Thus show that while the mean amplitude decays the energy remains constant. Using intensity and phase variables ? = I1/2ei? show that the model is simply a diffusion process for the phase.
Show that in terms of the quadrature operators X1 = a + a†, X2 = ?i(a ? a†),
the master equation (6.37) may be written


d? =i ? [X ,



X , ? ]


i ? [X ,



X , ? ]

dt 8 2 { 1
? 2r

} ? 8

1 { 2 }
? 2r

? 8 e

[X1, [X1, ? ]] ? 8 e?

[X2, [X2, ? ]]

where {,} is an anticommutator and we have taken N = sin h2r, M = sin hr cos h r for an ideal squeezed bath. Show that the first and second terms de- scribe damping in X1 and X2 respectively, while the third and fourth terms describe diffusion in X2 and X1, respectively.
Show that the homodyne conditional master equation for a driven simple har-
monic oscillator, damped into a zero temperature heat bath, has the same pure steady state as the unconditional case (Exercise 6.2).
If a cavity mode starts in a state for which the P-representation is Gaussian, show that under the conditional dynamics of the homodyne master in (6.15), the state continues to have a Gaussian P representation



References


1. F. Haake: (private communication)
2. W.H. Louisell: Quantum Statistical Properties of Radiation (Wiley, New York 1973)
3. M.C. Wang, G.E. Uhlenbeck: Rev. Mod. Phys. 17, 323 (1945)
4. M. Lax: In Brandeis University Summer Institute Lectures, ed. by M. Chretien, E.P. Gross,
S. Deser (Gordon and Breach, New York 1966) Vol. 2
5. M.D. Srinivas and E.B. Davies: J. Mod. Opt. 28, 981 (1981)
6. K. M?lmer, Y. Castin and J. Dalibard, J. Opt. Soc. Am. B, 10, 524 (1993)



Further Reading


Agarwal, G.S.: Quantum Statistical Theories of Spontaneous Emission and Their Relation to Other Approaches, Springer Tracts Mod. Phy., Vol. 70 (Springer, 1974)
Carmichael, H.J.: Statistical Methods in Quantum Optics 1, Master Equations and Fokker Planck Equations (Springer, 1999)