Open Quantum Systems with Applications to Precision
Measurements
byDavid Tieri
B.A., University of Colorado - Boulder, 2009
M.S., University of Colorado - Boulder, 2013

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillmentof the requirements for the degree ofDoctor of Philosophy
Department of Physics
2015This thesis entitled:Open Quantum Systems with Applications to Precision Measurements
written by David Tieri
has been approved for the Department of Physics

Murray Holland

John Bohn

Date

The final copy of this thesis has been examined by the signatories, and we find that both thecontent and the form meet acceptable presentation standards of scholarly work in the above
mentioned discipline.

iii
Tieri, David (Ph.D., Physics)
Open Quantum Systems with Applications to Precision Measurements
Thesis directed by Prof. Murray Holland

A spectrally pure coherent light source is an important component in precision measurement
applications, such as an atomic clock. The more spectrally pure the coherent light source, or the
narrower the linewidth of its power spectrum, the better for atomic clock experiments. A coherent
light light source, such as a laser, is intrinsically an open quantum system, meaning that it gains
and loses energy from an external environment.
The aim of this thesis is to study various open quantum systems in an attempt to discover ascheme in which an extremely spectrally pure coherent light source might be realized. Therefore,
this thesis begins by introducing the two main approaches to treating open quantum systems, thequantum master equation approach, and the quantum Langevin equation approach. In addition toderiving these from first principles, many of the solution methods to these approaches are given andthen demonstrated using computer simulations. These include the quantum jump algorithm, thequantum state diffusion algorithm, the cumulant expansion method, and the method of c-number
Langevin equations.
Using these methods, the theory of the crossover between lasing and steady state superradiance is presented. It is shown that lasing and steady state superradiance might be demonstrated in
the same physical system, but in different parameter regimes. The parameter space between these
two extreme limits is explored, and the benefits and drawbacks of operating a system at a given
set of parameters, i.e. to achieve the most spectrally pure light source, are discussed.
We also consider the phase stability of a laser that is locked to a cavity QED system comprisedof atoms with an ultra-narrow optical transition. Although the atomic motion introduces Doppler
broadening, the standing wave nature of the cavity causes saturated absorption, which can be used
to achieve an extremely high degree of phase stabilization. The inhomogeneity introduced by finite

iv
atomic velocities can also cause optical bistability to disappear, resulting in no regions of dynamic
instability that would otherwise restrict operational parameters in the experiment to be tuned
outside of the optimum region where the minimum linewidth occurs.

Dedication

This thesis is dedicated to my Grandma, Mom and Dad, Michael, Granny and Bob, Theo,UJ and Aunt Laura, Auntie Lori, Aunt Cathy and both of my Grandfathers.

vi

Acknowledgements

I would like to thank Jinx Cooper for teaching me to never accept that something is true justbecause someone else says it is. I know it took a lot of convincing. Thank you Murray Holland
for providing me with such amazing opportunities throughout my physics career. I would like toalso thank Minghui Xu for being a good friend; I can honestly say I would not have gotten a PhD
without you. To Dominic Meiser and Brandon Peden, thanks for being mentors to me, and forspending so much of your time teaching me stuffs. Thanks Jose D’Ianco for being an awesome
person, Ron Pepino, Rajiv Bhat, Meret Krämer for being so welcoming when I first started in
the group, Andrew Sykes for always being there to chat over a beer, Brandon Ruzic for beinga good friend through it all, especially all semester last semester. Ryan Wilson, Victor Colussi,
Jon Pfeiffer, you guys are awesome too. Thanks Jon for switching rooms with me when we lived
together and I didn’t like my room.
I would also like to thank Mary Nelson for teaching me Calculus in her office at 7am every
morning for half a year, Ben Safdi for getting me started in physics, and staying up until 3am
helping me do my homework, and then starting on his. Thank you to Noel Clark for telling me
to go home and have a beer after a hard exam, to Mark Ablowitz for taking the time to tell me”Its ok, you’re doing a good job, buddy”, when I was totally stressed out in his office ten minutes
before a homework trying to get the answers out of him. Thank you Will Forsythe for inspiring meeach week, Kim Browne, Jeff Walters, Marika, Jimie and Scott, Ryan at Robb’s Music, Andrew
and Amanda, and everyone from Pekoe. Thanks to Bihui and Athreya, and to Jan Thomsen andhis group for a great collaboration. And thanks to everyone else I forgot to thank.

viiContents

Chapter
1 Introduction

1

1.1

Background . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

2 Preliminaries

7

2.1

Quantization of the Electromagnetic Field . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Jaynes Cummings Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.3

7

2.2.1

Jaynes Cummings Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.2

The Rotating Wave Approximation . . . . . . . . . . . . . . . . . . . . . . . . 14

2.2.3

Rabi Flopping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

Wiener-Khinchin Theorem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

3 Master Equation

19

3.1

Density Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

3.2

The Master Equation for a Damped Harmonic Oscillator . . . . . . . . . . . . . . . . 22
3.2.1

Derivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2.2

Analytic Solutions for Photon Number . . . . . . . . . . . . . . . . . . . . . . 28

3.2.3

The Quantum Regression Theorem . . . . . . . . . . . . . . . . . . . . . . . . 29

3.3

The Master Equation for a Damped Two Level Atom

. . . . . . . . . . . . . . . . . 31

3.4

Quantum Jumps . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

viii

3.5

3.4.1

The Dyson Expansion of the Master Equation

3.4.2

Definition of Quantum Trajectories . . . . . . . . . . . . . . . . . . . . . . . . 35

3.4.3

Unraveling the Master Equation . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.4.4

The Quantum Jumps Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.4.5

Computer Simulation of Quantum Jumps . . . . . . . . . . . . . . . . . . . . 40

Quantum State Diffusion. . . . . . . . . . . . . . . . . 33

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

3.5.1

Master Equation of a Homodyne System . . . . . . . . . . . . . . . . . . . . . 42

3.5.2

Computer Simulation of the Quantum State Diffusion Equation . . . . . . . . 44

3.6

The Master Equation for a General Cavity QED System . . . . . . . . . . . . . . . . 46

3.7

Cumulant Expansion to Second Order . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.8

3.7.1

Closed Set of One-Time Expectation Values . . . . . . . . . . . . . . . . . . . 48

3.7.2

Closed Set of Two-Time Expectation Values . . . . . . . . . . . . . . . . . . . 52

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56

4 Langevin Equations

57

4.1

Derivation of the Quantum Langevin Equations . . . . . . . . . . . . . . . . . . . . . 58

4.2

Adiabatic Elimination and the Superradiance Master Equation . . . . . . . . . . . . 61

4.3

Fokker-Planck Equations for the P Function . . . . . . . . . . . . . . . . . . . . . . . 63

4.4

Equivalence of Langevin and Fokker-Planck Equations . . . . . . . . . . . . . . . . . 65

4.5

C-number Langevin Equations For Coherent State Amplitudes . . . . . . . . . . . . 68

4.6

A More Practical Derivation of Quantum Langevin Equations . . . . . . . . . . . . . 69

4.7

4.6.1

The Generalized Einstein Relations . . . . . . . . . . . . . . . . . . . . . . . . 69

4.6.2

Obtaining the Drift Terms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4.6.3

Obtaining the Diffusion Coefficients . . . . . . . . . . . . . . . . . . . . . . . 75

A More Practical Derivation of the C-number Langevin Equations . . . . . . . . . . 77
4.7.1

Solving for the C-Number Diffusion Coefficients . . . . . . . . . . . . . . . . . 77

4.7.2

Numerical Simulation of C-number Langevin Equations . . . . . . . . . . . . 80

ix
4.8

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

5 The Crossover from Lasing to Steady State Superradiance

85

5.1

Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

5.2

Solution Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
5.2.1

Langevin Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

5.2.2

2nd Order Cumulant Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

5.3

Characterization of the Crossover . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

5.4

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.5

5.4.1

Simulations with small atom numbers . . . . . . . . . . . . . . . . . . . . . . 98

5.4.2

Simulations with large atom numbers . . . . . . . . . . . . . . . . . . . . . . 102

5.4.3

Effect of Cavity and Atomic Level Instabilities on the Spectrum . . . . . . . 104

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

6 Laser Stabilization using Saturated Absorption in a Cavity QED System

108

6.1

Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

6.2

Discussion of Steady State Solutions . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

6.3

Shot Noise Limited Laser Stabilization . . . . . . . . . . . . . . . . . . . . . . . . . . 117

6.4

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 124

7 Conclusion

125

Bibliography

128

x

FiguresFigure3.1

The quantum jumps algorithm, used to find â† (t)â(t) for a damped harmonic
oscillator with n̄ = 1, and α(0) = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

3.2

The quantum state diffusion equation is solved for random trajectories which areused to find â† (t)â(t) for a damped harmonic oscillator with n̄ = 1, and α(0) = 2.

4.1

45

The c-number Langevin equations are used to find â† (t)â(t) for a damped harmonic
oscillator with n̄ = 1, and α(0) = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.2

The c-number Langevin equations are used to find g (1) (τ ) for a damped harmonic
oscillator with n̄ = 1, and α(0) = 2. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

4.3

The c-number Langevin equations for symmetric ordering are used to find the inverDEsion Jˆz (t) for resonance fluorescence. . . . . . . . . . . . . . . . . . . . . . . . . . 81

4.4

The c-number Langevin equations for symmetric ordering are used to find the Fouriertransform of g (1) (τ ) for resonance fluorescence. . . . . . . . . . . . . . . . . . . . . . 82

5.1

Comparison of the different solution methods in the superradiance, crossover, andlasing regions for N = 40 and C = 1.

5.2

. . . . . . . . . . . . . . . . . . . . . . . . . . 99

Solutions using the various methods in the superradiance, crossover, and lasing regions, with and without dephasing, for N=10000 and C = 0.1 . . . . . . . . . . . . . 101

5.3

Comparison of the linewidth and intracavity intensity for a system operated in the
crossover, lasing, and far lasing regions. . . . . . . . . . . . . . . . . . . . . . . . . . 103

xi
5.4

Instability in the atom-cavity system frequency ω with respect to the cavity frequency
ωc and atomic frequency ωc as a function of crossover parameter ξ . . . . . . . . . . 105

6.1

Schematic of the cavity QED experiment . . . . . . . . . . . . . . . . . . . . . . . . . 110

6.2

Development of the extra absorption and dispersion features in the transmission andphase shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

6.3

Measured phase shift of the cavity transmitted field when scanned across the atomic
resonance, for several input probe laser powers . . . . . . . . . . . . . . . . . . . . . 118

6.4

Measured phase shift of the cavity transmitted field when scanned across the atomic
resonance for several numbers of atoms . . . . . . . . . . . . . . . . . . . . . . . . . . 119

6.5

Linewidth at the optimum input intensity, as a function of included Doppleron orders121

6.6

Intracavity intensity and Slope at resonance as a function of input intensity, andStabilization linewidth as a function of N C0 , for several temperatures. . . . . . . . . 123

Chapter 1

Introduction

1.1

BackgroundOur understanding of the interaction between light and matter is important, both fundamen-tally, and because it allows for the development of applications that are integral to the operationof our society. For instance, it was this understanding that allowed for the invention of the laser in
1960, which was then called ”a solution looking for a problem” [82]. Since then, lasers have become
ubiquitous in our society, with a vast scope of applicability, ranging from medical applications, suchas cancer diagnosis and treatment [25, 26, 36], eye and general surgeries [21], and medical imaging[33], to industrial applications, such as cutting and welding [5]. Lasers are also used in nuclear
fusion research to deliver a tremendous amount of energy to an extremely small fuel cell, in an
attempt to achieve conditions where nuclear fusion [53, 54] will happen. If achieved, this internal
confinement nuclear fusion would provide an inexaustible source of clean energy.
In fundamental biology research, optical tweezers, which rely on a highly focused laser andcan hold microscopic particles in place [3], have opened the doors to new experiments, including thestudy of motor proteins [78], which are the mechanism behind transport phenomenon in neurons.
The information obtained from optical tweezer experiments could allow for the development of a
new generation of drugs and treatments to fight disease related to motor proteins [61].
Our understanding of the interaction between light and matter also has allowed for thedevelopment of the atomic clock. The idea was proposed in 1945 by Rabi [4, 72], and by 1955, an
atomic clock based on a microwave transition in Cesium was built [29, 28, 30], which surpassed in2accuracy and precision, the best quartz based clocks of the day [20, 47]. This led in 1967 to theredefinition of the second, in terms of the frequency between two atomic levels in Cesium [72].
Since then, atomic clocks have increased in accuracy and precision by about an order ofmagnitude every decade [46], and a vast number of new technologies, including GPS, have followed.
GPS is a space-based system that uses multiple synchronized atomic clocks to provide time andlocation information anywhere on Earth, and it has revolutionized navigation in all aspects of oursociety.
Recently, a new generation of atomic clocks based on optical transitions has been created[83, 41, 44, 8]. These new optical atomic clocks were made possible by the new techniques oftrapping and cooling of atoms with lasers, and by the invention of the optical frequency comb. An
optical frequency comb is a pulsed laser that acts like reduction gears in the frequency domain,
allowing optical frequencies to be geared down to measurable microwave frequencies. Optical atomic
clocks have now surpassed the cesium standard in precision and accuracy by orders of magnitude,
and still are becoming even more precise and accurate [8]. In addition to ”being another solution
looking for a problem”, more precise and accurate clocks could be used to detect potential slow
changes to the fundamental constants of nature [7], or to confirm the existence of gravitational
waves predicted in Einstein’s general theory of relativity [6].
Currently, the precision and accuracy of these optical clocks are limited by the linewidth ofthe laser used to interrogate the atoms through the Dick effect [18, 71]. Here, the cyclical operationof these clocks causes an aliasing effect that translates the interrogation laser’s phase noise into the
feedback control loop that keeps the laser on the atomic resonance. The noise in the interrogation
laser is a bottleneck which limits the achievement of better precision and accuracy, and a narrower
linewidth interrogation laser is needed to increase the precision and accuracy of these clocks.
Much effort has already been invested in making a narrow linewidth interrogation laser. The
current state-of-the-art is to send laser light into an ultra-stable reference cavity to reduce thelaser’s linewidth [42, 23, 91, 51, 50]. This method has made for orders of magnitude reduction
in interrogation laser linewidth, but it is limited in further reducing this linewidth by the cavity

3length fluctuations due to the thermal vibrations of the cavity mirrors. Even when using thesecavity stabilized lasers, optical clocks still suffer from the Dick effect, so that further reducing
the linewidth of the interrogation laser would directly allow more accurate and precise clocks.
Reduction of a laser’s linewidth by this method, past the current state of the art, is very difficult.
Because of this difficulty, alternative methods of producing a narrow linewidth laser havebeen considered. It was proposed, that in addition to a laser, coherent light can be produced in
an atom-cavity system that is operated in an exotic parameter regime. In this regime, the cavity
decay rate is much larger than the collective atomic decay rate, which is smaller than the decayrate of independent atoms [66, 65, 64]. In this so called ”bad cavity” or superradiance regime, thelinewidth can actually be smaller than the atomic transition itself, due to collective effects. In fact,
this light could be used directly as the frequency measure for a clock, circumventing the need for
an interrogation laser all together [13]. This type of coherent light source has been explored onlyin the extreme superradiance limit, and the crossover between superradiance and the more familiar
lasing limit has only recently been explored [81].
Another proposed method of producing a narrow linewidth laser is to place atoms with anextremely narrow transition inside an optical cavity, and to stabilize a laser to an atomic transition,
instead of to an empty cavity resonance [63]. Such a system, which again lies in the bad cavity, or
superradiance parameter regime, is insensitive to thermal cavity length fluctuations, and therefore
circumnavigates the main bottleneck of further reducing the linewidth of empty reference cavity
systems. An experimental effort to realize such a system has been made [86], but in contrast to theoriginal proposal [63], the atoms had a non-negligible velocity, making Doppler effects important.
Still, it was shown theoretically how these systems could be optimized to allow for linewidths
comparable, or even narrower, than the state of the art reference cavity systems [80].
In this thesis, we will explore in detail the topic of the crossover theory from lasing to superradiance, and demonstrate the potential for a crossover system to be used in precision measurement
applications [81]. We will also explore the theory of the stabilization of a laser to narrow linewidth
atoms inside a cavity, including the effects due to atomic motion [86, 80, 15].

4
Before these topics are explored, the background material necessary for their understanding
will be discussed. Specifically, these aforementioned systems are inherently open, meaning that they
are able to lose and gain energy. Therefore, the equations that describe open quantum systems
will be derived. Exact solution methods to these equations will be discussed. However, thesemethods are only able to describe small systems, with ∼ 10 atoms. A solution method does exist
which exploits an underlying symmetry in the derived equations, and therefore drastically reducescomplexity, allowing for larger systems to be analyzed [90]. Still, this method is limited to systems
having ∼ 100 atoms, and therefore cannot be used to directly describe many realistic experiments,
which can have ∼ 105 or more atoms. To describe these larger experimental systems, approximate
solution methods must be used. Therefore, we discuss several approximate solution methods. In
addition to their ability to describe the specific systems previously mentioned, these approximate
solution methods are able to capture the essential physical phenomenon in a wide class of devices,
including the synchronization between two different ensembles of atoms [89].

1.2

Overview
In Chapter 2, the quantization of the electromagnetic field will be discussed. This will provenecessary in order to correctly describe the Hamiltonian dynamics of our open quantum systems
when few photons are present. Also, it will be necessary in order to describe spontaneous irreversible
processes, such spontaneous emission. The power spectrum, and its relation to two-time correlation
functions will also be discussed. This will allow the spectral properties of a system, such as the
linewidth of the power spectrum, to be calculated.
In Chapter 3, the equation for open quantum systems in the Schrödinger picture of quantum
mechanics, known as the master equation, will be derived. Exact solution methods to the masterequation, such as the quantum jump and quantum state diffusion methods, will be discussed.
In addition, computer simulations of these methods are demonstrated. The cumulant expansion
method, which yields an approximate set of coupled equations for expectation values of systemobservables, is then discussed.

5
Chapter 4 will begin by deriving the quantum Langevin equations, which describe open
quantum systems in the Heisenberg picture of quantum mechanics. Although these equations canprove useful for some analytical work, for computer simulations, they are very difficult to work with.
Therefore the Fokker-Planck equation for a quazi-probability distribution describing the system,
and its equivalent set of c-number (complex number) Langevin equations will be discussed. Thesec-number Langevin equations can be simulated on a computer with relative ease. Once this proofof-principle procedure for deriving the c-number Langevin equations is shown, a more practical
way to derive c-number Langevin equations is given. This new procedure becomes more and moreadvantageous as system complexity increases. Computer simulations of the c-number Langevinequations for several systems will be demonstrated.
Chapter 5 will use several of the methods of the previous chapters to present the theory of thecrossover between lasing and steady state superradiance. Although lasing and steady state superradiance seem to be disjoint phenomena, it will be shown that they should actually be thought of as
the extreme limits of a single and more general phenomenon. In a laser, the phase information thatallows for coherence is stored in the light field. In steady state superradiance, the phase informationthat allows for coherence is stored in the atomic ensemble. In the in-between parameter region, or
crossover region, phase information that allows for coherence can be stored in the light field andthe atomic ensemble simultaneously. The properties of systems that lie in the superradiance, lasing
and crossover parameter regions will be compared.
The rate at which energy is re-supplied to a typical laser allows the system to operate justabove the lasing threshold. It will be seen that a crossover system operating at this same rate ofre-supplying energy will be much farther above threshold, which can allow the linewidth of thepower spectrum to be orders of magnitude smaller than the laser. Such a system also exhibits
a much larger intracavity intensity than a superradiant system, so that the crossover system ismore experimentally accessible. The linewidth and intracavity intensity will then be shown to beinsensitive to atomic dephasing in all parameter regions, so long as the rate of atomic dephasing issmaller than the rate of energy being resupplied to the system.

6
It will also be demonstrated that a system in the superradiance region is insensitive tofluctuations in cavity frequency, while a system in the lasing region is insensitive to fluctuations in
the atomic frequency. A crossover system is less sensitive to fluctuations in cavity frequency thana laser system, and less sensitive to fluctuations in atomic frequency than a superradiant system.
This is important, since fluctuations in cavity length are the cause of the current limit on thefurther reduction of the linewidth in today’s best ultrastable lasers.
In Chapter 6, we will consider the phase stability of a local oscillator (or laser) locked to acavity QED (Quantum Electrodynamics) system comprised of atoms with an ultra-narrow opticaltransition. The atoms will be modeled as being cooled to millikelvin temperatures and then released
into the optical cavity. Although the atomic motion introduces Doppler broadening, the standing
wave nature of the cavity causes saturated absorption features to appear, which are much narrower
than the Doppler width. These features can be used to achieve an extremely high degree of phasestabilization, competitive with the current state-of-the-art. A comparison between the developed
theory, and the results of an experiment on the system that we consider will be given. Finally, it will
be seen that the inhomogeneity introduced by finite atomic velocities can cause optical bistability todisappear, resulting in no regions of dynamic instability and thus enabling a new regime accessible
to experiments where optimum stabilization may be achieved.
Finally, chapter 7 will give a conclusion of the presented material.

Chapter 2

Preliminaries

The equations that describe open quantum systems are comprised of both coherent Hamiltonian dynamics, and incoherent dynamics. In the systems that we will consider, the Hamiltonianparts will describe atoms interacting reversibly with an electromagnetic field, meaning that energy
is transferred between the two, but is not lost to an external environment. If there are only a few
photons present in the field, a quantized electromagnetic field will be necessary for an accurate
description. A quantized electromagnetic field will also be required to describe the incoherent partsof the systems, in which energy is lost to an external environment. These parts are derived bycoupling the system to a continuum of quantized electromagnetic field modes in vacuum.
In addition, the spectral properties of the light produced by a system will be of interest to us.
This aspect concerns the Wiener-Khinchin theorem. This theorem relates the power spectrum toa two-time correlation function of field variables. In many cases the two-time correlation function
can be obtained with relative ease.

2.1

Quantization of the Electromagnetic Field
We begin by formulating the quantized theory of the electromagnetic radiation and demon-

strating the discrete photon excitations that this implies. Maxwell’s equations in vacuum and8without sources are,∇ · E(r, t) = 0(2.1)∇ · B(r, t) = 0(2.2)∂B(r, t)∂t
∂E(r, t)∇ × B(r, t) = µ0 0.
∂t∇ × E(r, t) = −

(2.3)
(2.4)

Eq. 2.2 allows B(r, t) to be written as the curl of a vector potential A(r, t),B(r, t) = ∇ × A(r, t).

(2.5)

In the Coulomb gauge, where ∇ · A(r, t) = 0, E(r, t) can be written as,
E(r, t) = −dA(r, t)
.
dt

(2.6)

Substituting Eqns. 2.5 and 2.6 in Eq. 2.4 yields,∇2 A(r, t) = µ0 0∂ 2 A(r, t)
.
∂t2

(2.7)

The general solution to Eq. 2.7 inside a cubical volume of side length L, with boundary conditions
corresponding to standing wave modes is,A(r, t) =X
Ak uk (r)e−ıωk t + A∗k u∗k (r)eıωk t ,

(2.8)

k

where,k=2πn , n = (nx , ny , nz ) , ni = 0, ±1, ±2, ± . . . ,
L

(2.9)

and,
ωk = c |k| .

(2.10)

 3/2
2
sin(k · r)ek ,L(2.11)The mode functions uk (r) are given by,
uk (r) =9where ek is a unit polarization vector orthogonal to k. These mode functions satisfy,
ω2∇2 uk (r) = − 2k uk (r),c(2.12)√
where c = 1/ µ0 0 is the speed of light in a vacuum. Also, they form a complete orthonormal
basis, i.e.,Zu∗k (r)uk0 (r)dr = δkk0 ,(2.13)

V

where the integral is over the cubical volume V , and δkk0 is the Kronecker delta. Note that for
simplicity we have neglected the discrete polarization index, but that can be accounted for byletting the index k → kλ and also summing over λ, the discrete polarization index.
The energy H stored in the electromagnetic field is
1
H=
2

Z




1
dV 0 (E(r, t) · E(r, t)) +(B(r, t) · B(r, t)) .
µ0(2.14)

From Eqns. 2.5 and 2.6, E(r, t) and B(r, t) are therefore given by,E(r, t) = ıXωk Ak uk (r)e−ıωk t − A∗k u∗k (r)eıωk t ,(2.15)


Ak (∇ × uk (r))e−ıωk t + A∗k (∇ × u∗k (r))eıωk t .

(2.16)k

B(r, t) =XkEqns. 2.15 and 2.16 may be substituted into Eq. 2.14, and after evaluating the integrals using
Eq. 2.13, Eq. 2.14 becomes,
H=

X

Hk ,

(2.17)kwhere
Hk = 0 V ωk2 (Ak A∗k + A∗k Ak ) .

(2.18)

After defining the variables qk and pk as,
Ak =
A∗k =

1
√

(ωk qk + ıpk )2ωk 0 V
1√(ωk qk − ıpk ),2ωk 0 V(2.19)

Eq. 2.18 becomes,
1
Hk = (pk + ωk2 qk ).
2(2.20)

10This Hamiltonian has the form as a harmonic oscillator with position qk and its conjugate momentum pk , with unit mass. Thus each mode of the electromagnetic field can be formally mapped intoa harmonic oscillator degree of freedom.
The quantization procedure is now preformed in the typical way of quantum mechanics bydemanding that the variables qk and pk become quantum mechanical operators,
qk → q̂k , pk → p̂k ,(2.21)[q̂k , p̂k ] = ı~δk,k0 .

(2.22)

that obey the commutation relations,

Then, Eq. 2.20 becomes the Hamiltonian for the quantum harmonic oscillator,
1Ĥk = (p̂k + ωk2 q̂k ).
2(2.23)

This quantization procedure turns each variable Hk into an operator Ĥk that acts on an infinite
dimension Hilbert space. Therefore, the total Hamiltonian Ĥ is a sum over each one of these Hilbert
spaces with each individual space corresponding to a mode of the electromagnetic field. Formally,
to sum over Hilbert spaces means to construct the following:Ĥ = Ĥk1 ⊗ 1̂ ⊗ 1̂ ⊗ 1̂ ⊗ · · ·+1̂ ⊗ Ĥk1 ⊗ 1̂ ⊗ 1̂ ⊗ · · ·+··· ,(2.24)

where 1̂ is the identity operator on each Hilbert space. With this implied meaning, Eq. 2.24 iswritten simply as,
Ĥ =

X

Ĥk .

(2.25)k

It is convenient to work with annihilation operators âk and creation operators â†k which are
dimensionless and defined by,âk =â†k =

1(ωk q̂k + ıp̂k ),2~ωk1√(ωk q̂k − ıp̂k ),2~ωk√(2.26)

11
where the commutation relations for âk and â†k follow from Eq. 2.22, and are given by,hiâk , â†k = δk,k0 .

(2.27)The Hamiltonian for each mode indexed by different values of k can now be written as,1Ĥk = ~ωk âk â†k + â†k âk ,2(2.28)


1†.
Ĥk = ~ωk âk âk +2(2.29)

or by using Eq. 2.27,

The eigenbasis of Eq. 2.29 are the quantum harmonic oscillator number states, which obey,√âk |nik =√â†k |nik =n |n − 1ik
n + 1 |n + 1ik .

(2.30)

with âk |nik = 0 and integer n ≥ 0. It is now obvious that the above quantization procedure isequivalent to letting,ss~âk , A∗k →2ωk 0 V

Ak →~â†k .
2ωk 0 V(2.31)

Therefore, the electric and magnetic field operators, Ê(r, t) and B̂(r, t) are given by,Ê(r, t) = ı

X

r

k


~ωk 
âk uk (r)e−ıωk t − â†k u∗k (r)eıωk t ,20(2.32)~ âk (∇ × uk (r))e−ıωk t + â†k (∇ × u∗k (r))eıωk t .
2ωk 0(2.33)sB̂(r, t) =XkEqns. 2.32 and 2.33 are operators in the Heisenberg picture, since they contain the explicit timedependence. Specifically, they include the time dependence of the free evolution of the field underEq. 2.25. This can be seen by considering the frame transformation,
ı

ı

ı

ı

e ~ Ĥt âk e− ~ Ĥt = âk e−ıωk t
e ~ Ĥt â†k e− ~ Ĥt = â†k eıωk t .

(2.34)

12
Therefore, in the Schrödinger picture the fields may be written as,r
X ~ωk Ê(r, t) = ıâk uk (r) − â†k u∗k (r) ,
20
k
s

X
~ 
B̂(r, t) =âk (∇ × uk (r)) + â†k (∇ × u∗k (r)) .
2ωk 0(2.35)(2.36)k

For a single mode cavity field of frequency ωc = ck, oriented in the z direction satisfying standing
wave boundary conditions, Eq. 2.11 becomes,ru(z) =

2sin(kz),V(2.37)

where V is the volume of the cavity. Since we will be mainly focused on the properties of Ê(z, t)rather than B̂(z, t), it will be convenient to us to put Ê(z, t) in its simplest form. We therefore letâ → −ıâ in the expressions for Ê(z, t) and B̂(z, t). Then, for the case of linear polarization in thex direction, Eqns. 2.32 and 2.33 become,rx~ωc  −ıωc tâe+ â† eıωc t sin(kz)Ê(z, t) =
V 0
|x|
r


~ωc µ0yB̂(z, t) = −ıâe−ıωc t − â† eıωc t cos(kz) ,V|y|(2.38)(2.39)√
yx
where |x|and |y|are unit vectors in the direction x and y directions, and the fact that c = 1/ µ0 0has been used. The Hamiltonian Ĥ for a single mode cavity field is,1†Ĥc = ~ωc â â +.
2(2.40)

Since the zero energy definition is arbitrary, the factor of ~ωc /2 can be dropped from Eq. 2.40
without loss of generality, resulting in,Ĥc = ~ωc â† â.

(2.41)In the Schrödinger picture, where the time dependence of Ĥc will appear on to the state vector of
the system under consideration, instead of on â, Eqns. 2.38 and 2.39 are,
r

~ωc 
x
Ê(z, t) =â + â† sin(kz)
V 0
|x|
r


~ωc µ0yB̂(z, t) = −ıâ − â† cos(kz) .
V|y|(2.42)(2.43)

13

2.2

Jaynes Cummings Model
The Jaynes Cummings model describes a single two level atom interacting with a single

quantized cavity field mode. Besides being a prototypical cavity QED system, it is used to describe
the Hamiltonian part of the open quantum system that will later be considered. We will find that,in the Jaynes Cummings model, the atom oscillates between its excited and ground states, similar tothe case where the cavity field is treated classically. However, it differs from the classical case, sincethese oscillations can occur even when there are no photons present. Although these oscillations arereversible, this model provides good motivation for describing the irreversible spontaneous processes
such as spontaneous emission, which are caused by the coupling of a system to the many modes ofthe electromagnetic field in vacuum.

2.2.1

Jaynes Cummings Hamiltonian

The Hamiltonian describing the free evolution of a two level atom is given by,
Ĥa =

~ωa
(|ei he| − |gi hg|) ,2(2.44)

where |ei and |gi are the excited and ground states of the atom, which form an orthonormal basis,
and where ~ωa is the energy difference between the two levels.
The energy of an atom with dipole moment operator d̂ interacting with a single mode electromagnetic field in a cavity Ê(R, t) is described by,ĤI = −d̂ · Ê(R, t),(2.45)

where R is the position of the center of mass of the atom. It is valid to evaluate the electricfield Ê(R, t) at R since we are interested in electromagnetic fields with wavelengths much larger
than atomic dimensions. Therefore, in this approximation, known as the dipole approximation, theelectric field is treated as constant over the dimensions of an atom. After inserting the identity

14
operator 1̂ = (|ei he| + |gi hg|) on both sides of d̂ in Eq. 2.45, it becomes,ĤI= − (|ei he| + |gi hg|) d̂ (|ei he| + |gi hg|) · Ê(R, t)= − he| d̂ |gi |ei hg| + hg| d̂ |ei |gi he| · Ê(R, t),(2.46)

where he| d̂ |ei = hg| d̂ |gi = 0. This is because under transformation by the parity operator Π̂, thedipole operator becomes Π̂† d̂Π̂ = −d̂, and since Π̂ |ei = − |ei and Π̂ |gi = − |gi, it follows thathe| d̂ |ei = − he| d̂ |ei = hg| d̂ |gi = − hg| d̂ |gi = 0. After substituting Eq. 2.38, Eq. 2.46 becomes,ĤI= ~g (|ei hg| + |gi he|) âe−ıωc t + â† eıωc t .

(2.47)

Here, the coupling rate g between the atom and field is given by,rωcxsin(kz0 ) he| d̂ ·
|gi ,
g≡−
~V 0
|x|

(2.48)

where z0 is the position of the center of mass of the atom.

2.2.2

The Rotating Wave Approximation

We move into an interaction picture that puts the time dependence of the free evolution of
the atom on the atomic operators, defined by,
ı
ı
^
|ei
hg| ≡ e ~ Ĥa t |ei hg| e− ~ Ĥa t = |ei hg| e−ıωa t
ı
ı
^
|gi
he| ≡ e ~ Ĥa t |gi he| e− ~ Ĥa t = |gi he| eıωa t ,(2.49)

where the tilde over the operators indicates that they are in an interaction picture. After substitution of Eqns. 2.49, Eq. 2.47 becomes,
ĤI




^
^
= ~g |eihg|eıωa t + |gi
he|e−ıωa t âe−ıωc t + â† eıωc t

^^= ~g â|ei
hg|eı(ωa −ωc )t + â† |ei
hg|eı(ωa +ωc )t^^+â|gihe|e−ı(ωa +ωc )t + â† |gi
he|e−ı(ωa −ωc )t .

(2.50)

Close to resonance, i.e. when ωc ≈ ωa , terms with ωa + ωc in the exponential oscillate much more
rapidly than terms with ωa − ωc in the exponential. These fast oscillating terms average to zero on

15
a time scale that is much faster than the time in which the slowly oscillating terms would average
to zero, and therefore the fast oscillating terms can be dropped. In this approximation, called therotating wave approximation, Eq 2.50 becomes,
ĤI



^
^
= ~g â|ei
hg|eı(ωa −ωc )t + â† |gi
he|e−ı(ωa −ωc )t .

(2.51)

Moving back to the Schrödinger picture, where the time dependence due to Eq. 2.40 is not on theoperator â and the time dependence due to Eq. 2.44 is not on the operator |gi he|, we have,ĤI= ~g â |ei hg| + â† |gi he| .

(2.52)

Adding together the Eqns. 2.41, 2.44 and 2.51, we arrive at the Jaynes-Cummings Hamiltonian,
ĤJC

2.2.3

= ~ωc â† â +~ωa(|ei he| − |gi hg|) + ~g â |ei hg| + â† |gi he| .
2(2.53)

Rabi Flopping

The Schrödinger equation with Eq. 2.53 as the Hamiltonian,ı~d|ψ(t)i = ĤJC |ψ(t)i ,dt(2.54)

can be solved by first noting that Eq. 2.53 only couples the state |ei |ni, where the atom is in the
excited state and the field has n excitations, to the state |gi |n + 1i, where the atom is in the groundstate and the field has n + 1 excitations. Therefore, ĤJC can be split up into independent two-state
manifolds, which can be described by the basis {|ei |ni , |gi |n + 1i}. The manifold of Eq. 2.53 withn excitations may be represented as the 2 × 2 matrix in the reduced subspace,



√
1
ω
+
nω
g
n
+
1
c

 2 a
ĤJC =~˙ √− 12 ωa + (n + 1)ωc .
g n+1(2.55)The time dependent state vector |ψ(t)i of the system is given by solving the Schrödinger equation
with this Hamiltonian,|ψ(t)i = c1 (t) |ei |ni + c2 (t) |gi |n + 1i ,(2.56)

16
where c1 (t) = (he| hn|) |ψ(t)i and c2 (t) = (hg| hn + 1|) |ψ(t)i. Projecting Eq. 2.56 on to |ei |ni and|gi |n + 1i, and using Eq. 2.54, yields,


√
d
1
ı c1 (t) =ωa + nωc c1 (t) + g n + 1c2 (t)
dt
2


√
1
d
− ωa + (n + 1)ωc c2 (t) + g n + 1c1 (t).
ı c2 (t) =dt
2(2.57)

We now assume exact resonance, ωa = ωc , and define an interaction picture as,c˜1 (t) = eıωc (n+1/2)t c1 (t)c˜2 (t) = eıωc (n+1/2)t c2 (t).

(2.58)In this interaction picture, Eqns. 2.57 become,√dc̃1 (t) = −ıg n + 1c̃2 (t)dt√dc̃2 (t) = −ıg n + 1c̃1 (t)dt(2.59)The solutions to Eqns. 2.59 with initial conditions c1 (0) = 1 and c2 (0) = 0 are,√c̃1 (t) = cos(g n + 1t)√c̃2 (t) = −ı sin(g n + 1t),(2.60)

so that the probability of the system being in the excited and ground states, respectively, are,√|c̃1 (t)|2 = cos2 (g n + 1t)√|c̃2 (t)|2 = sin2 (g n + 1t).

(2.61)

Eqns. 2.61 tell us that the dynamics of the system is to oscillate, or flop, between the excited√and ground state with Rabi frequency g n + 1. It is interesting to note that even when there are
initially no excitations in the field i.e. n = 0 to stimulate an excited atom to emit its photon, thesystem will exhibit Rabi flopping. Although the expectation value of the electromagnetic field in
the vacuum state |0i is zero,rh0| Ê(z, t) |0i =
= 0,



~ωc
x
h0| â + â† |0i sin(kz)V 0|x|(2.62)

17
the expectation value for the intensity of the field in the vacuum state is non-zero,
2


2
~ωc
h0| â + â† |0i sin2 (kz)V 0~ωcsin2 (kz).
V 0h0| Ê(z, t) |0i ==(2.63)

It can be seen from Eqns. 2.62 and 2.63 that even though the electromagnetic field in vacuum is
on average zero, it fluctuates. Therefore, the Rabi flopping that occurs when there are initiallyno excitations in the field is caused by fluctuations in the vacuum electromagnetic field. This is aconsequence of the quantum mechanical treatment of the electromagnetic field; an excited atom ina classical field that contains no energy would never decay.
Rabi flopping between excited and ground state is a reversible process, unlike spontaneousemission in which the photon leaves and never comes back. To correctly model the irreversible
dissipative processes that are observed in nature, such as spontaneous emission, the system must
be coupled to a continuum of modes. Models that include irreversible dissipative processes will bederived in Chapters 3 and 4.

2.3

Wiener-Khinchin Theorem
We will be frequently interested in information regarding the frequency content of the lightthat is produced by an open quantum system. The open quantum system that we consider may
for instance consist of atoms interacting with an electromagnetic field inside a cavity. The spectral
information about the cavity field, which is described by the annihilation operator â(t), can beobtained from the power spectral density S(ω), which is defined by,E
1D †âT (ω)âT (ω) ,T →∞ T

S(ω) ≡ lim(2.64)

where,Z T /2âT (ω) ≡dte−iωt â(t).

(2.65)−T /2

Eq. 2.64 tells us how the total energy that is contained in the cavity field is distributed as afunction of frequency ω. The power spectral density, described by Eq. 2.64, will often consist of a

18
distribution of frequencies around a central frequency. The width of this distribution, called thelinewidth, denoted by ∆ν, is measure of the spectral purity of the light emitted by the consideredsystem.
Since many of the solution methods to the equations of open quantum systems do not give
direct access to â(t), it will often be difficult to calculate the power spectral density from Eq. 2.64
directly. However, an alternative method to obtain Eq. 2.64 is possible. This is demonstrated bysubstituting Eq. 2.65 into Eq. 2.64, which yields,1
S(ω) = lim
T →∞ TZ T /2

Z T /2dtDE0
dt0 e−iω(t −t) â† (t)â(t0 ) .

After changing variables to τ ≡ t0 − t, Eq. 2.66 becomes,(Z
Z T /2
DE
0
1
S(ω) = limdτdte−iωτ â† (t)â(t + τ )T →∞ T−T−T /2−τ
)Z T Z T /2−τ
D
E
−iωτ
†
+
dt
dte
â (t)â(t + τ ).
0(2.66)

−T /2

−T /2

(2.67)−T /2

We assume that the system under consideration is in steady state, so that â† (t)â(t0 ) is independentof t, and only depends on τ . Then, Eq. 2.67 becomes,
Z
D
E
1 T
S(ω) = limdτ e−iωτ â† (t)â(t + τ ) (T − |τ |) .
T →∞ T −T
Z ∞DE=dτ e−iωτ â† (t)â(t + τ ) ,(2.68)−∞

or since S(ω) is real,Z ∞DEdτ e−iωτ â† (t + τ )â(t)
−∞
Z ∞
D
E
−iωτ
†
= 2Re
dτ e
â (t + τ )â(t) .

S(ω) =(2.69)

0

Eq. 2.69 is known as the Wiener-Khinchin Theorem. It says that the power spectrum is the Fouriertransform of the two-time correlation function â† (t + τ )â(t) . Therefore, the power spectrum S(ω)
and its linewidth ∆ν, can be obtained from â† (t + τ )â(t) , which in many cases can be directly
accessed.

Chapter 3

Master Equation

To correctly incorporate dissipation into a quantum mechanical system, the system is assumed
to be small, and weakly coupled to a much larger reservoir. Since we only care about keeping trackof system variables, the combined system and reservoir are evolved to second order in perturbation
theory, and then the reservoir variables are averaged over. The resulting equation, called the masterequation, keeps track of system variables only, and allows the system to lose or gain energy fromits weak coupling to the large reservoir.
This approach is in the Schrödinger picture, where the time dependence is put on the statevector of the combined system and reservoir. Because the reservoir will be averaged over, thesystem will evolve into a statistical mixture of states, and it will no longer be able to be describedby a single state vector. Instead, it will be necessary to describe the system by a reduced density
operator, which can account for the system evolving into a statistical mixture of states.

3.1

Density Operator
Operators corresponding to system observables do not act on the reservoir part of the total

Hilbert space. Therefore, if some sort of a partial expectation value over the reservoir part of thethe state vector describing the combined system and reservoir could be performed, the complexityof the problem would be drastically reduced from the start. Unfortunately, a partial expectation
value of a state vector can not be taken, but an analogous operation can be performed on a density

20
operator, which is defined asχ̂(t) ≡ |Ψ(t)i hΨ(t)| ,(3.1)

where |Ψ(t)i is the state vector for the combined system and reservoir.
Let θ̂SR be an arbitrary operator of the system and reservoir Hilbert spaces. The trace ofχ̂(t) and θ̂SR is given by,hiXX
T r χ̂(t)θ̂SR =hs| hr| χ̂(t)θ̂SR |si |ri
s

=

s

=

r

XX
XX
s

hs| hr|Ψ(t)i hΨ(t)| θ̂SR |si |rirhΨ(t)| θ̂SR |si |ri hs| hr|Ψ(t)ir!
= hΨ(t)| θ̂SRXX
s|si |ri hs| hr| |Ψ(t)ir= hΨ(t)| θ̂SR |Ψ(t)i ,(3.2)
where the states |si are a basis for the system Hilbert space, and the states |ri are a basis for the
reservoir Hilbert space. This tells us that performing the trace of χ̂(t) and θ̂SR is equivalent totaking the expectation value of θ̂SR .
If we are only interested observables associated with the system, we can preform partial traceover the reservoir Hilbert space, and arrive at a reduced density operator for the system only. Let
θ̂S , be an arbitrary system operator. The trace of χ̂(t) and θ̂S is,hiXX
T r χ̂(t)θ̂S =hs| hr| χ̂(t)θ̂S |si |ri
s

r

=

XX

=

X

s

hs| hr| χ̂(t) |ri θ̂S |si

r

!

=

hs|

X

hr| χ̂(t) |ri θ̂S |si

s

r

X

hs| ρ̂(t)θ̂S |si ,

s

(3.3)

21
where,
ρ̂(t) =

X

hr| χ̂(t) |rir= T rR [χ̂(t)] ,(3.4)
where ρ̂(t) is the reduced density operator of the system. From Eq. 3.3 it can be seen that if the
expectation value of a system operator is desired, it is valid to first take the partial trace of χ̂(t)over the reservoir Hilbert space to get ρ̂(t) and then later take the remaining partial trace of ρ̂(t)and θ̂ to get the total trace. In the next section, the equation of motion for a combined system andreservoir will be considered, and a partial trace over the reservoir Hilbert space will be preformed.
This will leave us with an equation of motion for the reduced density operator of the system. It isa drastic simplification to not have to keep track of the reservoir Hilbert space.
A consequence of performing a partial trace over the reservoir Hilbert space is that we will
no longer have complete information about combined system and reservoir. This causes the thereduced density operator to be a mixed state, or a statistical ensemble of states. It will therefore not
be possible to rewrite the reduced density it in terms of state vectors. Therefore, the reformulationof the problem in terms of density matrices is not only for convenience, it is essential.
d
The equation of motion for χ̂(t) that is equivalent to the Schrödinger equation, ı~ dt|Ψ(t)i =Ĥ(t) |Ψ(t)i, is called the Von Neumann equation. It is derived by taking a time derivative of χ̂(t)and substituting in the Schroödinger equation,ddχ̂(t) =(|Ψ(t)i hΨ(t)|)
dt
dt


d |Ψ(t)id hΨ(t)|=hΨ(t)| + |Ψ(t)i
dt
dt


1
=
Ĥ(t) |Ψ(t)i hΨ(t)| − |Ψ(t)i hΨ(t)| Ĥ(t)
ı~
i
1 h
=
Ĥ(t), χ̂(t) .
ı~
(3.5)

22

3.2

The Master Equation for a Damped Harmonic Oscillator
To derive the master equation [67], we will need to make two approximations. First, theinteraction between system and reservoir is assumed to be weak, so that it is sufficient to evolvethe system to leading order in perturbation theory. This is the Born approximation. Also, thefrequency bandwidth of the reservoir is assumed to be very large, so that in the time domain,
reservoir correlations decay away very fast. In the limiting case of infinite bandwidth, the reservoir
is delta correlated in time. This is the Markov approximation.
In order to best illuminate the approximations and the physics involved, the master equation
is derived for a simple system instead of in full generality. Once our foot is in the door, it is easy togeneralize to more complicated systems, since the derivation for those systems follow in an almost
identical manner.

3.2.1

Derivation

Consider a quantum harmonic oscillator of frequency ωc , which is coupled to many otherquantum harmonic oscillators of frequencies ωj . This could be, for instance, used to model the thecoupling of a cavity mode to all other modes of the electromagnetic field in a system cavity QED
system. The Hamiltonian of the combined system and reservoir Ĥ, is given by,Ĥ = ĤS + ĤR + ĤSR ,(3.6)

where ĤS is the Hamiltonian for the system,
ĤS = ~ωc â† â,(3.7)

ĤR is the Hamiltonian for the reservoir,
ĤR = ~X

ωj r̂j† r̂j

(3.8)

j

and ĤSR is the Hamiltonian for the coupling between the system and reservoir,
ĤSR = ~Xj
λj â† r̂j + λ∗j âr̂j† ,

(3.9)

23
where â is the annihilation operator in the system Hilbert space and r̂j is the annihilation operatorof the jth mode in the reservoir Hilbert space, and λj is the coupling between the system and thejth mode in the reservoir.
At time t = t0 , it is assumed that there are no correlations between the system and reservoir,
and therefore χ̂(t0 ) is of the separable form,χ̂(t0 ) = ρ̂(t0 ) ⊗ R̂,(3.10)

where ρ̂(t) is the reduced density operator of the system, and R̂ is the reduced density operator of
the reservoir, which are defined respectively as,ρ̂(t) ≡ trR {χ̂(t)} ,(3.11)

R̂ ≡ trS {χ̂(t)} .

(3.12)

andThe reservoir density operator, R̂, has no time dependence because we are assuming that thereservoir is broadband, so that in the time domain, it returns to its steady state value on a timescale much faster than system dynamics. Therefore at later times we can also assume that χ̂(t)will be seperable,χ̂(t) = ρ̂(t) ⊗ R̂.

(3.13)The reservoir is assumed to be in thermal equilibrium at temperature T , which is described by,−ĤR

e kB T Ĥ  .
R̂ =− R
trR e kB T

(3.14)

The first step in our derivation is to move into an interaction picture, so that the timedependence of ĤS and ĤR can be rotated out. This is done by the following transformation:
˜ = e ~i (ĤS +ĤR )t χ̂(t)e− ~i (ĤS +ĤR )t .
χ̂(t)(3.15)

It is worth noting that Eq. 3.14 implies that R̂ is unaffected by this transformation, since anoperator commutes with an exponential of itself, so that,˜R̂ = R̂

(3.16)

24
To find the time evolution in the interaction picture, the derivative of both sides of Eq. 3.15 is firstfound,
 i
i
1 
˜˙
χ̂(t)
= −
ĤS + ĤR e ~ (ĤS +ĤR )t χ̂(t)e− ~ (ĤS +ĤR )t
i~

 

 i
i
1
− ~i (ĤS +ĤR )t˙e− ~ (ĤS +ĤR )t . (3.17)+ χ̂(t)ĤS + ĤR
+e ~ (ĤS +ĤR ) χ̂(t)ei~
Upon substitution of Eq. 3.5 into Eq. 3.17, it becomes,

1 
˜
˜˙
ĤS + ĤR χ̂(t)
χ̂(t)
= −
i~
 h
i i

i
1
1 ˜ 
ĤS +ĤR )(~+e
Ĥ(t), χ̂(t) e− ~ (ĤS +ĤR )t + χ̂(t)ĤS + ĤR .
i~i~(3.18)

Finally, we now use Eq. 3.6 and simplify, and arrive at,



1 
˜ + 1 χ̂(t)˜˜˙ĤS + ĤR χ̂(t)ĤS + ĤR
χ̂(t)
= −
i~
i~



 i
i
1
+e ~ (ĤS +ĤR )ĤS + ĤR + ĤSR (t) χ̂(t) − χ̂(t) ĤS + ĤR + ĤSR (t) e− ~ (ĤS +ĤR )t
i~
 i

i
1
Ĥ
+
Ĥ
= e ~ ( S R)
ĤSR (t)χ̂(t) − χ̂(t)ĤSR (t) e− ~ (ĤS +ĤR )ti~ii
1 i (ĤS +ĤR )ĤSR (t)e− ~ (ĤS +ĤR ) e ~ (ĤS +ĤR ) χ̂(t)
=
e~
i~
!
i
− ~i (ĤS +ĤR ) ~i (ĤS +ĤR )−χ̂(t)eeĤ (t) e− ~ (ĤS +ĤR )t
SR

=

i
1 h˜
˜
ĤSR (t), χ̂(t),i~(3.19)

where,ii
˜ĤSR (t) ≡ e ~ (ĤS +ĤR )t ĤSR (t)e− ~ (ĤS +ĤR )t .

(3.20)

Next, both sides of Eq. 3.19 are integrated,˜ − χ̂(t˜ 0) = 1
χ̂(t)
i~

Z t

h
i
˜
˜ 0) ,
dt0 ĤSR (t0 ), χ̂(t(3.21)

t0

and then substitute this back into itself to obtain,˜ − χ̂(t˜ 0) = 1
χ̂(t)
i~

Z t
Z t0
i
h
h
ii
1
˜
˜
˜
0 ˜
0
˜ 00 ) . (3.22)
dt ĤSR (t ), χ̂(t0 ) − 2dtdt00 ĤSR (t0 ), ĤSR (t00 ), χ̂(t
~ t0
t0
t0

Z t

0

h

25
After taking the partial trace over the reservoir of both sides of Eq. 3.22, it is given by,
Z t0
Z t
iio
h
nh
1
˜
˜
0
˜ 00 )˜˜,
dt00 trR ĤSR (t0 ), ĤSR (t00 ), χ̂(t
ρ̂(t) − ρ̂(t0 ) = − 2
dt
~ t0
t0

(3.23)

where the first term on the right hand side of Eq. 3.22 vanishes since the trace over the reservoirof r̂j and r̂j† are zero. The vanishing of this term is why Eq. 3.23 was expanded up to second order˜
in ĤSR . The Born approximation is now made, by neglecting terms higher than this leading order˜
in ĤSR .
We assume that the integration time interval τ ≡ t − t0 is short compared to the time it takes
the reservoir to change, but long compared to the time it takes the system to change. Under thisassumption, it is valid to approximate,˜ 0 + τ ) − ρ̂(t
˜ 0)
ρ̂(t
d˜
).
≈ ρ̂(τ
τdt(3.24)

˜
Under this approximation, Eq 3.23 is turned into a coarse-grained equation of motion for ρ̂(t),
Z τ0
Z τ
nh
h
iio
d ˜
1
˜
˜
0
˜
dτ 00 trR ĤSR (τ 0 ), ĤSR (τ 00 ), χ̂(t),(3.25)
ρ̂(τ ) = − 2
dτ
dτ
~ τ 0
0
where the independent variable has been changed from t to τ , and the time origin has been shiftedby t0 .
In order to proceed, the form of ĤSR in the interaction picture must be found. To do this,the time dependence of the system operator â and reservoir operators r̂j in the interaction picture
are found,
i
i
˜
â(t)
= e ~ (ĤS +ĤR )t âe− ~ (ĤS +ĤR )tii= e ~ ĤS t âe− ~ ĤS t= âe−iωc t ,(3.26)iir̂˜j (t) = e ~ (ĤS +ĤR )t r̂j e− ~ (ĤS +ĤR )tii= e ~ ĤR t r̂j e− ~ ĤR t= r̂j e−iωj t .

(3.27)

26
Eqns. 3.26 – 3.27, and their Hermitian conjugates are substituted into Eq. 3.9, so that ĤSR
in the interaction picture is given by,X˜ĤSR (τ ) = ~λj ei(ωc −ωj )τ â† r̂j + λ∗j e−i(ωc −ωj )τ âr̂j† ,j= ı~ â† F̂ (τ ) − âF̂ † (τ ) ,(3.28)

where the reservoir operator F̂ (τ ) is defined as,F̂ (τ ) ≡ −ıX

λj r̂j ei(ωj −ωc )τ .

(3.29)j

After expanding the double commutator, and substituting in Eq. 3.28, Eq. 3.25 becomes,d˜ρ̂(t) =
dt

1
τ

Z τ

dτ 0

Z τ0

(
dτ 00 trRâ† F̂ (τ 0 ) − âF̂ † (τ 0 )˜â† F̂ (τ 00 ) − âF̂ † (τ 00 ) χ̂(t)

0

0






˜
− â† F̂ (τ 0 ) − âF̂ † (τ 0 ) χ̂(t)â† F̂ (τ 00 ) − âF̂ † (τ 00 )
=

1
τ

Z τ
dτ
0

0

Z τ0

)
+ H.C.

n
o
n
o
0
00
† ˜
0
† 00
˜
dτ 00 â† â† ρ̂(t)tr
R F̂ (τ )F̂ (τ )R̂ − â âρ̂(t)trR F̂ (τ )F̂ (τ )R̂

0

n
o
n
o
† 0
00
† 0
† 00
˜
˜
−ââ† ρ̂(t)tr
R F̂ (τ )F̂ (τ )R̂ + ââρ̂(t)trR F̂ (τ )F̂ (τ )R̂
o
o
n
n
0
† 00
†
˜
˜
F̂
(τ
)
R̂
F̂
(τ
)
−â† ρ̂(t)âtrR F̂ (τ 0 )R̂F̂ (τ 00 ) + â† ρ̂(t)âtr
R
n
o
n
o
† 0
† 00
†
˜
˜
+âρ̂(t)âtrR F̂ † (τ 0 )R̂F̂ (τ 00 ) − âρ̂(t)âtr
R F̂ (τ )R̂F̂ (τ )!
+ H.C.

(3.30)

Since R̂ is diagonal in the reservoir eigenbasis, terms with two F̂ (τ ) or two F̂ † (τ ) in thereservoir trace will be zero. However, terms with one F̂ (τ ) and one F̂ † (τ ) in the reservoir trace will
be non-zero. Consider, for example, the term proportional to,

Z τ
0

dτ 0

Z τ0
0

Z τ
Z τ0
n
o
n
o
X
0
00
dτ 00 trR F̂ † (τ 0 )F̂ (τ 00 )R̂
=
dτ 0
dτ 00
λ?i λj trR r̂i† r̂j R̂ e−i(ωi −ωc )τ ei(ωj −ωc )τ
0

0

Z τ

Z τ0

=
0

dτ 0

0

i,j

dτ 00

X

o
n
0
00
|λj |2 trR r̂j† r̂j R̂ ei(ωc −ωj )(τ −τ ) ,j(3.31)

27
where the last line follows again from the fact that R̂ is diagonal in the reservoir eigenbasis.
The sum is changed to an integral over ω by including a density of states D(ω),
Z τ0
Z τ
Z τ0
Z τ
o
n
o
n
X
0
00
dτ 00 trR r̂j† r̂j R̂ ei(ωc −ωj )(τ −τ )dτ 0=∆j |λj |2dτ 00 trR F̂ † (τ 0 )F̂ (τ 00 )R̂
dτ 0
0

0

0

0

j

Z B

2

Z τ

dωD(ω) |λ(ω)|

=

dτ

Z τ0

n
o
0
00
dτ 00 trR r̂† r̂R̂ ei(ωc −ω)(τ −τ )

0

0

0

0

(3.32)
where ∆j = 1 and B is the bandwidth of the reservoir. Since we assuming that the bandwidth of
the reservoir is very large, the reservoir relaxation rate τ1R will be much larger than the systemsnorelaxation rate τ1S . Then, for τ << τS , trR r̂† r̂R̂ will be approximately constant, so that it canbe factored outside of the integral,
Z τ
Z τ0
o
n
0
dτ
dτ 00 trR F̂ † (τ 0 )F̂ (τ 00 )R̂
0

0

Z B

Z τ0
n
oZ τ
0
00
†
0
dωD(ω) |λ(ω)| trR r̂ r̂R̂dτdτ 00 ei(ωc −ω)(τ −τ ) .

0

0

=

2

(3.33)

0

The resulting integral of the exponential will then destructively interfere when T ≡ τ 0 − τ 00 >> τR ,resulting in it averaging to zero. Therefore, after changing variables to T , the upper limit ofintegration on the T integral can be extended to T = ∞ with minimal error, so that Eq. 3.33
becomes,
Z τ0
Z τ
n
o
0
dτ 00 trR F̂ † (τ 0 )F̂ (τ 00 )R̂
dτ
0
0
Z B
Z ∞
n
oZ τ
=
dωD(ω) |λ(ω)|2 trR r̂† r̂R̂dτ 0dT ei(ωc −ω)T
0
0
0


Z B
n
oZ τ
2
†
0
=
dωD(ω) |λ(ω)| trR r̂ r̂R̂dτ πδ(ω − ωc ) − iP00where we have used,
Z t
Z t−t00 i(ωc −ω)(t−t0 )i(ωc −ω)T
lim
dt e
= lim
dT e
= πδ(ω − ωc ) − iP
t→∞ t
0

t→∞ 0

1
ω − ωc1ω − ωc, (3.34)

.

(3.35)

Neglecting the principal part, P (1/(ω − ωc )), which leads to a negligible frequency shift, Eq. 3.34
becomes,
Z τ
0

dτ 0

Z τ0
0

n
o
κτ n̄
dτ 00 trR F̂ † (τ 0 )F̂ (τ 00 )R̂=,2(3.36)

28
where,
κ = 2πD(ωc ) |λ(ωc )|2 ,(3.37)and n̄ is the average number of thermal excitations in the reservoir,non̄ = trR r̂† r̂R̂ .

(3.38)

After evaluation of the other reservoir trace terms, which follow in almost an identical manner tothe one demonstrated, Eq. 3.30 becomes the master equation for the damped harmonic oscillator,d˜ρ̂(t) =dt

1 † ˜†˜˜˜−â âρ̂(t)κ (n̄ + 1) − ââ† ρ̂(t)κn̄+ â† ρ̂(t)âκn̄+ âρ̂(t)âκ (n̄ + 1) + H.C.
2

 κ 

κ
†
†
†
˜ + ρ̂(t)â˜˜˜ + ρ̂(t)ââ˜˜= − (n̄ + 1) â† âρ̂(t)â − 2âρ̂(t)â− n̄ ââ† ρ̂(t)− 2â† ρ̂(t)â.
2
2(3.39)

3.2.2

Analytic Solutions for Photon Number

Eq. 3.39 can be used to calculate the average photon number â† â since in the Schrödinger
picture,d D † E† d ˜â â = tr â â ρ̂(t) .
dt
dt(3.40)

Now, Eq. 3.39 is substituted into Eq. 3.40 to obtainihd D † Eκ††˜ + â† âρ̂(t)â
˜
˜
â â
= − (n̄ + 1) tr â† ââ† âρ̂(t)â − 2â† ââρ̂(t)â
dt
2
h
i
κ
†
˜ + â† âρ̂(t)ââ˜˜− n̄tr â† âââ† ρ̂(t)− 2â† ââ† ρ̂(t)â2hiκ˜ + â† ââ† âρ̂(t)˜ − 2â† â† ââρ̂(t)˜= − (n̄ + 1) tr â† ââ† âρ̂(t)2hiκ˜ + ââ† â† âρ̂(t)˜ − 2ââ† ââ† ρ̂(t)˜− n̄tr â† âââ† ρ̂(t),2(3.41)

where in the last line, we have used the cyclic property of the trace,
h
i
h
i
h
i
tr ÂB̂ Ĉ = tr Ĉ ÂB̂ = tr B̂ Ĉ Â ,(3.42)

with Â, B̂, and Ĉ being arbitrary operators. Eq. 3.41 can be simplified by using the commutation
relations for creation and annihilation operators,hiâ, â† = ââ† − â† â = 1̂,(3.43)

29
where 1̂ is the identity operator. We find,hid D † Eκ˜ − 2â† â† ââρ̂(t)˜â â= − (n̄ + 1) tr 2â† â† â + 1 âρ̂(t)
dt
2
h



i
κ
˜ + ââ† ââ† − 1 ρ̂(t)˜ − 2ââ† ââ† ρ̂(t)˜− n̄tr ââ† − 1 ââ† ρ̂(t)2D E= −κ â† â + +κn̄.

(3.44)The solution to Eq. 3.44 is,
D

E
D E

â† â (t) = â† â (0)e−κt + n̄ 1 − e−κt .

(3.45)

Eq. 3.45 will prove useful in sections 3.4 and 3.5 when we will wish to compare the numericalsolution methods demonstrated in that section to an analytic solution. These numerical solution
can be used to treat more complex systems, where analytic solutions are not possible.

3.2.3

The Quantum Regression Theorem

Consider the two-time expectation value of two arbitrary system operators Â and B̂,
D

E
h i
i
i
Â(t + τ )B̂(t) = tr e ~ Ĥτ Â(t)e− ~ Ĥτ B̂(t)ρ̂ ,(3.46)

where Ĥ is the total Hamiltonian of system and reservoir. Using the cyclic property of the trace,
we findDÂ(t + τ )B̂(t)

E

h i
i
i
= tr e− ~ Ĥτ B̂(t)ρ̂e ~ Ĥτ Â(t)hi= tr Ŝ(τ )Â(t) ,(3.47)

where,i

i

Ŝ(τ ) ≡ e− ~ Ĥτ B̂(t)ρ̂e ~ Ĥτ .

(3.48)

Eq. 3.48 says that Ŝ(τ ) obeys the same equation of motion as ρ̂. Since ρ̂ obeys the master equation,
so does Ŝ(τ ). Taking the derivative with respect to τ of both sides of Eq. 3.47 yields,E
d DÂ(t + τ )B̂(t)
= tr
dτ





d
Ŝ(τ ) Â(t) .
dτ

(3.49)

30dThus, since dτŜ(τ ) is known, the equation of motion for the two-time correlation function
DE
Â(t + τ )B̂(t) can be calculated. This is the the quantum regression theorem.
To demonstrate the quantum regression theorem, â† (t + τ )â(t) is calculated for the damped
harmonic oscillator. We let Â → â† and B̂ → â in Eq. 3.49,



E
d D †
d
†
â (t + τ )â(t)=Ŝ(τ ) â (t)dτdτ

(3.50)

Using Eq. 3.39 for Ŝ(τ ) yields,
E
h
i
κ
d D †
â (t + τ )â(t)= − (n̄ + 1) tr â† âŜ(τ )â† + Ŝ(τ )â† ââ† − 2âŜ(τ )â† â†
dτ
2
h
i
κ
− n̄tr ââ† Ŝ(τ )â† + Ŝ(τ )ââ† â† − 2â† Ŝ(τ )ââ†2hiκ= − (n̄ + 1) tr â† â† âŜ(τ ) + â† ââ† Ŝ(τ ) − 2â† â† âŜ(τ )2ihκ− n̄tr â† ââ† Ŝ(τ ) + ââ† â† Ŝ(τ ) − 2ââ† â† Ŝ(τ ) ,2(3.51)

where in the last line, the cyclic property of the trace has been used. The commutation relations
for creation and annihilation operators are used to simplify, yielding,


i
E
h
d D †
κ
â (t + τ )â(t)= − (n̄ + 1) tr â† â† âŜ(τ ) + â† â† â + 1̂ Ŝ(τ ) − 2â† â† âŜ(τ )
dτ
2
h

i
κ
− n̄tr ââ† − 1̂ â† Ŝ(τ ) + ââ† â† Ŝ(τ ) − 2ââ† â† Ŝ(τ )
2
h
i κ
h
i
κ
= − (n̄ + 1) tr â† Ŝ(τ ) + n̄tr â† Ŝ(τ )
2
2
i
κ h †
= − tr â Ŝ(τ ) .
(3.52)
2
Eq. 3.48 is substituted with B̂ → â, and the cyclic property of the trace is again used. Then,
E
i
i
i
d D †
κ h
â (t + τ )â(t)= − tr â† e− ~ Ĥτ B̂(t)ρ̂e ~ Ĥτdτ2
κ h i Ĥτ † − i Ĥτ i= − tr e ~ â e ~ âρ̂
2
i
κ h †
= − tr â (t + τ )â(t)ρ̂(t)
2
E
κD †
= −
â (t + τ )â(t) .
2(3.53)The solution to Eq. 3.53 is,D

E DEâ† (t + τ )â(t) = â† (t)â(t) e−κτ /2 ,(3.54)

31
and from Eq. 3.45,
D

3.3

E D Eâ† (t + τ )â(t) = â† â (0)e−κt + n̄ 1 − e−κt e−κτ /2 .

(3.55)The Master Equation for a Damped Two Level Atom
The master equation for a two level atom coupled to a reservoir of harmonic oscillators

follows almost instantaneously. This system could, for instance, describe the spontaneous emissionof radiation from a two level atom due to its coupling to the many modes of the electromagnetic
field. Describing how two level atoms are damped by a reservoir will be essential to systems thatwill be treated later.
The Hamiltonian of the combined system and reservoir, Ĥ, is again given by Eq. 3.6, butnow, the system Hamiltonian ĤS is given by,
ĤS =

~ωa z
σ̂ ,
2

(3.56)

where ωa is the frequency difference between the two atomic levels. The Hamiltonian for the
coupling between the system and reservoir, ĤSR , is now given by,
ĤSR = ~Xλj σ̂ + r̂j + λ∗j σ̂ − r̂j† ,(3.57)jwhere r̂j is the annihilation operator of the jth mode in the reservoir Hilbert space, λj is coupling
rate between the system and the jth mode of the reservoir, and σ̂ − , σ̂ + , and σ̂ z are the Pauli
operators, which obey, + −σ̂ , σ̂= σ̂ z z +σ̂ , σ̂= 2σ̂ +
 − −σ̂ , σ̂= 2σ̂ − .

(3.58)

32
To find Eq. 3.57 in the interaction picture, σ̂ − in the interaction picture must be found,˜ − (t) = e ~i (ĤS +ĤR )t σ̂ − e− ~i (ĤS +ĤR )tσ̂ii= e ~ ĤS t σ̂ − e− ~ ĤS t= σ̂ − e−iωa t ,(3.59)

Eqns. 3.26 and 3.27, and their Hermitian conjugates are substituted into Eq. 3.57, so that
ĤSR in the interaction picture is given by,X˜ĤSR (τ ) = ~λj ei(ωc −ωj )τ σ̂ + r̂j + λ∗j e−i(ωc −ωj )τ σ̂ − r̂j† ,j= ı~ σ̂ + F̂ (τ ) − σ̂ − F̂ † (τ ) ,(3.60)

where F̂ (τ ) is defined identically to the way it was defined in Eq. 3.29. After making the replacements â → σ̂ − and â† → σ̂ + , Eq. 3.60 has the exact same form as Eq. 3.28, so the derivation is
identical, except now the decay rate for the atoms, γ, is given by,
γ = 2πD(ωa ) |λ(ωa )|2 ,(3.61)The master equation for the damped two level atom is therefore,


γ
d˜
+ −
+
˜ + ρ̂(t)σ̂˜˜ρ̂(t) = − (n̄ + 1) σ̂ + σ̂ − ρ̂(t)σ̂ − 2σ̂ − ρ̂(t)σ̂
dt
2

γ 
− +
−
˜ + ρ̂(t)σ̂˜˜− n̄ σ̂ − σ̂ + ρ̂(t)σ̂ − 2σ̂ + ρ̂(t)σ̂
.
2

3.4

(3.62)

Quantum Jumps
The master equation can be solved numerically, without any approximations, by using thequantum jumps algorithm. In this algorithm, many random pieces, or trajectories, of the Dyson
expansion of the master equation are found, and they are then averaged over. Because these
trajectories are chosen with the correct weights to reproduce the Dyson expansion of the masterequation, the correct evolution is given. The Dyson expansion of the master equation, how to findthe correct weights of specific trajectories, and how to preform the algorithm are described.

33
3.4.1

The Dyson Expansion of the Master Equation

A general master equation can be written in the form,
N n
i 1X
o
d
ıh
ρ̂(t) = − Ĥ0 , ρ̂(t) −Ŝk+ Ŝk− ρ̂(t) + ρ̂(t)Ŝk+ Ŝk− − 2Ŝk− ρ̂(t)Ŝk+ ,dt~
2(3.63)

k=1

where Ĥ0 describes the coherent evolution of the system, and
Ŝk =

√

γk ŝk ,(3.64)

is a jump operator, where ŝk is a system operator, and γk is a system decay rate.
We define an effective Hamiltonian asĤef f ≡ Ĥ0 −ı~ X †Ŝk Ŝk ,2(3.65)

kso that Eq. 3.63 becomes
 X
d
ı
†
Ŝk† ρ̂(t)Ŝk .
ρ̂(t) = −Ĥef f ρ̂(t) − ρ̂(t)Ĥef
f +
dt
~

(3.66)k

Eq. 3.66 can now be written in the formdρ̂(t) = Lˆρ̂(t),dt(3.67)

Lˆρ̂(t) = Lˆ0 ρ̂(t) + LˆJ ρ̂(t),(3.68)

ı
†Ĥef f ρ̂(t) − ρ̂(t)Ĥeff ,~X
LˆJ ρ̂(t) =Ŝk ρ̂(t)Ŝk† .

(3.69)

withwhereL̂0 ρ̂(t) = −k

In an interaction picture defined by the transformation,
˜ = e−Lˆ0 t ρ̂(t),ρ̂(t)(3.70)d˜ˆˆ ˜ρ̂(t) = e−L0 t LˆJ eL0 t ρ̂(t).
dt(3.71)

Eq. 3.67 reduces to34
Both sides of Eq. 3.71 are integrated from t0 = 0 to t0 = t to obtain,˜ = ρ̂(0)
˜ +
ρ̂(t)

Z t

ˆ 0
ˆ 0˜ 0dt0 e−L0 t LˆJ eL0 t ρ̂(t).

(3.72)

0

The above equation is substituted into its self to yield˜ = ρ̂(0)
˜ +
ρ̂(t)

Z t

0 −Lˆ0 t0

dt e

ˆ 0˜
LˆJ eL0 t ρ̂(0)
+

Z t
dt

0

Z t0

(3.73)

0

0

0

ˆ 0
ˆ 0 00ˆ 00 ˜ 00
dt00 e−L0 t LˆJ eL0 (t −t ) LˆJ eL0 t ρ̂(t)˜
This procedure can be repeated an infinite number of times to yield the Dyson expansion for ρ̂(t),
˜
˜ +
ρ̂(t)
= ρ̂(0)Z t

0 −Lˆ0 t0

dt e

ˆ 0˜
LˆJ eL0 t ρ̂(0)
+

Z t
dt

d1

+

dt00

0

0

˜ +
= ρ̂(0)

Z t0

Z t00

ˆ 0
ˆ 0 00ˆ 00 ˜dt00 e−L0 t LˆJ eL0 (t −t ) LˆJ eL0 t ρ̂(0)

0

0

0

Z t0

Z t

0

ˆ 0
ˆ 0 00ˆ 00 000
ˆ 000 ˜dt000 e−L0 t LˆJ eL0 (t −t ) LˆJ eL0 (t −t ) LˆJ eL0 t ρ̂(0)
+ ···

0

∞ Z t
X

Z t3

Z tn
dtn−1 · · ·

dtn

dt2
0

0

n=1 0

Z t2
dt1
0

ˆ
ˆ
ˆ
ˆ ˜
×e−L0 tn LˆJ eL0 (tn −tn−1 ) LˆJ · · · LˆJ eL0 (t2 −t1 ) LˆJ eL0 t1 ρ̂(0).

(3.74)

Applying the inverse transformation that was applied in Eq. 3.70 yields the Dyson expansion in
the Schrödinger picture,
Lˆ0 t

∞ Z t
X

Z tn

Z t3

Z t2

dtn
dtn−1 · · ·
dt2
dt1
0
0
0
n=1 0
ˆ
ˆ
ˆ
ˆ
×eL0 (t−tn ) LˆJ eL0 (tn −tn−1 ) LˆJ · · · LˆJ eL0 (t2 −t1 ) LˆJ eL0 t1 ρ̂(0),ρ̂(t) = eρ̂(0) +

(3.75)

˜ = ρ̂(0), by definition. After plugging the more explicit form of the jump Liouvillian into
where ρ̂(0)Eq. 3.75, it becomes,ρ̂(t) = eLˆ0 tρ̂(0) +∞X
X Z t
n=1 k1 ,··· ,kn

0

Z tn

Z t3
dtn−1 · · ·

dtn
0

Z t2
dt2

0

dt1
0

ˆ
(k ) ˆ
(k)(k ) ˆ(k ) ˆ
×eL0 (t−tn ) LˆJ n eL0 (tn −tn−1 ) LˆJ n−1 · · · LˆJ 2 eL0 (t2 −t1 ) LˆJ 1 eL0 t1 ρ̂(0),(3.76)
where(k )LˆJ i = Ŝki ρ̂(t)Ŝk†i .

(3.77)

35
3.4.2

Definition of Quantum Trajectories

Justifying the name later, we define a ”trajectory” to be the quantity inside all of the sums
and integrals in Eq. 3.76,
ˆ

ˆ

ˆ

ˆ

ρ̂t1 ,k1 ;t2 ,k2 ;··· ;tn ,kn (t) = eL0 (t−tn ) LˆJ n eL0 (tn −tn−1 ) LˆJ n−1 · · · LˆJ 2 eL0 (t2 −t1 ) LˆJ 1 eL0 t1 ρ̂(0). (3.78)
(k )

(k

)

(k )

(k )

A trajectory depends explicitly on the values of ki and the values of ti that appear in its definition.
In a trajectory, there are only two types of operations being performed. One is the free evolution
under Ĥef f ,ˆρ̂(ti ) = eL0 (ti −ti−1 ) ρ̂(ti−1 ),(3.79)

which is the solution to the equation,dρ̂(ti ) = Lˆ0 ρ̂(ti )
dti

1 
†
=
Ĥef f ρ̂(ti ) − ρ̂(ti )Ĥeff ,ı~(3.80)and the other is the application of a jump operator,(k)
LˆJ ρ̂(ti ) = Ŝk ρ̂(ti )Ŝk† .

(3.81)

If the initial condition for the density operator is chosen to be in a pure state,ρ̂(0) = |ψ(0)i hψ(0)| ,(3.82)

evolution under Ĥef f preserves the purity of the state and therefore is equivalent to a Schrödinger
type evolution,d1|ψ(ti )i = Ĥef f |ψ(ti )i .
dti
ı~(3.83)

Also, the application of a jump operator to a pure state preserves the purity of the state,(k)LˆJ |ψ(ti )i hψ(ti )| = Ŝk |ψ(ti )i hψ(ti )| Ŝk† ⇐⇒ Ŝk |ψ(ti )i .

(3.84)

36
Therefore, the trajectory can be rewritten in terms of state vectors asıı|ψt1 ,k1 ;t2 ,k2 ;··· ;tn ,kn (t)i = e− ~ Ĥef f (t−tn ) Ŝkn e− ~ Ĥef f (tn −tn−1 ) Ŝkn−1 · · ·ıı×Ŝk2 e− ~ Ĥef f (t2 −t1 ) Ŝk1 e− ~ Ĥef f t1 |ψ(0)i= Ût,tn Ŝkn Ûtn ,tn−1 Ŝkn−1 · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i ,(3.85)

whereıÛti ,ti−1 ≡ e− ~ Ĥef f (ti −ti−1 ) .

(3.86)

Eq. 3.85 says that a quantum trajectory is constructed by evolving an initial state |ψ(0)i with Ĥef f
from 0 to t1 , and then applying the jump operator Ŝk1 , evolving this resulting state from t1 to t2 ,applying the jump operator Ŝk2 , and repeating this procedure n times.

3.4.3

Unraveling the Master Equation

After taking the rest of the trace of Eq. 3.76, it can be written as,
1 =†Ût,0 |ψ(0)i +hψ(0)| Ût,0

+

X Z t
k1 ,··· ,kn

0

Z t
dt1
0

Z t3
dtn · · ·X†
Û Ŝ Û |ψ(0)i + · · ·hψ(0)| Ût†1 ,0 Ŝk†1 Ût,t
1 t,t1 k1 t1 ,0

k1

Z t2
dt2

0

dt1
0

†
Û Ŝ · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i + · · · (3.87)
× hψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ût,tn t,tn kn

The first term in Eq. 3.87 can be identified as the probability that no jumps will occur inthe interval [0, t]. The second term can be identified as the probability that exactly one jump will
occur in the interval [0, t], and likewise, the (n + 1)th term can be identified as the probability thatn jumps will occur in the interval [0, t].
In the (n + 1)th term, the sum over ki is a sum over all of the possible jump operators that
could be applied in the trajectory, Ût,ti−1 Ŝki−1 · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i. Furthermore, the integral
over dti is an integral over all the possible times that the ith jump operator is applied. Therefore,
the quantity,37†
P(t1 ,k1 ;t2 ,k2 ;··· ;tn ,kn ) ≡ hψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ût,t
Û Ŝ · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i (3.88)
n t,tn kn

can be interpreted as the probability of the specific jumps {Ŝk1 , Ŝk2 , · · · , Ŝkn } occurring at the
specific times {t1 , t2 , · · · , tn } in the interval [0, t].
Because of this, Eq. 3.85 is rewritten asq†· · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)ihψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ût,t|ψt1 ,k1 ;t2 ,k2 ;··· ;tn ,kn (t)i =nÛt,tn Ŝkn Ûtn ,tn−1 Ŝkn−1 · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)iq††hψ(0)| Ût1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ût,t
Û Ŝ · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i
n t,tn kn
q|ψt1 ,k1 ;t2 ,k2 ;··· ;tn ,kn (t)i=, (3.89)
P(t1 ,k1 ;t2 ,k2 ;··· ;tn ,kn ) phψt1 ,k1 ;t2 ,k2 ;··· ;tn ,kn (t)|ψt1 ,k1 ;t2 ,k2 ;··· ;tn ,kn (t)i×

Thus, we know the probability with which every single possible trajectory will come into the Dyson
expansion of the master equation.
The probability of getting a specific trajectory can be split up into a product of the probabilities of intermediate steps. This can be seen by rewriting it asP(t1 ,k1 ;t2 ,k2 ;··· ;tn ,kn ) =†Û Ŝ · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)ihψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ût,tn t,tn kn

hψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ŝkn · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i×hψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n Ŝkn · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)ihψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 · · · Ŝk†n−1 Ŝkn−1 · · · Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i×···×hψ(0)| Ût†1 ,0 Ŝk†1 Ût†2 ,t1 Ŝk†2 Ŝk2 Ût2 ,t1 Ŝk1 Ût1 ,0 |ψ(0)i×hψ(0)| Ût†1 ,0 Ŝk†1 Ŝk1 Ût1 ,0 |ψ(0)ihψ(0)| Ût†1 ,0 Ŝk†1 Ŝk1 Ût1 ,0 |ψ(0)ihψ(0)|ψ(0)i(3.90)

Therefore, a trajectory with the record (t1 , k1 ; t2 , k2 ; · · · ; tn , kn ) can be constructed in a piecewise
fashion, one time evolution and jump at time. It is beneficial to construct a trajectory this way,since it is fairly easy to obtain the probability distribution for getting the record {t1 , k1 ; · · · ; ti , ki }conditioned upon having gotten the record {t1 , k1 ; · · · ; ti−1 , ki−1 }.
To find this, it is noted that for a small time interval [0, τ ], the probability that an emission
occurs at time τ should be related to the probability that no emission occurs on the interval. The

38equation of motion for the norm of |ψ(τ )i isdhψ(τ )|ψ(τ )i =
dt
=
=
=
=

d hψ(τ )|d |ψ(τ )i|ψ(τ )i + hψ(τ )|
dt 
dt



1
1 †
Ĥef f |ψ(τ )ihψ(τ )| − Ĥef f |ψ(τ )i + hψ(τ )|
ı~
ı~

 

1 †
1
hψ(τ )| − Ĥef f +Ĥef f |ψ(τ )i
ı~
ı~
!
!
1
ı~ X †
1
ı~ X †hψ(τ )| −Ĥ0 +Ŝk Ŝk +Ĥ0 −Ŝk Ŝk |ψ(τ )i
ı~
2
ı~
2
k
k
X
−
hψ(τ )| Ŝk† Ŝk |ψ(τ )i .
(3.91)
kIntegrating both sides from t0 = 0 to t0 = τ ,Z τhψ(τ )|ψ(τ )i − hψ(0)|ψ(0)i = −

dt0

0

X

ψ(t0 ) Ŝk† Ŝk ψ(t0 ) ,(3.92)

kor since hψ(0)|ψ(0)i = 1,
Z τ

dt0

X

0

ψ(t0 ) Ŝk† Ŝk ψ(t0 ) + hψ(τ )|ψ(τ )i = 1.

(3.93)k

The above equation says that the loss of norm of |ψ(τ )i, produced by evolving with Ĥef f , is takennoup by the set of states Ŝk |ψ(τ )i , and that the sum of the norm of both is equal to one.

3.4.4

The Quantum Jumps Algorithm

We are now in a position to describe the quantum Monte Carlo algorithm that simulates the
Dyson expansion of the master equation. This algorithm can be performed as follows:(1) Start in the state |ψi i at time ti .
(2) Pick a random number r1 in the interval [0, 1].
(3) Evolve |ψi i under Ĥef f ,ı|ψ(t)i = e− ~ Ĥef f (t−ti ) |ψi i ,until hψ(t)|ψ(t)i = r1 , and call this time ti+1 .
(4) Choose another random number r2 on the interval [0, 1].

(3.94)

39(5) Set r2 equal to the cumulative distribution function of the conditional probabilities of adecay into mode k at time ti+1 ,
Pk=ki
r2 =

k=0

hψ(ti+1 )| Ŝk† Ŝk |ψ(ti+1 )i†k hψ(ti+1 )| Ŝk Ŝk |ψ(ti+1 )iP,(3.95)and solve for ki , the ith random mode of decay.
(6) Apply the jump operator corresponding to ki to the state |ψ(ti+1 )i, normalize, and call thisnew state|ψi+1 i =Ŝki |ψ(ti+1 )ihψ(ti+1 )| Ŝk†i Ŝki |ψ(ti+1 )i.

(3.96)(7) Repeat steps 1-6 until t ≥ T, the time at which the experiment ends, keeping track of thepieces of the trajectory |ψi i and construct a trajectory out of these pieces, |ψ1 i 0 ≤ t < t1





 |ψ2 i t1 ≤ t < t2
|ψR i =
..



.





 |ψ i tnn−1 ≤ t < tn ,(3.97)

where R = {t1 , k1 ; t2 , k2 ; · · · ; tn , kn }, is the record of the random times and modes of decayof the trajectory.
(8) The simulation the Dyson expansion of the master equation, and hence the full time evolution of the density operator under the master equation, is given by repeating steps 1-7
N times, and forming,N

1 X |ψR im hψR |m,
N →∞ NhψR |m |ψR im

ρ̂(t) = lim(3.98)

m=1

where |ψR im is the mth trajectory constructed by repeating steps 1-6 N times. The expectation value of any observable Ô can be found by,Nˆ =tr[Ôρ(t)]1 X hψR |m Ô |ψR imlim.
N →∞ NhψR |m |ψR im
m=1(3.99)

40
3.4.5

Computer Simulation of Quantum Jumps

Although analytical solutions exist for simple systems, such as the damped harmonic oscillator, for more complex systems, for instance where multiple atoms are coupled to a cavity field,
analytic solutions do not exist. To treat these more complex systems, the numerical solution methods must be relied on. Even so, it will be advantageous to demonstrate these methods for thesimple systems since they are easily extended to more complicated systems.
In Fig. 3.1, the results of a computer simulation which implements the quantum jumps algorithm for the damped harmonic oscillator is shown. We consider the average number of excitationâ† (t)â(t) in the damped mode with n̄ = 1 as a function of t/κ . For reference, the analyticalsolution for â† (t)â(t) , Eq. 3.45, can be seen in black. In a single trajectory, which is shown inblue, there is a precise number of excitations until a random jump occurs at a random time. Then,the number of excitations is increased or decreased by one, depending on which jump operator was
randomly chosen. When many of these random trajectories are averaged together, the correct evolution for â† (t)â(t) is approached, which is demonstrated in red, where 1000 random trajectories
are averaged over.

3.5

Quantum State Diffusion
The unraveling of the master equation into quantum jumps is only one of many possible

unravelings. It corresponds to a measurement scheme in which light from the system is directly
detected on a photodetector. Another possible measurement scheme is homodyne detection, in
which the light coming from the system is mixed with a strong laser called the local oscillator.
An unraveling of the master equation that corresponds to homodyne detection yields the quantum
state diffusion equation. It is basically the limit of the quantum jump method in which the mixed
light from the system and local oscillator has an infinite number of infinitesimally small jumps.

41<a †a>
4

3

2

1

1

2

3

4

5

kt

Figure 3.1: The quantum jumps algorithm is used to find â† (t)â(t) for a damped harmonic
oscillator with n̄ = 1, and α(0) = 2. The analytic expression for â† (t)â(t) is shown in black, asingle trajectory is shown in blue, and 1000 trajectories is shown in red.

42
3.5.1

Master Equation of a Homodyne System

Consider the master equation for a damped single mode cavity field of linewidth κ,odκn †ρ̂(t) = −â âρ̂(t) + ρ̂(t)â† â − 2âρ̂(t)â† .
dt
2(3.100)

Where â is the annihilation operator of the mode, and ρ̂(t) is the system’s density operator. If thisdamped cavity field is mixed with strong local oscillator of amplitude β through a beam splitterwith high transmittivity, the amplitude of the resulting field that is incident on a detector is given
by,
Ĉ =

√

κ (â + β) ,(3.101)

where we have assumed that the phase of the local oscillator is zero. In an interaction picture thatrotates at the frequency of the cavity field, the master equation describing this system is of theform of Eq. 3.63 with Sˆk → Ĉ and Ĥ0 dropped,od1n †Ĉ Ĉ ρ̂(t) + ρ̂(t)Ĉ † Ĉ − 2Ĉ ρ̂(t)Ĉ †ρ̂(t) = −dt2
= Lˆ0 ρ̂(t) + LˆJ ρ̂(t),(3.102)where,
ı
†Ĥef f ρ̂(t) − ρ̂(t)Ĥeff ,~(3.103)ı~ †Ĉ Ĉ,
2(3.104)

LˆJ ρ̂(t) = Ĉ ρ̂(t)Ĉ † .

(3.105)L̂0 ρ̂(t) = −withĤef f = −andThe Dyson expansion of Eq. 3.102 has the same form as Eq. 3.76, so that the analysis of theprevious section carries over directly into this section.
In a time interval ∆t, the photon flux that reaches the photodetector is
E
D
ED
E
D
Ĉ † Ĉ ∆t = κ∆t â† + β â + β DE≈ κ∆t β â† + â + β 2 ,(3.106)

43
where â† â has been dropped since β is large. We assume κ∆t ∼ β −3/2 , so that for largeβ, the number of photodetections n = κ∆tβ 2 ∼ β 1/2 is large, while the change of the systemκ∆tβ â† + â ∼ β −1/2 is small.
By Eq. 3.87, the probability that n jumps occur in a time interval ∆t, is given byZ ∆t
dtn · · ·Pn (∆t) =

Z t2

Z t3

0

dt1

dt2
0

0

†
hψ(0)| Ût†1 ,0 Ĉ † Ût†2 ,t1 Ĉ † · · · Ĉ † Û∆t,tÛ
Ĉ · · · Ĉ Ût2 ,t1 Ĉ Ût1 ,0 |ψ(0)i ,n ∆t,tn(3.107)

where,
Ûti ,ti−1ı= e− ~ Ĥef f (ti −ti−1 )1†= e− 2 Ĉ Ĉ(ti −ti−1 ) .

(3.108)To order O(β −3/2 ), Ûti ,ti−1 and Ĉ commute, so that,Pn (∆t) =

n

tn
†
Ĉ † Ĉ Ût,0 |ψ(0)i .
hψ(0)| Û∆t,0n!

(3.109)After substitution of Eqns. 3.101 and 3.108, Eq. 3.109 becomes,
2

Pn (∆t) = e−(κ∆tβ )(κ∆tβ 2 )n 
1 + (n − (κ∆tβ 2 ))2xβ −1 + O(β −3/2 ) ,n!

(3.110)where x = hψ(0)| (â† + â)/2 |ψ(0)i. From Eq. 3.110, it can be shown [87] that the random variable ncorresponding to this probability distribution can be approximated by a Gaussian random variableof the form,



 √
n = κ∆tβ 2 1 + 2xβ −1 + O(β −3/2 ) + κ∆W β 1 + O(β −1/2 ) ,(3.111)

where ∆W is a Wiener increment [34] which satisfies h∆W i2avg = ∆t.
The state of a specific trajectory after a time interval ∆t, is given by Eq. 3.85 which now
takes the form,|ψR (∆t)i = Û∆t,tn Ĉ Ûtn ,tn−1 Ĉ · · · Ĉ Ût2 ,t1 Ĉ Ût1 ,0 |ψ(0)i ,(3.112)where the subscript R denotes a specific record of jump at times R = {t1 , t2 , · · · , tn }. To order44O(β −3/2 ), Ûti ,ti−1 and Ĉ commute, so that Eq. 3.112 can be rewritten,|ψR (∆t)i = Û∆t,0 Ĉ n |ψ(0)iκ†= e− 2 (â +β )(â+β)∆t√nκ(â + β) |ψ(0)i .

(3.113)

Finally, Eq. 3.113 is expanded to O(β −3/2 ), Eq. 3.111 is substituted for n, and the limit of β → 0 is
taken, resulting in the quantum state diffusion equation for a trajectory corresponding to homodyne
measurement,|ψR (t + ∆t)i =




√
1 †
1 + κ∆t − â â + 2x(t)â + κ∆W (t)â |ψR (t)i .
2(3.114)

In numerical simulations, Eq. 3.114 must be normalized after each time step. The density operator
is then reconstructed by simulating N random trajectories, and forming,N

1 X |ψR im hψR |m,
N →∞ NhψR |m |ψR im

ρ̂(t) = lim(3.115)m=1

The expectation value of any observable Ô can be found by,N

X hψR | Ô |ψR imm
ˆ = lim 1
tr[Ôρ(t)].
N →∞ NhψR |m |ψR im(3.116)

m=1

3.5.2

Computer Simulation of the Quantum State Diffusion Equation

In Fig. 3.2, the results of a computer simulation which solves the quantum state diffusion
equation for random trajectories, for the damped harmonic oscillator is shown. We consider the
average number of excitation â† (t)â(t) in the damped mode with n̄ = 1 as a function of t/κ . For
reference, the analytical solution for â† (t)â(t) , Eq. 3.45, can be seen in black. In a single trajectory, which is shown in blue, there are many very small jumps, so that the number of excitations
randomly diffuses in time. When many of these random trajectories are averaged together, thecorrect evolution for â† (t)â(t) is approached, which is demonstrated in red, where 1000 random
trajectories are averaged over.

45<a †a>
4

3

2

1

1

2

3

4

5

kt

Figure 3.2: The quantum state diffusion equation is solved for random trajectories which are used tofind â† (t)â(t) for a damped harmonic oscillator with n̄ = 1, and α(0) = 2. The analytic expression
for â† (t)â(t) is shown in black, a single trajectory is shown in blue, and 1000 trajectories is shownin red.

46

3.6

The Master Equation for a General Cavity QED System
In a cavity QED system, in which energy is re-supplied to the atoms, energy can not be

supplied directly at the frequency of the considered transition. Otherwise, when a large numberof atoms become excited, stimulated emission would occur, and the inversion would never become
positive. However, if a higher third level is repumped, and this level spontaneously decays down,
the system can achieve positive inversion. This repumping process can be thought of as reverse
spontaneous emission, and can be modeled with the master equation for a two level atom withγ → w, σ̂ + → σ̂ − , and σ̂ − → σ̂ + , where w is now the repumping rate.
In addition to the damping of the excited state population due to spontaneous emission, it
is also possible that decoherence of the phase between excited and ground states can occur. Thisdephasing can occur, for instance, due to elastic collisions in an atomic cloud. It can be modeledby a system reservoir coupling Hamiltonian ĤSR of the form,
ĤSR = ~Xλj σ̂ z r̂j + λ∗j σ̂ z r̂j† ,(3.117)jwhere r̂j is the annihilation operator of the jth mode in the reservoir Hilbert space, λj is coupling
rate between the system and the jth mode of the reservoir, and σ̂ z is a Pauli operator, which obeys
Eqns. 3.58. This form of ĤSR models the dephasing between excited and ground states, since theaction of the operator σ̂ z is to change the relative phase between the excited and ground state byπ. Eq. 3.117 is the same as Eq. 3.57 with σ̂ + → σ̂ z . Therefore the derivation is preformed in thesame way as in Sec. 3.3, with γ → 1/2T2 , where 1/T2 is the dephasing rate. Then, the masterequation describing atomic dephasing is given by,


d˜
1
z z
z
˜ + ρ̂(t)σ̂˜˜ρ̂(t) = −(n̄ + 1) σ̂ z σ̂ z ρ̂(t)σ̂ − 2σ̂ z ρ̂(t)σ̂dt
4T2
1  z z˜
z z
z
˜
˜
−
n̄ σ̂ σ̂ ρ̂(t) + ρ̂(t)σ̂σ̂ − 2σ̂ z ρ̂(t)σ̂
4T2

 1 

1
z
z
˜
˜
˜
˜
=
(n̄ + 1) σ̂ z ρ̂(t)σ̂− ρ̂(t)n̄ σ̂ z ρ̂(t)σ̂− ρ̂(t).
2T2
2T2(3.118)Many cavity QED experiments are done at close to zero temperature, so that average numberof thermal excitations n̄ = 0. Putting it all together, the master equation for N two level atoms

47
of frequency ωa , decay rate γ, repumping rate w, and dephasing rate 2T1 2 , coupled with rate g to asingle mode cavity of decay rate κ is given by,
N n
i κn
o γX
o
d
ıh
†
†
†
ρ̂ = − Ĥ, ρ̂ −â âρ̂ + ρ̂â â − 2âρ̂â −σ̂j+ σ̂j− ρ̂ + ρ̂σ̂j+ σ̂j− − 2σ̂j− ρ̂σ̂j+
dt
~
2
2
j=1

w
−
2

N n
X

σ̂j− σ̂j+ ρ̂ + ρ̂σ̂j− σ̂j+ − 2σ̂j+ ρ̂σ̂j−

j=1

o

N

1 X z zσ̂j ρ̂σ̂j − ρ̂ ,+2T2

(3.119)

j=1

with,
Ĥ =

N
N

~ωa X z~g X  − †σ̂j + ~ωc â† â +σ̂j â + σ̂j+ â ,22
j=1

(3.120)j=1

where ρ̂ is the system density operator, and â is the annihilation operator for the cavity mode. The
jth atom is described by the operators,σ̂j− = |gij he|jσ̂j+ = |eij hg|jσ̂jz = |eij he|j − |gij hg|j= |eij he|j + |gij hg|j ,1̂j(3.121)

where |eij is the excited state, |gij is the ground state, and 1̂j is the identity operator for thejth atom. Then we can use the fact that |eij and |gij are orthonormal, i.e. he|i |gii = 0, andhe|i |eii = hg|i |gii = 1, to simplify calculations. The commutation relations between the atomic
operators are,hi
σ̂i+ , σ̂j− = σ̂iz δijhi
σ̂iz , σ̂j+ = 2σ̂i+ δijhi
σ̂i− , σ̂j− = 2σ̂i− δij ,(3.122)where δij is the Kronecker delta.

3.7

Cumulant Expansion to Second OrderThe quantum jumps and quantum state diffusion methods are only able to describe verysmall systems, with atom numbers N ∼ 10. This is because the basis used to describe the system

48
scales with atom number N as 2N . Many experiments done today have atom numbers N ∼105 or higher, and in order to describe these larger systems, approximation methods must be
employed. The cumulant expansion is one such approximate method that allows large systems
to be described. In this method, equations for expectation values of observables of interest will
be derived. These equations will couple to higher order moments, which will in turn couple toeven higher order moments. We choose to factorize all moments higher that second order, using
the cumulant expansion [56]. This will allow a closed set of equations for expectation values tobe derived. With the help of the quantum regression theorem, and again using the cumulant
expansion, equations for two-time expectation values can also be derived.

3.7.1

Closed Set of One-Time Expectation Values

We begin by deriving the equation of motion for hσ̂1z i, the inversion for a single atom of theatomic ensemble. Using Eq. 3.119,d zdhσ̂1 i = tr[σ̂1z ρ̂]
dt
dt
N
N
h 



i
 ıg X
ıωa X  z ztr σ̂1 σ̂j ρ̂ − σ̂jz σ̂1z ρ̂ −tr σ̂1z σ̂j− â† + σ̂j+ â ρ̂ − σ̂j− â† + σ̂j+ â σ̂1z ρ̂
= −
2
2
j=1

−

j=1

N
i
γX h z + −tr σ̂1 σ̂j σ̂j ρ̂ + σ̂j+ σ̂j− σ̂1z ρ̂ − 2σ̂j+ σ̂1z σ̂j− ρ̂
2
j=1

−

N
i
wX h z − +tr σ̂1 σ̂j σ̂j ρ̂ + σ̂j− σ̂j+ σ̂1z ρ̂ − 2σ̂j− σ̂1z σ̂j+ ρ̂
2
j=1

+

N

1 X  z z ztr σ̂j σ̂1 σ̂j ρ̂ − σ̂1z ρ̂ ,2T2
j=1

(3.123)

49
where the cyclic property of the trace has been used. After using Eq. 3.121 to simplify, Eq. 3.123
becomes,id z
ıg hhσ̂1 i = − tr −σ̂j− â† + σ̂j+ âρ̂ − σ̂j− â† + σ̂j+ âdt2γ  + −− tr σ̂1 σ̂1 ρ̂ + σ̂1+ σ̂1− ρ̂ + 2σ̂1+ σ̂1− ρ̂2w  − +− tr −σ̂1 σ̂1 ρ̂ − σ̂1− σ̂1+ ρ̂ − 2σ̂1− σ̂1+ ρ̂
2DE= ıg â† σ̂1− − âσ̂1+ − (w + γ) (hσ̂1z i − d0 ) ,w−γ. Since â† σ̂1−where d0 ≡ w+γ(3.124)

couples into Eq. 3.124, its equation must also be calculated.

Eq. 3.119 is used again to find,dd D † −E
= tr[â† σ̂1− ρ̂]â σ̂1dt
dt
Ni ıω hiıωa X h † − zc= −tr â σ̂1 σ̂j ρ̂ − σ̂jz â† σ̂1− ρ̂ −tr â† σ̂1− â† âρ̂ − â† ââ† σ̂1− ρ̂
2
2
j=1

−

ıg
2

N
X

i




h
tr â† σ̂1− σ̂j− â† + σ̂j+ â ρ̂ − σ̂j− â† + σ̂j+ â â† σ̂1− ρ̂

j=1

i
κ h
− tr â† σ̂1− â† âρ̂ + â† ââ† σ̂1− ρ̂ − 2â† â† σ̂1− âρ̂2NiγX h † − + −−tr â σ̂1 σ̂j σ̂j ρ̂ + σ̂j+ σ̂j− â† σ̂1− ρ̂ − 2σ̂j+ â† σ̂1− σ̂j− ρ̂
2
j=1

−

N
i
wX h † − − +tr â σ̂1 σ̂j σ̂j ρ̂ + σ̂j− σ̂j+ â† σ̂1− ρ̂ − 2σ̂j− â† σ̂1− σ̂j+ ρ̂
2
j=1

+

N
i
1 X h z † − ztr σ̂j â σ̂1 σ̂j ρ̂ − â† σ̂1− ρ̂ .
2T2

(3.125)j=1

With the use of the commutation relations for creation and annihilation operators, Eq. 3.43, andwith the help of Eqns. 3.121, Eq. 3.125 becomes,ihid D † −Eıωa h † −tr â σ̂1 ρ̂ + â† σ̂1− ρ̂ + ıωc tr â† σ̂1− ρ̂â σ̂1= −dt
2
Niıg X h † − +−tr â âσ̂1 σ̂j ρ̂ − ââ† σ̂j+ σ̂1− ρ̂2
j=1iκ h † − i γ h † − i w h † − i
1 h− tr â σ̂1 ρ̂ − tr â σ̂1 ρ̂ − tr â σ̂1 ρ̂ − tr â† σ̂1− ρ̂ .
2
2
2
T2

(3.126)

50
The coupling term requires closer attention, since here, the j 6= 1 terms in the sum are not zero:d D † −E
1 D † −Eκ γ w
â σ̂1â σ̂1=−ıωa + ıωc − − − −
dt
2
2
2
T2
N
i ıg Xhiıg h † − +† + −tr â† âσ̂1− σ̂j+ ρ̂ − ââ† σ̂j+ σ̂1− ρ̂− tr â âσ̂1 σ̂1 ρ̂ − ââ σ̂1 σ̂1 ρ̂ −
2
2
j6=1
D

E
1
κ γ w
â† σ̂1−=−ıδ − − − −
2
2
2
T2


ıg
ıg D † z E 1(3.127)â âσ̂1 + (1 + hσ̂1z i) + (N − 1) σ̂1+ σ̂2− ,+222where δ ≡ ωa − ωc , and the identical atom assumption has been invoked in the last term. It isclear now that it is the coupling term that causes nth order moments to couple to (n + 1)th ordermoments. We therefore use the cumulant expansion [56], dropping cumulants higher than second
order. Applying this to â† âσ̂1z yields,Dâ† âσ̂1z

E

D E
E
D
E
D E
â† â hσ̂1z i + â† σ̂1z hâi + hâσ̂1z i â† − 2 â† hâi hσ̂1z i
D E=â† â hσ̂1z i ,=D(3.128)where the last line follows from the fact that hâi = 0. Eq. 3.128 is substituted into Eq. 3.127 toobtain,d D † −E=â σ̂1
dt




κ γ w
1 D † −E−ıωa + ıωc − − − −â σ̂1
2
2
2
T2
D E

ıg
1
ıg
+
â† â hσ̂1z i + (1 + hσ̂1z i) + (N − 1) σ̂1+ σ̂2− .
2
2
2(3.129)The equations of motion of both â† â and σ̂1+ σ̂2− are now calculated, and the cumulant
expansion is used to expand moments higher than second order, dropping all cumulants of higher
than second order. Then, we arrive at the closed set of equations,
D
E

d z
hσ̂1 i = ıg â† σ̂1− − âσ̂1+ − (w + γ) (hσ̂1z i − d0 )dt
1 D † −Ed D † −Eκ γ w
â σ̂1=−ıδ − − − −â σ̂1
dt
2
2
2
T2
D E

1
ıg
ıg
†
z
z
+
â â hσ̂1 i + (1 + hσ̂1 i) + (N − 1) σ̂1+ σ̂2−
2
2
2


D
E

d
2
ıg
σ̂1+ σ̂2− = − w + γ +σ̂1+ σ̂2− − hσ̂1z i â† σ̂1− − âσ̂1+
dt
T2
2
D E ıN g DEd D † Eâ â= −κ â† â −â† σ̂1− − âσ̂1+ ,dt
2(3.130)51
w−γwith d0 ≡ w+γ. Eqns. 3.130 can be solved on a computer with a numerical differential equation

solver in fractions of a second. Even so, for most purposes we are only interested in the steadystate solution, which can be found by setting the time derivatives to zero, and solving the resulting
algebraic equations. These algebraic equations can also be easily solved on a computer, which iswhat is typically done in practice. However, for a broad strokes understanding of the physics, simple
analytic solutions are desirable, and can be obtained in the limit of large N . We set δ = T12 = γ = 0,
for simplicity, and combine Eqns. 3.130 with the time derivatives set to zero,DE
0 = ıg â† σ̂1− − âσ̂1+ − w (hσ̂1z i − d0 ) , κ wDE ıg D E
1
ıg
0 =
− −
â† σ̂1− +â† â hσ̂1z i + (1 + hσ̂1z i) + (N − 1) σ̂1+ σ̂2−
2
2
2
2
2

E
D
ıg
0 = −w σ̂1+ σ̂2− − hσ̂1z i â† σ̂1− − âσ̂1+
2
D E ıN g DE
(3.131)0 = −κ â† â −â† σ̂1− − âσ̂1+ .
2
Then, after defining,
C≡

N g2
w+κ

≡

κ
,
N (w + κ)(3.132)

we arrive at the quadratic,
0 = C (1 − ) (1 − hσ̂1z i)2 + (1 − C + 2C) (1 − hσ̂1z i) − 2C.

(3.133)For  = 0, (1 − hσ̂1z i) = 1 − C1 , so for   1, we add the quantity ξ to this,(1 − hσ̂1z i) = 1 −1+ξC(3.134)

where ξ is of order . Eq. 3.134 is substituted into Eq. 3.133, and only first order terms in  are
kept. In this way, ξ is found to be,ξ= C2 + 1,
C C −1(3.135)

From Eqns. 3.131, the rest of the steady state observables are easily solved for. All together, they

52
are,
hσ̂1z i =σ̂1+ σ̂2− =D Eâ† â
=DEâ† σ̂1−=

1
 (C 2 + 1)−
C
C (C − 1)1 zhσ̂ i (1 − hσ̂1z i)
2 1Nw(1 − hσ̂1z i)
2κiw(1 − hσ̂1z i) .
2g

(3.136)The optimum repumping value w = wopt , at which the maximum photon number is achieved,dis found by setting dwâ† â = 0 and solving for the w at which this occurs. Doing this, it is found,

0=

wopt
N
d D † Eâ â =
− 2
dw
2κ
g

(3.137)so thatwopt =

N g2.
2κ

(3.138)Substituting Eq. 3.138 into Eqns. 3.136 yields,
hσ̂1z iopt =D Eâ† â
=
opt

3.7.2

1
2
N 2g2
8κ2

1
σ̂1+ σ̂2− opt =
8
E
D
iN g
† −
=
â σ̂1.
8κ
opt(3.139)Closed Set of Two-Time Expectation Values

The cumulant expansion method can also be used to calculate â† (t + τ )â(t) using thequantum regression theorem demonstrated in Sec. 3.2.3. With Â → â† and B̂ → â in Eq. 3.49, andiinoting now that Ŝ(τ ) = e− ~ Ĥτ âρ̂e ~ Ĥτ , we arrive at,

 
E
d D †
d
â (t + τ )â(t)= trŜ(τ ) â†
dτ
dτ
h
i
= −ıωc tr â† â† âŜ(τ ) − â† ââ† Ŝ(τ )
N

i


ıg X h †  − †tr â σ̂j â + σ̂j+ â Ŝ(τ ) − σ̂j− â† + σ̂j+ â â† Ŝ(τ )2j=1iκ h † †− tr â â âρ̂ + â† ââ† ρ̂ − 2â† â† âρ̂ ,(3.140)2−

53
where the cyclic property of the trace has been used. The commutation relations for creation andannihilation operators, Eq. 3.43, are used to simplify,Ni κ hi ıg X
E
h
h
i
d D †
†
†
â (t + τ )â(t)= ıωc tr â Ŝ(τ ) − tr â Ŝ(τ ) +tr σ̂j+ Ŝ(τ )
dτ
2
2
j=1


D
E
D
E
κ
ıgN
=
ıωc −
â† (t + τ )â(t) +σ̂j+ (t + τ )â(t) ,22(3.141)

where we have assumed identical atoms in the last line. Since σ̂1+ (t + τ )â(t) couples into Eq. 3.141,
its equation of motion must also be found. We therefore calculate,dσ̂ + (t + τ )â(t)
dτ 1

 
d
Ŝ(τ ) σ̂1+
= tr
dτ
N
i
ıωa X h + z= −tr σ̂1 σ̂j Ŝ(τ ) − σ̂jz σ̂1+ Ŝ(τ )
2


j=1

−

−

ıg
2

N
X

i



h 
tr σ̂1+ σ̂j− â† + σ̂j+ â Ŝ(τ ) − σ̂j− â† + σ̂j+ â σ̂1+ Ŝ(τ )

j=1

N
i
γX h + + −tr σ̂1 σ̂j σ̂j Ŝ(τ ) + σ̂j+ σ̂j− σ̂1+ Ŝ(τ ) − 2σ̂j+ σ̂1+ σ̂j− Ŝ(τ )
2
j=1

−

N
i
wX h + − +tr σ̂1 σ̂j σ̂j Ŝ(τ ) + σ̂j− σ̂j+ σ̂1+ Ŝ(τ ) − 2σ̂j− σ̂1+ σ̂j+ Ŝ(τ )
2
j=1

+

N
i
1 X h z + z
tr σ̂j σ̂1 σ̂j Ŝ(τ ) − σ̂1+ Ŝ(τ ) ,2T2(3.142)j=1

where the cyclic property of the trace has been used. After simplifying using Eq. 3.121, Eq. 3.142
becomes,dσ̂ + (t + τ )â(t)dτ 1hi ıg hi= ıωa tr σ̂1+ Ŝ(τ ) − tr â† σ̂1+ σ̂1− − σ̂1− σ̂1+ Ŝ(τ )
2
i w h
i
i
γ h +
1 h
− tr σ̂1 Ŝ(τ ) − tr σ̂1+ σ̂j− σ̂j+ Ŝ(τ ) − tr σ̂1+ Ŝ(τ )
2
2
T2

D
E
D
E
γ w
1
ıg † z=ıωa − − −σ̂1+ Ŝ(τ ) −â σ̂1 Ŝ(τ )
2
2
T2
2


E
γ w
1
ıg D †=
ıωa − − −σ̂1+ (t + τ )â(t) −â (t + τ )σ̂1z (t + τ )â(t) .
2
2
T2
2(3.143)

54
To be consistent with the approximations of Sec. 3.7.1, the third order moment in Eq. 3.143 isfactorized and the third order cumulant is dropped, which yields,
D

E
D E
D
E
D E
D E
â† (t + τ )σ̂1z (t + τ )â(t)=â† â hσ̂1z i + â† σ̂1z hâi + hâσ̂1z i â† − 2 â† hâi hσ̂1z i
DE= hσ̂1z (t + τ )i â† (t + τ )â(t) .
(3.144)

Putting it all together, in an interaction picture that rotates at ωa , defined by,˜ − (t) = e−iωa t σ̂ − (t),σ̂11˜ = e−iωa t â(t)â(t)(3.145)

we have,EE ıgN DEκ D †d D †â (t + τ )â(t)= − ıδ +â (t + τ )â(t) +σ̂j+ (t + τ )â(t)
dτ
2
2
D
E
d
ıg
σ̂ + (t + τ )â(t) = − hσ̂1z (t + τ )i â† (t + τ )â(t)
dτ 1
2


1
γ w
+ +
σ̂1+ (t + τ )â(t) ,
−
2
2
T2

(3.146)

where δ ≡ ωa − ωc .
To solve Eqns. 3.146, we note that they are of the form,dA(τ ) = c1 A(τ ) + c2 B(τ )dτd
B(τ ) = c3 A(τ ) + c4 B(τ ),dτ(3.147)

which can be solved by taking the Laplace transform of both sides,
Z ∞


Z ∞
d
dτ e
A(τ )
=
dτ e−sτ {c1 A(τ ) + c2 B(τ )}
dτ


Z0 ∞
Z 0∞
d
dτ e−sτB(τ )=dτ e−sτ {c3 A(τ ) + c4 B(τ )} .
dτ
0
0
−sτ



(3.148)After integrating the left hand sides by parts, and defining,Z ∞A(s) =
Z0 ∞
B(s) =
0dτ e−sτ A(τ )dτ e−sτ B(τ ),(3.149)

55
Eq. 3.148 becomes,−A(0) + sA(s) = c1 A(s) + c2 B(s)−B(0) + sB(s) = c3 A(s) + c4 B(s).

(3.150)Eq. 3.150 is algebraically solved for A(s), yielding,A(s) ==(s − c4 ) A(0) + c2 B(0)2s − (c1 + c4 ) s + (c1 c4 − c2 c3 )(s − c4 ) A(0) + c2 B(0),s − ∆ν2 +s − ∆ν2 −(3.151)

where,q(c1 + c4 )2 − 4 (c1 c4 − c2 c3 )q11
(κ + w) ∓(−κ + w)2 − 4N g 2 hσ̂1z i.
22∆ν± = − (c1 + c4 ) ∓=(3.152)

By the method of partial fractions, Eq. 3.151 can be put into the form,
A(s) = Xs − ∆ν2 ++Ys − ∆ν2 −,(3.153)so that taking the inverse Laplace transform yields
D
E
∆ν+ τ
∆ν− τ
â† (t + τ )â(t) = A(τ ) = Xe− 2 + Y e− 2 .

(3.154)Since ∆ν+  ∆ν− , the second term decays fast to zero, so that the first term dominates. Then,since the Fourier transform of an exponential with decay rate ∆ν is a Lorentzian with a full widthhalf maximum linewidth ∆ν, the steady state linewidth ∆ν of a system described by Eq. 3.119 isgiven by,∆ν =11
(κ + w) −22q(−κ + w)2 − 4N g 2 hσ̂1z i.

(3.155)κ, we find a simple
By plugging Eq. 3.136 into Eq. 5.50, and again assuming N  w+κ

expression for the linewidth,∆ν =C2 + 1Ω2 κ.
C(C − 1) (w + κ)2(3.156)

56

3.8

Conclusion
We have seen that an open quantum system can be described by coupling it to a large

reservoir, and evolving the combined system and reservoir to second order in perturbation theory.
Then, after taking the partial trace over the reservoir Hilbert space, and making the Markov approximation, the master equation was obtained. For simple systems, such as the damped harmonic
oscillator, the master equation can be used to derive exact analytical solutions for expectation
values. With the help of the quantum regression theorem, analytical solutions for two-time expectation values can also be derived. We also introduced the quantum jumps and quantum state
diffusion methods which can be used to solve the master equation using a computer. The quantum jumps method unravels the master equation into random trajectories, which correspond to
the random measurement record of photons leaving the system and being detected. The quantum
state diffusion method unravels the master equation into random trajectories, which correspondto a randomly fluctuating homodyne photocurrent that results from mixing the photons from thesystem with a strong local oscillator. In both cases, the random trajectories are averaged overto obtain the correct evolution in the ensemble sense. Computer simulations using both of thesemethods were demonstrated in the system of the damped harmonic oscillator. For systems with alarge number of atoms, the quantum jumps and quantum state diffusion methods are not possible,
since the size of the Liouville space grows exponentially with atom number. For very large systems,
approximate methods such as the cumulant expansion method must be employed. We derived aclosed set of equations for expectation values by factorizing all moments higher than second orderthat appear in the equations. With the help of the quantum regression theorem, equations for
two-time expectation values were also derived.

Chapter 4

Langevin Equations

Another way to incorporate dissipation into quantum mechanics is through the quantum
Langevin equations. This formalism is completely equivalent to the master equation approach,but instead of the time dependence being on the state of the system, it is put on the operatorsthemselves. Thus, the Langevin equations are in the Heisenberg picture of quantum mechanics,
while the master equation is in the Schrödinger picture. In analogy to the master equation approach,dissipation occurs in the Langevin approach due to the coupling of the small system to a large
reservoir.
The quantum Langevin equations are operator equations, which are cumbersome and noteasily simulated on a computer. Therefore, for solutions methods that do not make approximations,the quantum master equation is a better starting point. However, the quantum Langevin equations
can be used to derive approximate equations that can describe large systems, where exact solution
methods are not possible. Therefore, after the quantum Langevin equations are derived fromfirst principles, the Fokker-Planck equation for a quazi-probability distribution which describes thequantum system, is derived. With the aid of stochastic calculus, this Fokker-Planck equation will
be shown to be equivalent to a set of c-number Langevin equations, which can be used to describe
large systems, and are easily simulated on a computer.
Once this correspondence between quantum operator Langevin equations and c-number
Langevin equations is demonstrated, a more practical route to obtain the c-number Langevinequations is given.

58

4.1

Derivation of the Quantum Langevin Equations
We will begin by deriving the quantum Langevin equations [67] for N two level atoms coupled

to a single cavity mode. The derivation will focus on the decay of the cavity only, for simplicity.
However, in analogy to the master equation derivation, the incoherent decay and repumping of
the atoms is derived in an almost identical manner as presented here. The Hamiltonian Ĥ of the
combined system and reservoir is given by,
Ĥ =



X
X  †
~ωa ˆ
~g  ˆ †Jz + ~ωc â† â +J− â + Jˆ+ â + ~νk b̂†k b̂k + ~λk b̂k â + b̂k â† ,
2
2
k

(4.1)

k

which describes N two level atoms of frequency ωa coupled with frequency g to a single mode cavityof frequency ωc , which is coupled with rates λk to a large reservoir of other harmonic oscillatorsof frequency νk . The cavity mode is described by the annihilation operator â, and the reservoir is
described by the annihilation operators b̂k , while the N two level atoms are described by,
Jˆz =

N
X

σjz ,

j=1

Jˆ± =

N
X

σj± ,

j=1

where σj− , σj+ , σjz are the Pauli matricies for the jth atom.
In the Heisenberg picture, the equation of motion of â is given by,
â˙ =

i
ih
Ĥ, â~"! #gˆ † X†= iωc â â + J− â +λk b̂k â , â
2
k
h
i
h
iX
gh † i ˆ†= iωc â â, â + i â , â J− + i â† , â
λk b̂k2
kXgˆ= −iωc â − i J− − iλk b̂k ,
2
†

k

(4.2)

59
and the equation of motion of b̂k is given by,iihĤ, b̂k
~"
!
#
X †
X †
= i
νl b̂l b̂l +λl b̂l â , b̂k

˙
b̂k =

l

= i

X

l

h
i
X h † iνl b̂†l b̂l , b̂k + iâ
λl b̂l , b̂kll= −iνk b̂k − iλk â.

(4.3)

As an aside, it is now apparent that if the decay of the atoms was accounted for, it would notchange the form of Eq. 4.2 or Eq. 4.3, but would change the equations associated with the atoms.
Eq. 4.3 is now integrated,−iνk tb̂k (t) = b̂k (0)e

Z t
− iλk

0

dt0 e−iνk (t−t ) â(t0 ),

(4.4)

0

and substituted into Eq. 4.2 to obtain,gâ˙ = −iωc â − i Jˆ− − i2X−iνk tλk b̂k (0)e

Z t
− iλk

dt e


â(t )

X

λk b̂k (0)e−iνk t .

0 −iνk (t−t0 )

0

0

k

X
g
= −iωc â − i Jˆ− −
λ2k
2

Z t

0

dt0 e−iνk (t−t ) â(t0 ) − i

0

k

(4.5)

k

An interaction picture is defined by,˜ = eiωc t â,â˜Jˆ− = eiωc t Jˆ− ,so that in this interaction picture, Eq. 4.5 becomes,X˜˙ = −i g J˜
ˆ− −
â
λ2k
2
k

Z t

0
˜
dt0 e−i(νk −ωc )(t−t ) â(t0 ) + F̂,

(4.6)

0

where,X˜F̂ = −iλk b̂k (0)e−i(νk −ωc )t .
k

(4.7)

60
In the second term on the right hand side of Eq. 4.6, the sum over k is changed to an integral
over ν by including a density of states,
Z t
X
X Z t
0
0
dt0 e−i(νk −ωc )(t−t ) â(t0 )dt0 e−i(νk −ωc )(t−t ) â(t0 ) =
∆kλ2k
λ2k
0

k

0

k

Z B
=

dνD(ν) |λ(ν)|2

Z t

0

0

Z B

Z t

=

dνD(ν) |λ(ν)|20

dt0 e−i(ν−ωc )(t−t ) â(t0 )dτ e−i(ν−ωc )τ â(t − τ ).

(4.8)

0

0

In analogy to the derivation of the master equation, it is assumed that the reservoir has a verylarge bandwidth B, so that the reservoir relaxation rate τ1R will be much larger than the systems
relaxation rate τ1S . Then, for τ << τS , â(t − τ ) will be approximately constant, so that it can beevaluated at τ = 0 and factored outside of the integral. The resulting integral of exponentials willthen destructively interfere when τ >> τR , resulting in it averaging to zero. Therefore the upper
limit of that integral can be extended to τ = ∞ with minimal error. In this Markov approximation,
Eq. 4.8 becomes,
Z ∞
Z ∞
X Z t
2
0 −iνk (t−t0 ) ˜ 02
dt eâ(t ) =dνD(ν) |λ(ν)| â(τ )dτ e−i(ν−ωc )τ
λk
k

0

0

0

Z ∞
dνD(ν) |λ(ν)|2 â(t)dτ e−i(ν−ωc −iη)τ
lim
η→0 0
0



Z ∞
1
2
=
dνD(ν) |λ(ν)| â(t) πδ(ν − ωc ) + iP
. (4.9)ν − ωc0Z ∞=where D(ν) is a density of states, and Eq. 3.35 has been again used. Neglecting the principal partof Eq. 4.9, which leads to a small frequency shift related to the Lamb shift, and substituting it into
Eq. 4.6 results in,˜˜˙ = −i g J˜ˆ− − κ â
˜ + F̂,â2
2(4.10)

κ = 2πD(ωc ) |λ(ωc )|2 .

(4.11)

where,For use in later calculations, using Eq. 4.7, we calculate,
D E
E
X D˜
˜
F̂
= −i
λk b̂k (0) e−i(νk −ωc )tk= 0,(4.12)

61DE
˜
since b̂k (0) = 0, because it is the trace of an annihilation operator. We also calculate,D˜˜F̂ † (t)F̂ (t0 )

E

D
E
0
˜
˜
λ∗k λk0 b̂†k (0)b̂k (0) ei(νk −ωc )t e−i(νk0 −ωc )t

=

XX
k

k0

=

X

|λk |2 n̄k ei(νk −ωc )(t−t )

0

k

Z ∞

2

≈ D(ωc ) |λ(ωc )| n̄(ωc )0dνei(νk −ωc )(t−t )−∞= κn̄(ωc )δ(t − t0 ),(4.13)

where the discrete sum has again been approximated as an integral over a density of states, andthe slowly varying variables have been pulled outside the integral, which is valid under the Markov
approximation described above. Also, the fact that,Z ∞0dνei(νk −ωc )(t−t ) = 2πδ(t − t0 ),(4.14)−∞

has been employed. In the same way, it can be shown that,DE˜ ˜F̂ (t)F̂ † (t0 ) = κ (n̄(ωc ) + 1) δ(t − t0 ),(4.15)

D
E D
E
˜ ˜
˜
˜
F̂ (t)F̂ (t0 ) = F̂ † (t)F̂ † (t0 ) = 0.

(4.16)

andThe quantities defined in Eqns. 4.13, 4.15 and 4.16 are called the diffusion coefficients, and areneeded to characterize the system, as we will see later.

4.2

Adiabatic Elimination and the Superradiance Master Equation
As a brief aside, we are now in a position to derive the superradiance master equation. Thisis done by first integrating Eq. 4.10, obtaining,
Z t

 g

κ
0
˜
˜
dt0 e− 2 (t−t ) −i Jˆ− + F̂
2
0
Z t
Z t
κ
κ
κ
g
0 ˜
0 ˜
−2t
˜
= â(0)e−i
dt0 e− 2 (t−t ) Jˆ− +dt0 e− 2 (t−t ) F̂.
2 0
0

−κ
t
˜
˜
2 +
â(t)
= â(0)e(4.17)κ

Assuming the system is in the steady state, t is large enough so that e− 2 t → 0. This allows
the first term in Eq. 4.17 to be dropped, but not the other terms because they have exponentials

62
that depend on t − t0 which can be small. The second term in Eq. 4.17 is now integrated by parts
multiple times, and under the steady state assumption,
Z t

κ
2˜
0 ˜
dt0 e− 2 (t−t ) Jˆ− = Jˆ− −

κ

0

 2
2
d ˜ˆ
J− + · · ·κdt(4.18)

Eq. 4.18 tells us that if,
1 d ˜ˆ˜J− ,Jˆ− >>
κ dt(4.19)

it is valid to make the adiabatic approximation, keeping only the first order term in Eq. 4.18. Under
this approximation, Eq. 4.17 becomes˜
˜ = −i g J˜ˆ− + Ĝ,â(t)κ(4.20)

where
˜
Ĝ =

Z t

κ
0 ˜
dt0 e− 2 (t−t ) F̂.

(4.21)

0

An upper limit on the adiabatic condition, Eq. 4.19, in terms of the parameters of the problem
˜
can be obtained by finding the equation of motion for Jˆ− ,g ˜ h ˜ˆ ˜ˆ i˜˙Jˆ− = i â
J+ , J−2˜J˜ˆz= ig â g
˜ ˜
˜= ig −i Jˆ− + Ĝ Jˆzκg 2 ˜ˆ ˜ˆ˜=
J− Jz + ig G̃Jˆz ,κ(4.22)˜
substituting this into Eq. 4.19, and noting that the largest value of Jˆz occurs when the Blochvector points to one of the poles, having the value N . Eq. 4.19 then says,
κ >> N Cγ,(4.23)

where
C=

g2
.
κγ

(4.24)The master equation for a system of atoms interacting with a damped cavity mode, which has
been adiabatically eliminated is then given by inserting Eq. 4.20 into Eq. 3.119. Since the master

63
D E
˜
˜
equation is course grain averaged over the timescale of the reservoir, it is valid to set Ĝ → ĜRin Eq. 4.20 before inserting it into Eq. 3.119, and still consistent with the approximations of theED˜is given by,master equation. The expectation value Ĝ(t)
R

D

E
˜
Ĝ(t)

Z t
=

R

D
E
κ
0
˜
dt0 e− 2 (t−t ) F̂ (t0 )

0

= 0.

(4.25)

Therefore, we arrive at the superradiance master equation,od˜Cγ n ˜ˆ ˜ˆ˜ ˜˜ ˜ρ̂ = −
J+ J− ρ̃ + ρ̃Jˆ+ Jˆ− − 2Jˆ− ρ̃Jˆ+dt
2Noγ Xn + −˜ ˜ + −˜ +σ̂j σ̂j ρ̂ + ρ̂σ̂j σ̂j − 2σ̂j− ρ̂σ̂
−
j
2
j=1
N

−

o
w Xn − +˜ ˜ − +˜ − .
σ̂j σ̂j ρ̂ + ρ̂σ̂j σ̂j − 2σ̂j+ ρ̂σ̂j2

(4.26)

j=1

4.3

Fokker-Planck Equations for the P Function
Since Eq. 4.10 is an equation of quantum operators, it cannot simply be used in numericalsimulations. We therefore seek a method that can be used for numerical simulation. To this end, theFokker-Planck equation for a quazi-probability distribution describing the system is derived fromthe master equation. It is then demonstrated that the resulting Fokker-Planck equation can beequivalently described in terms of a set of c-number Langevin equations, which can be easily solved
using numerical simulations. The method is demonstrated for the damped harmonic oscillator.
The master equation for the damped harmonic oscillator is,
h
i γ

 γ 

ρ̂˙ = −ıω0 â† â, ρ̂ + (n̄ + 1) 2âρ̂â† − â† âρ̂ − ρ̂â† â + n̄ 2â† ρ̂â − ââ† ρ̂ − ρ̂ââ†
2
2
h
i γ



†
†
†
†
2âρ̂â − â âρ̂ − ρ̂â â + γ n̄ âρ̂â† + â† ρ̂â − â† âρ̂ − ρ̂ââ† . (4.27)= −ıω0 â â, ρ̂ +2
It is assumed that ρ̂ can be expanded as,Zρ̂ =d2 αP (α, t) |αi hα| ,(4.28)

64
where P (α, t) is a quazi-probability function called the P function, and |αi is a coherent state.
Eq. 4.28 is substituted into Eq. 4.27, yielding,Z∂d α |αi hα| P (α, t) =
∂t
2

Z

n


d2 αP (α, t) − ıω0 â† â |αi hα| − |αi hα| â† âγ
+2â |αi hα| â† − â† â |αi hα| − |αi hα| â† â2o+γ n̄ â |αi hα| â† + â† |αi hα| â − â† â |αi hα| − |αi hα| ââ† .
(4.29)

To simplify Eq. 4.29, the fact that,12†|αi = e− 2 |α| eαâ |0i ,(4.30)

is used to show that,∂|αi hα| =∂α∂  −|α|2 αâ†α∗ âee |0i h0| e∂α=â† − α∗ |αi hα| ,(4.31)and that also,∂|αi hα| =∂α∗∂  −|α|2 αâ†
α∗ â
e
e
|0i
h0|
e
∂α∗

= |αi hα| (â − α) .

(4.32)

Eqns. 4.31 and 4.32 can then be used to simplify the various terms in Eq. 4.29, which then becomes,
Z∂d α |αi hα| P (α, t) =
∂t
2

Z

n γ
 ∂
d2 αP (α, t) −+ ıω0 α
2
∂α
γ

2 o
∂
∂
−
− ıω0 α∗ ∗ + γ n̄
|αi hα| .
2
∂α
∂α∂α∗

(4.33)

Eq. 4.33 is now integrated by parts, assuming that P (α, t) vanishes sufficiently rapidly at infinity
so that the surface terms can be dropped. Then,Z∂d α |αi hα| P (α, t) =∂t2Zd2 α |αi hα|

n γ

 ∂

α
∂α
 ∂
γ
∂2 o
∗
+
− ıω0
α
+
γ
n̄
P (α, t).
2
∂α∗
∂α∂α∗
2

+ ıω0

(4.34)

65
Eq. 4.34 implies that,
n γ
 ∂
γ
 ∂
∂
∂2 o
∗
α
+
γ
n̄
P (α, t).
P (α, t) =
+ ıω0
α+
− ıω0
∂t
2
∂α
2
∂α∗
∂α∂α∗

(4.35)

This is the Fokker-Planck equation for the damped harmonic oscillator in the P representation.
Eq. 4.35 can be put into a more convenient form by defining α = x + ıy and substituting,nγ∂
P (α, t) =
∂t
2

4.4








∂
∂
∂
∂
n̄ ∂ 2∂2 ox+y + ω0y−x ++P (α, t).
∂x
∂y
∂x
∂y
2 ∂x2 ∂x2

(4.36)

Equivalence of Langevin and Fokker-Planck Equations
It is now demonstrated a Fokker-Planck equation for a probability distribution has an equiv-

alent Langevin equation for the random variable of that probability distribution. To do this, wewill need to take a brief aside into Stochastic calculus [34].
A Langevin equation is a stochastic differential equation of the form,dx= a(x, t) + b(x, t)ξ(t),dt(4.37)

where ξ(t) is a rapidly fluctuating random function of time. We assume that ξ(t) is idealized in thesense that we assume that ξ(t) and ξ(t0 ) are statistically independent. Mathematically this is thestatement that the time average of ξ(t) and ξ(t0 ) are delta correlated,
ξ(t)ξ(t0 ) = δ(t − t0 ).

(4.38)

hξ(t)i = 0,(4.39)We also assume thatsince any non-zero mean can be absorbed into the definition of a(x, t). This idealized random
process is called white noise.
It can be shown thatZ t0

dt0 ξ(t0 ) = W (t),(4.40)

66
where W (t) is the Wiener process. Since the Wiener process is not differentiable, Eq. 4.37 does not
actually exist. However, the integral equation
Z t
Z t
0
0
x(t) − x(0) =dt a(x, t ) +ξ(t0 )dt0 b(x, t0 )
0
0
Z t
Z t
=
dt0 a(x, t0 ) +dW (t0 )b(x, t0 ),0(4.41)0

where
dW (t) ≡ dW (t + dt) − dW (t) = ξ(t)dt(4.42)

does exist, but comes with its own difficulties, namely that the last term is not well defined. Thiscan be seen by looking at the integral before dt is allowed to be come infinitesimal. Consider
Sn =

n
X

G(τi ) {W (ti ) − W (ti−1 )} ,(4.43)

i=1

where τi is some value between ti and ti−1 . The crux of the whole situation is that the value of Sn
is dependent on the value of τi that we choose. This issue, which does not occur with an ordinary
Riemann sum defined integral, is caused by the discontinuous nature of W (t). The choice τi = ti−1 ,defines what is known as the Ito stochastic integral. Other choices, such as τi = 12 (ti + ti−1 ), which
defines the Stratonovich stochastic integral, will not be considered here. By using the definition of
the Ito integral, it can be shown that in the Ito interpretation of a stochastic integral,dW (t)2 = dt,
dW (t)2+N= 0,

where N is a positive integer. This causes the chain rule from ordinary calculus not to hold. To
find the Ito calculus chain rule, consider the differential of a function f (x(t)) of a stochastic variablex(t),df [x(t)] = f [x(t) + dx(t)] − f [x(t)]1= f 0 [x(t)]dx(t) + f 00 [x(t)]dx(t)2 + . . .
2= f 0 [x(t)] {a(x, t)dt + b(x, t)dW (t)}1
+ f 00 [x(t)] {a(x, t)dt + b(x, t)dW (t)}2 + . . .
2(4.44)

67
To first order in dt, this gives,
1
df [x(t)] = f 0 [x(t)] {a(x, t)dt + b(x, t)dW (t)} + f 00 [x(t)]b(x, t)dW (t)22
1
= f 0 [x(t)] {a(x, t)dt + b(x, t)dW (t)} + f 00 [x(t)]b(x, t)dt.
2(4.45)

This is known as Ito’s formula.
The equivalence between a stochastic differential for a stochastic variable x(t) and a FokkerPlanck equation for the stochastic variable’s probability distribution P (x(t)) can now be demonstrated. Consider the time development of the average value of an arbitrary function f (x(t)),
dhdf (x(t))ihf (x(t))i =
dt

 dt
1 00
0
=
f [x(t)]a(x, t) + f [x(t)]b(x, t) ,2

(4.46)

where the average value of the term multiplied by dW (t) is zero because dW (t) = ξ(t)dt andhξ(t)i = 0. The average value of a stochastic variable is found by integrating the stochastic variablemultiplied by its conditional probability density P (x0 , t0 |x0 , t0 ) from time x0 to x. Following thisprocedure, Eq. 4.46 becomes,dhf (x(t))i =
dt

Z




1 00
dx f (x)a(x, t) + f (x)b(x, t) P (x, t|x0 , t0 ).
20Integrating by parts and discarding the surface terms results in,


Z

d
1 d2 d2
hf (x(t))i = dxf (x) − [a(x, t)P (x, t|x0 , t0 )] +
b(x,
t)
P
(x,
t|x
,
t
)
.
0 0
dt
dx
2 dx2

(4.47)(4.48)The time development of the average value of an arbitrary function f (x(t)) can also be written asZddhf (x(t))i = dxf (x) P (x, t|x0 , t0 ).
(4.49)
dt
dt
Equating Eqns. 4.48 and 4.49 and remembering that f (x) is an arbitrary function of x results in,dd
1 d2 
P (x, t|x0 , t0 ) = − [a(x, t)P (x, t|x0 , t0 )] +b(x, t)2 P (x, t|x0 , t0 ) .
2
dt
dx
2 dx

(4.50)

Eq. 4.50 is the Fokker-Planck equation for the conditional probability density P (x, t|x0 , t0 ) for
the stochastic variable x(t). This shows the equivalence between a stochastic differential equation,
Eq. 4.37, for a stochastic variable x(t) and a Fokker-Planck equation for the probability distributionof that variable.

68

4.5

C-number Langevin Equations For Coherent State Amplitudes
Eq. 4.50 is now used to find the stochastic differential equations for α(t) that is equivalent

to Eq. 4.36. Since P (α, t) is complex, there are two stochastic equations; one for the real part ofP (α, t),dx = −

γ
2



r



r

x + ω0 y dt +n̄
dW1 (t),2(4.51)n̄dW2 (t).
2(4.52)and one for the imaginary part of P (α, t),dy = −

γ2y − ω0 x dt +These equations can be rewritten as one complex Langevin equation asdα = −

γ2

√+ ıω0 αdt + n̄dη(t),(4.53)

wheredη(t) =

(dW1 (t) + ıdW2 (t))√.
2

Eq. 4.53 can be used to generate trajectories for α. These trajectories can then be used to calculate
various expectation values of interest, even n time expectation values because,(â† (t1 ) · · · â† (tm )â(tm+1 ) · · · â(tn )Z=

d2 αP (α) {α∗ (t1 ) · · · α∗ (tm )α(tm+1 ) · · · α(tn )}
N

=

1 X ∗
αi (t1 ) · · · αi∗ (tm )αi (tm+1 ) · · · αi (tn ),N(4.54)

i=1

where αi (t) is a single trajectory, and N is the number of trajectories.
The results of a computer simulation that solves Eq. 4.53 and uses Eq. 4.102 to find â† (t)â(t)and g (1) (τ ) = Re â† (tss + τ )â(tss ) are shown in Fig. 4.1 and Fig. 4.2, respectively, where tss isthe time in which the system reaches steady state.
In both figures, the analytical solutions, Eqns. 3.45 and 3.55 are shown in black. In a single
trajectory, which are shown in blue in both figures, the system is undergoing a random walk, sothat the number of excitations randomly diffuses in time. When many of these random trajectories

69
are averaged together, the correct evolution for â† (t)â(t) and g (1) (τ ) are approached, which isshown in red, where 1000 random trajectories are averaged over. The figures show good agreement
between the analytical results and the numerical results from the above Langevin method.

4.6

A More Practical Derivation of Quantum Langevin Equations
In theory, the quantum Langevin equations, along with the corresponding diffusion coeffi-

cients, for an arbitrarily complicated system can be derived as they were in section 4.1. In practice,
however, there is a much simpler way to derive them, which proves very valuable with increasingsystem complexity. Langevin equations consist of two separate terms, called drift and diffusion
terms. The drift terms consist of the Hamiltonian dynamics, and the decay due to the coupling of
the system to a large reservoir. The diffusion terms contain the fluctuations caused by the couplingof the system to a large reservoir. The generalized Einstein relations provide a way to obtain the
diffusion terms if the drift terms and all of the second order moments are known. The drift terms,
and the second order moments can be calculated using the master equation that describes thesystem. Then, the diffusion terms can be found using the Einstein relations. We begin by deriving
the generalized Einstein relations.

4.6.1

The Generalized Einstein Relations

Consider the Langevin equations for a general system described by the operators Âµ ,dÂµ (t) = D̂µ (t) + F̂µ (t),dt

(4.55)

where D̂µ (t) are the drift terms and F̂µ (t) are the diffusion terms, or noise operators. We assume
that the noise operators obeyDEF̂µ (t) = 0,(4.56)

EF̂µ (t)F̂ν (t0 ) = 2Mµν (t)δ(t − t0 ),(4.57)

andD

where Mµν (t) are called the diffusion coefficients.

70"a †a#
4

3

2

1

2

4

6

8

10

Κt

Figure 4.1: The c-number Langevin equations are used to find â† (t)â(t) for a damped harmonic
oscillator with n̄ = 1, and α(0) = 2. The analytic expression for â† (t)â(t) is shown in black, asingle trajectory is shown in blue, and an average over 1000 trajectories is shown in red.

71

g!1"!Τ"
2.0
1.5
1.0
0.5

1

2

3

4

5

ΚΤ

#0.5



Figure 4.2: The c-number Langevin equations are used to find g (1) (τ ) = Re â† (tss + τ )â(tss )for a damped harmonic oscillator with n̄ = 1, and α(0) = 2, and tss is the time in which the systemreaches steady state. The analytic expression for g (1) (τ ) is shown in black, a single trajectory is
shown in blue, and an average over 1000 trajectories is shown in red.

72
We substitute Eq. 4.55 into the identity,
Z t

dt0

Âµ (t) = Âµ (t − ∆t) +t−∆tdÂµ (t0 ),dt(4.58)

multiply by F̂ν (t) on the left, and then arrive at,
D

E Z t
Âµ (t)F̂ν (t) = Âµ (t − ∆t)F̂ν (t) +
E

D

0

dt

D

0

0

E

D̂µ (t )F̂ν (t) + F̂µ (t )F̂ν (t) .

(4.59)

t−∆t

The first term on the right hand side of Eq. 4.59 is zero on physical grounds. This is becausea system operator at a given time can have no correlation with a reservoir operator at a laterEDtime, hence this term is factorizable, and F̂µ (t) = 0. A similar line of reasoning shows thatDE
D̂µ (t0 )F̂ν (t) = 0, unless t = t0 , so the lower limit of integration on this term can be changed tot, and hence this term is zero. Then,E Z t
DÂµ (t)F̂ν (t) =ED
dt0 F̂µ (t0 )F̂ν (t) .

(4.60)

t−∆t

Using Eq. 4.57, combined with the fact that a delta function lying on one of the limits of integration
causes the value of the integral to be half of what it would be otherwise, results in,DEÂµ (t)F̂ν (t) = Mµν (t).

(4.61)In an analogous manner to the method just described, it can be shown that,DE
F̂µ (t)Âν (t) = Mµν (t).

(4.62)

Now consider,
 


d 
d 
Âµ (t) Âν (t) + Âµ (t)Âν (t)
dt
dt
D
E D
E D
E D
E
=
D̂µ (t)Âν (t) + Âµ (t)D̂ν (t) + F̂µ (t)Âν (t) + Âµ (t)F̂ν (t) ,(4.63)

Ed D
Âµ (t)Âν (t)=dt

(4.64)
where in the last line Eq. 4.55 and linearity has been used. Substituting Eqns. 4.61 and 4.62 into
Eq. 4.64 and solving for Mµν (t), results in,
2Mµν (t) =

E D
E D
E
d D
Âµ (t)Âν (t) − Âµ (t)D̂ν (t) − D̂µ (t)Âν (t) .
dt(4.65)

73
Eq. 4.65 is called the generalized Einstein relation. It shows how to calculate the diffusion coefDE
ficients Mµν (t) if the drift terms and the equations of motion for Âµ (t)Âν (t) are known. The
generalized Einstein relation proves invaluable when the system under consideration is complex andthe exact form of the noise operators is unknown.

4.6.2

Obtaining the Drift Terms

We specialize to the system of N damped two level atoms with decay rate γ and frequency
ωa at zero temperature, driven by a strong classical field that couples to the atoms with Rabi
frequency Ω. By Eq. 3.119, this system is described by,
N 

i γX
d
ıh
σ̂j+ σ̂j− ρ̂(t) + ρ̂(t)σ̂j+ σ̂j− − 2σ̂j− ρ̂(t)σ̂j+ ,ρ̂(t) = − Ĥ, ρ̂(t) −dt~
2(4.66)

j=1

with
Ĥ =

N

N

j=1

j=1

~ωa X z ~Ω X − ıωc tσ̂j +(σ̂j e+ σ̂j+ e−ıωc t ),2
2

(4.67)

where,Jˆz =NXσ̂jz ,

j=1

Jˆ± =

N
X

σ̂j± ,

j=1

where σ̂j− , σ̂j+ , σ̂jz are the Pauli matricies for the jth atom.
We seek equations for Jˆ− , Jˆ+ , Jˆz that are of the form of Eq. 4.55,d ˆJµ (t) = D̂µ (t) + F̂µ (t),dt

(4.68)

where µ = {−, +, z}. From Eq 4.56, the drift terms are,
D



E
d
D̂µ (t)=Âµ (t)dtd= tr Âµ ρ̂(t) ,dt(4.69)

74
where the last line is valid in the Schrödinger picture. Using Eq. 4.66 we find,






d
σ̂j− ρ̂(t)

= tr

dt



=

d
σ̂j− ρ̂(t)dt−ıωa −γ  D − E ıΩ −ıωc t zσ̂j + eσ̂j ,22

(4.70)

so that,
D


N 
1 X
− d
σ̂j ρ̂(t)
N
dt
j=1


ıΩ −ıωc t ˆγ ˆ
J− + eJz ,=−ıωa −
2
2

E
D̂− (t)=(4.71)and therefore,

γ ˆ
ıΩ
D̂− (t) = −ıωa −J− + e−ıωc t Jˆz .
22(4.72)

Likewise, it can be shown that,
γ ˆ
ıΩ
ıωa −
J+ − eıωc t Jˆz ,
2
2




D̂z (t) = −γ 1̂ + Jˆz + ıΩ Jˆ− eıωc t − Jˆ+ e−ıωc t ,D̂+ (t) =

(4.73)(4.74)

where,
N

1̂ =

1 X
1̂j .
N

(4.75)

j=1

We move into an interaction picture by defining,˜Jˆ− = eıωc t Jˆ− ,˜Jˆ+ = e−ıωc t Jˆ+ ,(4.76)and combine these with Eqns 4.72 – 4.74 to obtain,ıΩγ  ˜ˆ
J− + Jˆz ,22γ  ˆ˜ıΩ ˆ=ıδ −J+ − Jz ,
2
2




˜
˜
= −γ 1̂ + Jˆz + ıΩ Jˆ− − Jˆ+ ,

˜
D̂− =
˜
D̂+
˜
D̂z
where δ = ωa − ωc .



−ıδ −

(4.77)

75
4.6.3

Obtaining the Diffusion Coefficients

Now that the drift terms are known, Eq. 4.65 can be used to find the diffusion coefficients.
We first find,d D ˆ ˜ˆ E D ˆ ˜ E D ˜ ˜ˆ E
Jz J+ − Jz D̂+ − D̂z J+
dt

 
ıΩ ˆ
1 d D ˆ˜ E
γ  ˜ˆ
ˆ=J+ − Jz
J− − Jz
ıδ −
N dt
2
2


D

 E
˜
˜
˜
− −γ 1̂ + Jˆz + ıΩ Jˆ− − Jˆ+ Jˆ+ ,2M̃z+ =

(4.78)

where Eqns. 4.77, and the fact that,˜Jˆz Jˆ+ =N

N

1 X X z ˜+σ̂i σ̂jN2
i=1 j=1N=

1 X ˜+
σ̂j
N2

=

1 ˜ˆ
J+ ,
N

j=1

(4.79)

have been used. Also, to obtain Eq. 4.79, Eqns. 3.121 have been employed. From Eq. 4.68, it is
D E D E˜˜deasy to see that dt
Jˆ+ = D̂+ , so that after substitution of Eq. 4.77, and further simplification,
Eq. 4.78 becomes,
2M̃z+




γ  ˜ˆıΩ ˆ
1γ  ˜ˆıΩıδ −J− − Jz −ıδ −J− − 1̂
=
2
2
N
2
2






1
ıΩ
˜
˜
−
1̂ − Jˆz−γ Jˆ+ + Jˆ+ +N
2E
2γ D ˆ
˜
=
J+ .
N
1
N



(4.80)

Likewise, it can be shown,
2M̃−z =
2M̃zz =
2M̃−+ =

2γ D ˜ˆ E
J− ,
N

D E
2γ
˜
1 + Jˆz ,Nγ.
N(4.81)

76
Putting it all together, the quantum Langevin equations for N two level atoms driven by a strongclassical field in an interaction picture that rotates at the field frequency are,d ˜γ  ˜ˆ
ıΩ
˜
Jˆ− =
−ıδ −
J− + Jˆz + F̂−
dt
2
2


d ˜
γ
ıΩ
˜
˜
Jˆ+ =
ıδ −
Jˆ+ − Jˆz + F̂+
dt
2
2




d ˆ
˜
˜
˜
Jz = −γ 1̂ + Jˆz + ıΩ Jˆ− − Jˆ+ + F̂z ,dt(4.82)

where the noises are correlated according to,DE˜˜F̂µ (t)F̂ν (t0 ) = 2M̃µν δ(t − t0 ),(4.83)

with,−2M̃µν =



−

+

0

γ
N



+
0

z 0

0
D E
˜
2γ
Jˆ+
N

z
D E
˜
2γ
Jˆ−



N




0


D E
˜
2γ
1 + Jˆz.

(4.84)

N

It will be convenient in the next section to work with real quantities, so we define,˜˜˜Jˆx ≡ Jˆ+ + Jˆ−


˜
˜
˜
Jˆy ≡ −i Jˆ+ − Jˆ−˜˜˜F̂x ≡ F̂+ + F̂−


˜
˜
˜
F̂y ≡ −i F̂+ − F̂− .

(4.85)

Using these definitions, and Eqns. 4.82 and Eqns. 4.84, it is then trivial to show that,
γ ˜
d ˜
˜
˜
Jˆx = − Jˆx − δ Jˆy + F̂x ,
dt
2
d ˜
γ ˜
˜
˜
Jˆy = δ Jˆx − Jˆy − ΩJˆz + F̂y ,
dt
2


d ˆ
˜
˜
Jz = −γ 1̂ + Jˆz + ΩJˆy + F̂z ,dt(4.86)

where the noises are correlated according to,DE˜˜F̂µ (t)F̂ν (t0 ) = 2M̃µν δ(t − t0 ),(4.87)

77
with,
x
x
2M̃µν =



y
− iγ
N

γ
N



iγ
y

 D EN D E˜˜z Nγ
Jˆx + i Jˆy

γ
N

− iγ
N

D

E
D E
˜
˜
Jˆx + i Jˆy

z
D E
˜ˆ
˜
Jx − i Jˆy

D E
D E
˜ˆ 
˜ˆ
iγ
Jx − i Jy N
D E 
˜ˆ2γ.
N 1 + Jz
γ
N

D

E

(4.88)

4.7

A More Practical Derivation of the C-number Langevin Equations

4.7.1

Solving for the C-Number Diffusion Coefficients

There is also a more practical way to derive c-number Langevin equations than the method
described in Sections 4.3 - 4.5. We assume that there exists a set of c-number equations that havea correspondence with the quantum Langevin equations,dγ
Jx = − Jx − δJy + Fx ,dt
2dγ
Jy = δJx − Jy − ΩJz + Fy ,dt
2dJz = −γ (1 + Jz ) + ΩJy + Fz ,dt

(4.89)

with,hFµ (t)i = 0,(4.90)

where µ = {−, +, z}. Also the noises are correlated according to,Fµ (t)Fν (t0 ) = 2Mµν (t)δ(t − t0 ),(4.91)

where Mµν is a 3 × 3 an unknown matrix of c-numbers, as a function of time. The entire problem
thus boils down to finding the c-number diffusion matrix elements Mµν that make the c-number
Langevin equation have a correspondence with the quantum Langevin equation. To make thiscorrespondence, we force the second order moments from the quantum equations to be equal to
the second order moments from the c-number equation. However, since c-numbers commute, and78
quantum operators do not, we must choose an ordering for the quantum equations before thecorrespondence can be made.
To illustrate this point, consider the following c-number second order moment, constructedfrom Eq. 4.89,
   dd
Jx Jy + Jx
Jy
dt
dt
D γ
 E D Eγ=
− Jx − δJy + Fx Jy + Jx δJx − Jy − ΩJz + Fy
2
2


d
hJx Jy i =dt= −γ hJx Jy i + δ (− hJy Jy i + hJx Jx i) − Ω hJx Jz i + 2Mxy ,where Mxy =D˜ ˜F̂x JˆyE=D(4.92)E˜ ˜
Jˆx F̂y . Since c-numbers commute, and quantum operators do not,Eq. 4.92 can correspond with equal validity to the quantum second order moment,   d ˜˜˜ d ˜ˆ
Jˆx Jˆy + Jˆx
Jy
dt
dt
D γ
 E D 
E
γ ˜
˜
˜
˜ ˜
˜
˜
˜
=
− Jˆx − δ Jˆy + F̂x Jˆy + Jˆx δ Jˆx − Jˆy − ΩJˆz + F̂y
2
2
D
E
 D
E D
E
D
E
˜
˜
˜
˜
˜
˜
˜
= −γ Jˆx Jˆy + δ − Jˆy Jˆy + Jˆx Jˆx − Ω Jˆx Jˆz + 2M̃xy ,

E
d D˜
˜
Jˆx Jˆy=dt

(4.93)or its anti-ordered counterpart,   d ˜˜˜ˆ d ˜ˆˆˆ
Jy Jx + Jy
Jxdtdt
 E D  γ
E
D
γ ˜
˜ ˜
˜
˜
˜
˜
˜
=
δ Jˆx − Jˆy − ΩJˆz + F̂y Jˆx + Jˆy − Jˆx − δ Jˆy + F̂x
2
2
D
E
 D
E D
E
D
E
˜ ˜
˜ ˜
˜ ˜
˜
= −γ Jˆy Jˆx + δ − Jˆy Jˆy + Jˆx Jˆx − Ω Jˆz Jˆx + 2M̃yx ,

E
d D˜
˜
Jˆy Jˆx=dt(4.94)or a symmetric combination of the two,
 D
E 1D
E
1 ˜
˜
˜
˜
Jˆx Jˆy +Jˆy Jˆx
2
2
D
E
 D
E
D
E 
D
E
˜ ˜
˜ ˜
˜ ˜
˜
= −γ Jˆx Jˆy + δ − Jˆy Jˆy + Jˆx Jˆx− Ω Jˆx Jˆz + M̃xy + M̃yx . (4.95)

Ed D ˜ˆ ˜
Jx Jˆy
=
dt
s

d
dt

s

s

s

s

Here, we choose to make the symmetric correspondence, so that,
2Mxy ≡ M̃xy + M̃yx
= −
= 0.

iγ iγ
+
N
N
(4.96)

79
We choose the symmetric correspondence for several reasons. First, in this case, and in many othercases, it allows for the simplest diffusion matrix 2Mµν . Second, in the symmetric correspondence,
2Mµν is, by construction, automatically a symmetric matrix. This fact allows an easy algorithm
to be performed to find the noises that reproduce 2Mµν for numerical simulations. Finally, again
by construction, there will not be any commutator pieces left to be absorbed into the diffusion coefficients when the symmetric correspondence is made. To demonstrate this, consider the following
c-number second order moment,   dd
Jy Jy + Jy
Jy
dt
dt
D

E D 
E
γ
γ
=
δJx − Jy − ΩJˆz + Fy Jy + Jy δJx − Jy − ΩJˆz + Fy22dhJy Jy i =dt= 2δ hJx Jy i − γ hJy Jy i − 2Ω hJz Jy i + 2Myy ,(4.97)and the following quantum second order moment,   d ˜˜˜ˆ d ˜ˆˆˆJy Jy + Jy
Jy
dt
dt

E D 
E
D
γ ˜
γ ˜
˜
˜
˜ ˜
˜
˜
=
δ Jˆx − Jˆy − ΩJˆz + F̂y Jˆy + Jˆy δ Jˆx − Jˆy − ΩJˆz + F̂y
2
2
D
E
D
E
D
E
D
E
D
E
˜ ˜
˜ ˜
˜ ˜
˜
˜
= δ Jˆx Jˆy + δ Jˆy Jˆx − γ Jˆy Jˆy − Ω Jˆz Jˆy − Ω Jˆy Jˆz + 2M̃yy . (4.98)

E
d D˜
˜
Jˆy Jˆy
=
dt

D
E
D
E
˜
˜
Since, for instance, hJz Jy i can not correspond to both Jˆy Jˆz and Jˆz Jˆy simultaneously, commutation relations must be used to put Eq. 4.98 in terms of one or the other before a correspondence
can be made. The resulting left over commutator piece will have to be absorbed into the definitionof 2Myy . However, in the symmetric correspondence, since,E
1 d D ˜ˆ ˜ˆ E 1 d D ˜ˆ ˜ˆ E
d D˜˜Jˆµ Jˆν =Jµ Jν +Jν Jµ ,
dt
2 dt
2 dt
s

(4.99)

every term will automatically have its symmetric counterpart, so that there will be no left overcommutator pieces to absorb into Mµν . Therefore, in the symmetric correspondence, we can
automatically set,
2Mµν = M̃µν + M̃νµ .

(4.100)80
In the symmetric correspondence, therefore, the c-number diffusion matrix is,xx
2Mµν =



γ
N



y
 0

z Nγ hJx i

y

z

0

γ
N hJx i

γ
N








2γ
N (1 + hJz i) .
γ
N hJy iγ
N hJy i(4.101)

Eqns. 4.89 and Eq. 4.101 are the c-number Langevin equations for N two level atoms driven by astrong field, corresponding to symmetric ordering. These c-number Langevin equations correspondto a Fokker-Planck equation for a Wigner quazi-probability distribution.

4.7.2

Numerical Simulation of C-number Langevin Equations

The final problem that must be overcome before a numerical simulation can take place is how
to choose the noises in Eqns. 4.89 so that Eq. 4.101 is reproduced. Since Eq. 4.101 is symmetric, it
is easy to diagonalize. In this diagonal basis, the noises are simply the square root of independent
real Wiener processes. Then the inverse transformation is applied, and in the original basis, thenoises will be linear combinations of the independent real Wiener processes, and they will reproduce
Eq. 4.101.
Now that this final problem has been solved, Eqns. 4.89 with noises that satisfying Eq. 4.101
can be used to generate trajectories. These trajectories can then be used to calculate various
symmetric expectation values of interest,NJˆµ (t1 ) · · · Jˆν (tm ) s =1 X
Jµ (t1 ) · · · Jν (tm ),N(4.102)i=1

where Jµ (t) is an operator for a single trajectory, and N is the number of trajectories.

D

The results of a computer simulation that solves Eqns. 4.89 subject to Eq. 4.101 to findEhDEiJˆz (t) and S(ω), the normalized Fourier transform of g (1) (τ ) = Re Jˆ+ (tss + τ )Jˆ− (tss )areshown in Fig. 4.3 and Fig. 4.4, respectively. Here, tss is the time in which the system reaches
steady state.

81

#Jz$

1.0

0.5

2

4

6

8

10

Γt

"0.5

"1.0

FigureD
E4.3: The c-number Langevin equations for symmetric ordering are used to find the inversionJˆz (t) for resonance fluorescence for N = 100 atoms, with Ω = 5γ and δ = 0. The analyticDEˆ
expression for Jz (t) for N = 1 is shown in black, a single trajectory is shown in blue, and an
average over 1000 trajectories is shown in red.

82

S"Ω#

1.0

0.8

0.6

0.4

0.2

!10

!5

5

10

Ω!Γ

Figure 4.4: The c-number Langevin equations for symmetric ordering are used to find S(ω), thenormalized Fourier transform of g (1) (τ ) for resonance fluorescence for N = 100 atoms, with Ω = 5γ
and δ = 0. The analytic expression for S(ω) for N = 1 is shown in black, a single trajectory is
shown in blue, and an average over 1000 trajectories is shown in red.

83
In both figures, the analytical solutions for single atom resonance fluorescence [85] are shownin black. In a single trajectory, which are shown in blue in both figures, the system is undergoinga random walk, so that the number of excitations randomly diffuses in time. When many ofDEthese random trajectories are averaged together, the correct evolution for Jˆz (t) and S(ω) areapproached, which is shown in red, where 1000 random trajectories are averaged over.

4.8

Conclusion
We have seen how Langevin equations for quantum operators are derived. The approxi-

mations that go into this derivation are equivalent to the approximations used in the derivationof the master equation. In the superradiance limit, the cavity field described by the quantum
Langevin equations can be adiabatically eliminated, and used to derive the superradiance master equation. Since Langevin equations of quantum operators are difficult to solve numerically, a
Fokker-Planck equation for a quazi-probability distribution was derived. Here, we used the normal
ordered quazi-probability distribution, called the P-function, which can be used to calculate normal ordered moments. This distribution was shown to be equivalent to a set of c-number Langevinequations by using methods of Ito stochastic calculus. These c-number Langevin equations weresolved for many trajectories, which can be used to find general n-time expectation values. Since
the derivation of the Fokker-Planck equation becomes quite complicated for systems of increasing
complexity, a method to more practically derive c-number Langevin equations was described. This
involved deriving the generalized Einstein relations, which allow the diffusion, or noise terms ofa set of quantum Langevin equations to be known, if the diffusion terms are known. The drift
terms and second ordered moments were found using a master equation, and the Einstein relations
were then used to find the diffusion terms. Then, the method to find a correspondence to a setof c-number equations for a given ordering was described. We specialize here to use a symmetriccorrespondence, because in this ordering, the form of the equations are best suited for the numericalalgorithm that we describe. Although this method may seem like unnecessary machinery to solve
such a simple problem, it really proves worth its weight in gold when more complex systems are84
analyzed involving large numbers of atoms.

Chapter 5

The Crossover from Lasing to Steady State Superradiance

The bulk of this chapter is being prepared for submission to Physical Review A. [81]Since its first demonstration in 1960 [60], the laser has had a profound impact on fundamental
science research and numerous areas of society. Lasers are integral to many fields of research andtechnological applications, ranging from atomic and molecular physics research, atomic clocks, theglobal positioning system, nuclear fusion research, biology research, medicine, and consumer electronics. Although many different types of lasers exist, with their characteristic parameters (suchas intensity, power, linewidth, physical size) spanning many orders of magnitude, all lasers share
a common conceptual foundation. A laser is a cavity quantum electrodynamics (QED) system
consisting of a gain medium inside an optical cavity. We will oftentimes refer to the gain medium
as “atoms” for brevity. Lasers typically operate in the good cavity regime of cavity QED where
the linewidth of the cavity is much narrower than the bandwidth of the gain medium. The atoms
generate a coherent electromagnetic field in the cavity by means of stimulated emission [73]. Stimulated emission is a quantum mechanical interference effect in which the presence of a large numberof photons in a particular mode of a light field increases the probability that an atom will emit
into that mode. The energy emitted into the cavity field has to be replenished by some repumping
mechanism to achieve steady state operation. In a laser, the macroscopic phase information thatis associated with the coherence of the generated radiation is encoded in the light field.
Around the same time as the laser was first demonstrated, the effect of superradiance was86
predicted [19], and soon thereafter, experimentally demonstrated. Superradiance is a quantum
mechanical interference effect in which correlations between atoms lead to collective emission. Superradiance has most commonly been considered as a pulsed phenomenon. Atoms in an ensemble
are prepared in the excited state. Spontaneous emission into one or a few spatial modes is thenenhanced via growth of atom-atom correlations. However, it has been known for some time thatsuperradiance can also occur in steady state [66, 65, 64, 22] by placing the atomic ensemble inside acavity. In contrast to lasers, superradiance in steady-state occurs in a cavity with a much broader
linewidth than the atomic linewidth. This regime is referred to as the bad-cavity limit of cavity
QED [84, 57, 14, 39, 40]. The radiation produced in steady state superradiance is also coherent.
However, in contrast to a laser, the coherence is encoded in the atomic medium. Progress has
recently been made towards the experimental verification [10, 9] of this proposal.
An important application of lasers is as a stable local oscillator for optical atomic clocks
and precision spectroscopy. These lasers rely on stabilization against reference cavities. The mostadvanced such lasers reach linewidths of < 0.1Hz corresponding to quality factors of Q > 1015 [16].
The principal limiting factor in the way of further improvement of these local oscillators is thermal
vibrations of the dielectric coatings in the cavity mirrors [51]. To overcome this technical challenge,
researchers have proposed an alternate approach using an active system based on steady statesuperradiance on a clock transition to create an even more stable light source [66, 14]. However,this proposal has challenges of its own. First of all, in spite of the enhancement that occurs due tosuperradiance, the produced intensity is orders of magnitude lower than for a conventional laser.
Second of all, perturbations of atomic transition frequencies can potentially lead to phase andfrequency perturbations in the generated field.
In this chapter we develop a unified theory for lasers and steady state superradiance. In thisunified theory, lasers and steady state superradiance are the extreme limits of a continuous crossoverand the theory continuously interpolates between the two. This allows us to directly compare andcontrast lasers, steady state superradiant systems, and systems in the crossover region using thesame language. Our analysis thus further clarifies the qualitative and quantitative differences

87
between lasing and steady state superradiance. From the perspective of applications the unified
theory enables us to determine the optimal system for ultra stable local oscillators and precision
measurement applications.
We analyze the model using different levels of approximation: An exact method using SU(4)operators, a semi-classical method based on c-number Langevin equations, a semi-classical method
based on the cumulant expansion method, a quantum phase diffusion model for the field amplitude, and a mean field model. The different approaches provide insight into different aspects ofthe problem. Highly simplified models like the mean field equations and phase diffusion yield aqualitative understanding of the general characteristics of systems throughout the crossover. By
comparison between the approximations we can differentiate between truly critical physical effects
and less important details. We find that fluctuations and correlations are essential for the noiseproperties of the system (e.g. the linewidth of the generated light) but that the fluctuations andcorrelations can be modeled semi-classically. Comparison with the exact SU(4) method for small
numbers of atoms shows that c-number Langevin equations provide an accurate description of thesystem. Due to their much smaller computational complexity, we are then able to use the c-number
Langevin equations to quantitatively study much larger systems relevant for experiment.
The rest of this chapter is organized as follows. In section 5.1 we summarize the physicalmodel upon which our analysis is based. In 5.2 we discuss several approximation methods. We
compare the approximations with one another to determine their accuracy and to evaluate theirability to capture the various physical signatures. In section 5.3 we define a crossover parameter
which characterizes the relative importance of stimulated emission to collective atomic effects in acavity QED system. In section 5.4 we discuss our results on the crossover.

88

5.1

Model
We model our system as a collection of N two-level atoms inside a single mode optical cavity

using the quantum Born-Markov master equation to describe the open quantum system,
i
d
1 h
ρ̂ =
Ĥ, ρ̂ + L̂ [ρ̂] ,
dt
i~

(5.1)

where,
Ĥ =NN
~ωa X z~Ω X  † −σ̂j + ~ωc â† â +â σ̂j + σ̂j+ â ,
2
2
j=1

(5.2)

j=1

and L̂ [ρ̂] denotes the Liouvillian superoperator.
The Hamiltonian Ĥ describes the coherent evolution of the coupled atom cavity system,
where ωa is the frequency between the two levels of each of the atoms, and ωc is the frequency ofthe cavity mode. The Pauli spin matrices for the atoms are σ̂j+ ,σ̂j− and σ̂jz , and â is the annihilation
operator of the cavity mode. The atom-cavity coupling rate is Ω. That we take this to be the samefor all items implies assumption of the equivalence of the atom-cavity coupling at the positions of
all of the atoms. This could be achieved by strongly confining the ensemble in a optical latticepotential at the antinodes of the cavity modefunction.
The incoherent evolution is described by the Liouvillian L̂ [ρ̂],κ †â âρ̂ + ρ̂â† â − 2âρ̂â†2N
γ X + −−σ̂j σ̂j ρ̂ + ρ̂σ̂j+ σ̂j− − 2σ̂j− ρ̂σ̂j+
2L̂ [ρ̂] = −

j=1
N

−


w X − +σ̂j σ̂j ρ̂ + ρ̂σ̂j− σ̂j+ − 2σ̂j+ ρ̂σ̂j− ,
2
j=1

+

N

1 X z zσ̂j ρ̂σ̂j − ρ̂ ,
2T2

(5.3)

j=1

where ρ̂ is the system’s density matrix, κ is the decay rate of the cavity, γ is the incoherent decayrate of the atoms, w is the incoherent repumping rate of the atoms, and T12 is the inhomogeneousdephasing rate.

89

5.2

Solution Methods
Eq. 5.1 becomes increasingly complex with increasing atom number N , since the associatedHilbert space scales as 2N . Exact simulation methods, such as the quantum jump and quantum
state diffusion methods, are limited to N ∼ 10 with the cavity field basis truncated to just a few
photons.
The SU(4) solution method [43] [88] exploits an underlying permutation symmetry in Eq. 5.1
to drastically reduce the number of basis states needed to describe the problem, allowing systems
with N ∼ 100 to be treated. However to treat this many particles, the cavity field basis mustalso be truncated to just a few photons. In order to treat systems which have a larger number ofphotons, we combine the SU(4) method with the quantum jump method [17] [27] [69]. Even withthe ability to describe problems with larger photon number, N ∼ 100 is still too few particles to
describe many realistic experimental situations.
We therefore also describe the c-number Langevin equation method [75, 55], which approximate the quantum Langevin equations that are equivalent to Eq. 5.1. These equations do not scalewith N in the same way, and therefore can often be used to treat systems with large N where an
exact treatment is not possible. Another method that does not scale with N in the same way is themean field solutions of the quantum Langevin equations, along with an analytic solution for the
linewidth of the power spectrum, which assumes Gaussian phase fluctuations and neglects amplitude fluctuations [39, 40]. Finally, the cumulant expansion method [56] does not scale exponentiallywith N , and can therefore also be used to treat large systems. These methods are now describedin detail.

90
5.2.1
5.2.1.1

Langevin Theory
Quantum Langevin Equations

The Quantum Langevin equations that are equivalent to Eq. 5.1 are are derived using themethod demonstrated in Chapter 4.6, and are given by,d1
iN Ω −
Ŝ + F̂ a ,â = − (κ + 2iωc )â −
dt
2
2
d −
1
iΩ z
Ŝ = − (Γ + 2iωa ) Ŝ − +âŜ + F̂ − ,
dt
2
2




d z
z
† −
+
Ŝ = −(w + γ) Ŝ − d0 + iΩ â Ŝ − âŜ
+ F̂ z ,
dt

(5.4)
(5.5)
(5.6)

where δ = ωa − ωc . We have defined the collective operators,
N

Ŝ

−

=

Ŝ z =

1 X −
σ̂k ,
N
1
N

k=1
N
X

σ̂kz ,

k=1+−
where Γ ≡ w + γ + T22 , d0 = w−γ
w+γ and where Ŝ and Ŝ are Hermitian conjugates of one another,Ŝ + = (Ŝ − )† . The noise operators F̂ µ have zero mean and their second order correlations are givenbyDEF̂ µ (t)F̂ ν (t0 ) = 2Dµν δ(t − t0 ) .

(5.7)

The diffusion matrix elements Dµν can be calculated using the Einstein relations [67],†2Daa = κ

D E
1
1 
w+
1 + Ŝ z
2D
=
N
T2

D E
1
1 
2D−+ =
γ+
1 − Ŝ z
N
T2
D
E
2w
2γ D + EŜ +Ŝ
2D+z = −2Dz+ =
N
N
D
E
D
E
2γ
2w
2D−z =
Ŝ −2Dz− = −Ŝ −NND E 2w 
D E
2γ
2Dzz =
1 + Ŝ z+
1 − Ŝ z .
N
N
+−

(5.8)

91
5.2.1.2

c-number Langevin equations for numerical simulations

The quantum Langevin equation are operator valued stochastic differential equations. As
such they are not suited for practical computations. To obtain practical equations we construct asemi-classical theory by replacing the operators in the quantum Langevin equations by c-numbers,d1iN Ω −
a = − (κ + 2iωc )a −S + F a,
dt
2
2
d −
1
iΩ z
S = − (Γ + 2iωa ) S − +aS + F − ,
dt
2
2


d z
S = −(w + γ) (S z − d0 ) + iΩ a† S − − aS + + F z ,
dt

(5.9)
(5.10)(5.11)

where the omission of the hats over the variables signifies that they are c-numbers. The noise
terms F a , F − , and F z should be interpreted according to the rules of Ito calculus. The noises arecorrelated according toF µ (t)F ν (t0 ) = 2D µν δ(t − t0 ) .

(5.12)

It is easier to construct the semi-classical equations by introducing real variables according
toq̂ =Ŝ x =1 †â + â ,2
1 +Ŝ + Ŝ − ,
2


1  †
â − â ,
2i

1  +
y
Ŝ =
Ŝ − Ŝ − .
2ip̂ =(5.13)(5.14)The equations of motion in terms of these variables are
d
q
dt
d
p
dt
d x
S
dt
d y
S
dt
d z
S
dt

= −κq − 2ωc p − N ΩS y + F q ,(5.15)= −κp + 2ωc q + N ΩS x + F p ,(5.16)= −ΓS x − 2ωa S y + ΩpS z + F x ,(5.17)= −ΓS y + 2ωa S x − ΩqS z + F y ,(5.18)= −(w + γ) (S z − d0 ) + 2Ω (qS y − pS x ) + F z .
(5.19)

92
The correspondence between the semi-classical and quantum mechanical Langevin equations
is established by requiring that they produce identical equations for first and second moments ofthe system operators. Comparison of the first moments (i.e. the expectation values of the equationsof motion) leads to
hF q i = hF p i = hF x i = hF y i = hF z i = 0 ,(5.20)

Comparison of the second moments allows us to find the classical diffusion matrix elements
D µν from the quantum mechanical ones. To make this procedure well defined we have to choosea specific ordering of the quantum mechanical operators. We choose to make the correspondence
using symmetric ordering defined by the symmetric expectation valueDÂµ ÂνEs=

1 D µ ν E D ν µ E,Â Â + Â Â
2(5.21)

where Âµ and Âν are generic system operators. We point out that in this formulation, the classicalLangevin equations are equivalent to a Fokker-Planck equation for the Wigner quasi-probability
distribution. A tedious but straightforward calculation yields
κ4Γ
2D xx =
4Nκ
4
Γ
2D yy =
4N
−w + γ D x E
2D xz = 2D zx =ŜN−w + γ D y E
2D yz = 2D zy =
Ŝ
N
D E
2 
2D zz =
(w + γ) + (−w + γ) Ŝ z .
N

2D qq =

2D pp =

(5.22)

We solve the stochastic differential equations (5.19) by means of the explicit second orderweak scheme [52]. An ensemble of trajectories can be evolved simultaneously and the expectation
values appearing in Eq. 5.22 can be calculated from the ensemble. This allows the additive formof the explicit second order weak scheme to be used, which is simpler than the general form.
The noises F µ are found by means of
Fµ =

Xp
T
λν Mµ,νfν ,ν

(5.23)

93
where Mµ,ν is the orthogonal matrix that diagonalizes the diffusion matrix, λν are its eigenvalues,
and f ν are independent normalized Wiener processes. It was found empirically that the symmetrically ordered diffusion matrix is positive definite when the system is above the first threshold(defined in Sec. 5.2.1.3). Below this threshold, the symmetrically ordered diffusion matrix is not
positive definite, and divergent trajectories can occur. Typically, an ensemble of 1000 trajectories
can be used to achieve convergence to within a few percent.

5.2.1.3

Mean Field Equations

By taking expectation values of the semi-classical Eqns. 5.15–5.19 we obtain the so called
mean field equationsd
1
iN Ω −
a0 = − (κ + 2i(ωc − ω))a0 −
S ,
dt
2
2 0
d −
1
iΩ
S0 = − (Γ + 2i(ωa − ω)) S0− +a0 S0z ,
dt
2
2

d z
S0 = −(w + γ) (S0z − d0 ) + iΩ a∗0 S0− − a0 S0+ .
dt(5.24)(5.25)(5.26)The noise terms drop out because they have zero mean. The mean field equations capture many of√the most important features of the physical system because the noise terms scale as N while theexpectation values scale as N . In the limit of large numbers of atoms the noise terms are therefore
less important. Note that Eqns. 5.24–5.26 are written in the reference frame rotating at frequency
ω.
Eqns. 5.24–5.26 can be written in steady-state by setting the left hand sides to zero and thensolved in closed form. We findS0z =(κ + 2i(ωc − ω))(Γ + 2i(ωa − ω))N Ω2(5.27)

for the steady state inversion. From this, the oscillation frequency of the atom-cavity coupledsystem w can be determined using the condition that S0z must be real, givingω=κωa + Γωc.
κ+Γ(5.28)

94
When δ = ωa − ωc  Γ, κ the inversion can be simplified to
S0z ≈

1
,
C

(5.29)

2Ω
where C ≡ NκΓ
is the generalized many-atom cooperativity parameter (generalized since here theeffective linewidth Γ includes the incoherent repumping w as well as the single atom linewidth γ
and dephasing 1/T2 ).
In the same limit of small detuning we find for the steady state photon number|a0 |2 =N (w + γ)1(d0 − ) .
2κC

(5.30)

1
N Ω2
−γ− .
2κ
T2

(5.31)The photon number has a maximum atw = wopt =In the limit where the collective decay rate CΓ is much larger than the single atom rates γ and T12
we find a simple expression for the maximum photon number,(|a0 |2 )opt =N 2 Ω2.
8κ2(5.32)The zeros of the intra cavity photon number Eq. 5.30 determine the thresholds of the system.
At the first threshold
w1 = γ,(5.33)

energy is supplied to the system at a high enough rate to sustain a macroscopic field amplitudein the cavity. The emergence of a coherent macroscopic field amplitude is accompanied by theformation of a collective atomic dipole. A second threshold occurs at
w2 =

N Ω2
.
κ

(5.34)At this point, S0z is close to unity, and the noise due to the strong repumping prevents the formationof either a macroscopic coherent state in the cavity or a macroscopic dipole in the atomic ensemble.

95
5.2.1.4

Phase Diffusion Linewidth

It is also possible to derive an analytic expression for the linewidth from the quantum
Langevin equations, Eqns. 5.4 – 5.6 by assuming that only phase fluctuations of â are important, and neglecting amplitude fluctuations. To this end, we follow closely to the derivation done in[39, 40]. Eq. 5.4 is differentiated, and substituted into it Eqns. 5.4 – 5.5 and the integral of Eq. 5.6,
to arrive at an equation for â alone,2¨ = − 1 (κ + Γ)â˙ − κΓ â + N Ω âŜ z + F̂ ,
â
2
4
4

(5.35)

where
Ŝ

z

Z t
=

0 −(w+γ)(t−t0 )dt e(w + γ) + F̂ z
0


2
d †
†
† a
a†
−
(â â) + κâ â − â F̂ − F̂ â,
N dt(5.36)and,F̂ =Γ a iN Ω −˙F̂ −F̂ + F̂ a .
2
2(5.37)

The annihilation operator â is decomposed according to,â = (a0 + ρ̂)eiφ̂ .

(5.38)

Above threshold, amplitude fluctuations are expected to be small, so thatDE
DEâ† (t)â(0) = a20 ei(φ̂(t)−φ̂(0)) .

(5.39)

After substituting Eq. 5.38 into Eq. 5.35, we take the imaginary part to first order in productsof operators, and find,11
˙¨φ̂ = − (κ + Γ)φ̂ + Im[F̂ ],2
a0(5.40)

where a factor of e−iφ has been absorbed into F̂ . Eq. 5.40 is then integrateed, assuming that (κ+Γ)
is large, to arrive at,2
φ̂(t) − φ̂(0) =a0 (κ + Γ)Z t


Γ a iN Ω −
dt Im
F̂ −
F̂ .
2
2
0
0



(5.41)

96
Since F̂ a and F̂ − are Gaussian, we can useDi(φ̂(t)−φ̂(0))

eE= e− 2 h(φ̂(t)−φ̂(0)) i .
12(5.42)

Therefore, we use Eq. 5.41, along with Eqns. 5.8 to find,D2(φ̂(t) − φ̂(0))

E


=

Γ
Ω2 κ
(C + 1)2(Cd0 − 1) (w + γ) (κ + Γ)2t,(5.43)

so that the linewidth ∆ν given by,
∆ν =

5.2.2

(C + 1)ΓΩ2 κ.
2(Cd0 − 1) (w + γ) (κ + Γ)2(5.44)

2nd Order Cumulant Theory

The cumulant theory is derived in the same way as in Chapter 3.7. Eq. 5.1 is used to findequations for expectation values of system observables, which are truncated to second order in acumulant expansion, resulting in the closed set of equations,
o
E
nD
d z
hσ̂1 i = −(w + γ) (hσ̂1z i − d0 ) + iΩ â† σ̂1− − âσ̂1+dt
DE iΩ  D E
d D † −E1â σ̂1= − (Γ + κ − 2iδ) â† σ̂1− +â† â hσ̂1z i
dt
2
2

1
+ −
z
+ (hσ̂1 i + 1) + (N − 1) σ̂1 σ̂2
2
o
nD
E
Ω
d
σ̂1+ σ̂2− = −Γ σ̂1+ σ̂2− + hσ̂1z i â† σ̂1− − âσ̂1+
dt
2i
D
E
D
E
o
d
N Ω nD † − Eâ† â = −κ â† â +â σ̂1 − âσ̂1+ .
dt
2i(5.45)(5.46)(5.47)(5.48)

This second order cumulant expansion, along with the quantum regression theorem, can be used in
the same way as in chapter 3.7.2, to find a closed set of equations for two-time expectation values,


 

κ
iN Ω
†
†
d  â (τ )â(0)   − 2 + iδ2
  â (τ )â(0) (5.49)
=.
iΩhσ̂1z idτ++Γσ̂1 (τ )â(0)−2σ̂1 (τ )â(0)− 2These equations are written in a frame rotating at the atomic frequency. They are important, since
the Fourier transform of â† (τ )â(0) is the spectrum of the light inside our cavity QED system.

97
By using Laplace transform methods, it is found that the solution of the â† (τ )â(0) equation,of Eq. 5.49 is dominated by a single exponential with decay constant,11∆ν = (κ + Γ − 2iδ) −22q(−κ + Γ + 2iδ)2 + 4N Ω2 hσ̂1z i,(5.50)and hence the spectrum is a Lorentzian with linewidth given by Eq. 5.50.

5.3

Characterization of the Crossover
The crossover from superradiance to lasing is characterized by a transition from coherence

encoded in an atomic ensemble to coherence encoded in the light field. The key parameter in
identifying the regime is the ratio of the number of photons to the number of atoms when thesystem is operated at its most efficient point in terms of output intensity and spectral frequency
width. With this motivation, we introduce a crossover parameter asξ≡(|a0 |2 )opt,
N

(5.51)

that is, the dimensionless ratio of the intracavity photon number to the number of atoms. The
parameter ξ quantifies the relative importance of stimulated emission to collective atomic spontaneous emission. If ξ  1, the system is in the bad cavity or superradiant regime. If ξ  1 thesystem is in the good cavity or laser regime. In the crossover or intermediate region, ξ ∼ 1, andthe system possesses features of both.
The mean field equations allow us to rewrite this expression in an alternate way that illuminates the role of the system parameters in determining the operating regime. From Eq. 5.32,
which is applicable in the typical case where wopt  γ, 1/T2 , Eq. 5.51 can be rewritten as an
approximately equivalent definition
ξ≈

N Ω2
.
8κ2

(5.52)The interpretation of this is that the crossover is also characterized by the ratio of the collective
coupling between the many atom ensemble and the cavity mode, i.e. N Ω2 , to the decay rate κ ofthe cavity mode.

98

5.4

Results
In this section we present numerical results for field intensity, inversion, line widths, andcorrelations throughout the crossover from lasing to superradiance. We begin by considering thedifferent solution methods for small numbers of atoms were we can compare with the exact SU(4)Monte-Carlo results. This comparison shows that the semi-classical theory gives an accurate description of first and second moments of the system operators. This is an interesting result in itself
besides verifying the semi-classical methods because it demonstrates that the system is not strongly
correlated and we don’t need to account for genuinely quantum fluctuations.
With the validity of the semi-classical method established we then consider systems withlarge numbers of atoms. These systems are experimentally relevant and allow us to more clearly
understand the behavior of the system in the large N limit.
For simplicity, in all of these simulation results, we consider the parameter regime in which
the single atom rates γ and 1/T2 are negligible compared to the cavity loss rate κ and use Eq. 5.52
to define ξ. The results of the c-number Langevin simulations are an average over 1000 trajectories.
In each trajectory, the trajectory values were used in the diffusion matrix to compute the noise.
This method is justified by its agreement to the method that uses expectation values in the diffusion
matrix, and also by its agreement with the exact SU(4) results.

5.4.1

Simulations with small atom numbers

In order to determine the validity of the approximate solution methods, we begin by comparing them to the exact solutions. We use N = 40 atoms for the comparison. This is a smallenough to still be tractable by exact SU(4) Monte-Carlo simulations and at the same time it islarge enough to expect the approximate solution methods to be reasonably accurate.
Specifically, Fig. 5.1 shows several observables obtained using the mean field Langevin method,
the phase diffusion method, the c-number Langevin method, the cumulant expansion method, andexact SU(4) Monte-Carlo simulations for three different values of the crossover parameter: ξ = 0.2,

Superradiance99+<! z >
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
0.0 0.5 1.0 1.5 2.0
w/w opt+-<! 1 ! 2 >
0.12
0.10
0.08
0.06
0.04
0.02
0.00
0.0 0.5 1.0 1.5 2.0
w/w opt

<a a>
8

(2)

6

!" /#
1

4

0.1

2

10

0
0.0

0.5

1.0 1.5
w/w opt-2

0.0 0.5 1.0 1.5 2.0
w/w opt

2.0

Lasing

Crossover+<! z >
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
0.0 0.5 1.0 1.5 2.0
w/w opt<! 1 ! 2 >
0.12
0.10
0.08
0.06
0.04
0.02
0.00
0.0 0.5 1.0 1.5 2.0
w/w opt<! z >
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
0.0 0.5 1.0 1.5 2.0
w/w opt<! 1 ! 2 >
0.12
0.10
0.08
0.06
0.04
0.02
0.00
0.0 0.5 1.0 1.5 2.0
w/w opt

(a)

+

-

<a a>
40

G (0)
2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.0 0.5 1.0 1.5 2.0
w/w opt
(2)

!" /#
1

30
20

0.1

10

10

-2

0
0.0 0.5 1.0 1.5 2.0
w/w opt

0.0 0.5 1.0 1.5 2.0
w/w opt

G (0)
2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.0 0.5 1.0 1.5 2.0
w/w opt
(2)

+

-

(b)

0
0.0 0.5 1.0 1.5 2.0
w/w opt

0.0 0.5 1.0 1.5 2.0
w/w opt

G (0)
2.0
1.8
1.6
1.4
1.2
1.0
0.8
0.0 0.5 1.0 1.5 2.0
w/w opt

(c)

(d)

(e)

+

!" /#
1

<a a>
150
100

0.1

50

10

-2

Figure 5.1: Comparison of the different solution methods in the superradiance (ξ = 0.2), crossover(ξ = 1), and lasing (ξ = 5) regions for N = 40 and C = 1. The analytic Langevin solution is shownin red (light gray), the 2nd order cumulant solution is shown in blue (dark gray), The exact SU(4)solution is shown by grey triangles, and the c-number Langevin simulation results are shown by
black circles. The observables considered are (a) the inversion hσ̂ z i, (b) the correlation between
atoms σ̂1+ σ̂2− , (c) the intracavity photon number â† â , (d) the linewidth ∆ν, and (e) the intensity
correlation function G(2) (0).

100
ξ = 1, and ξ = 5. These values of ξ put the system in the superradiance, crossover, and lasing
parameter regions, respectively. We note that the SU(4) method becomes more computationallyintensive as ξ is increased, since there are more photons as ξ is increased, and hence more basis
states need to be tracked. Therefore, for ξ = 5, we used a method that combines the SU(4) approach with the quantum jump method [81]. The results using this method that are shown are an
average over 500 trajectories.
Fig. 5.1 shows excellent agreement between c-number Langevin and the exact SU(4) theory
in all parameter regions for all of the considered observables. Therefore, the c-number Langevintheory can be relied upon for larger atom numbers inaccessible to direct numerical simulation.
Figures 5.1 (a) and (c) show that the mean field equations are accurate near the collective
emission peak, w/wopt = 1, but they are less accurate outside that region. Figure 5.1 (d), showsthat the phase diffusion model for the linewidth also agrees in the region around w/wopt = 1, but
disagrees outside that region, where the phase diffusion approximation breaks down. Although
they do not quantitatively agree with the exact SU(4) method, the analytic solutions obtained by
the mean field and phase diffusion models capture the correct qualitative behavior of the system.
The 2nd order cumulant solution works well qualitatively in all parameter regions, but asseen in Fig. 5.1 (b) and (c) in the superradiance row, and in column (d), there is a quantitative
disagreement between this theory and the SU(4) and c-number Langevin theories.
This quantitative disagreement is due to the fact that there is a difference between a Gaussian
theory, one in which all moments factorize into products of no higher than second order moments,
and a theory that uses Gaussian noise. This difference can be seen by considering the third ordermoment â† âσ̂iz . According to the cumulant theory, â† âσ̂iz= â† â hσ̂iz i, a fact which wasessential to use in deriving the closed the set of equations used to find the spectrum, Eqns. 5.49.
In the Langevin theory, however, â† âσ̂iz 6= â† â hσ̂iz i, a fact we have checked with our Langevinsimulations. Even though the difference between â† âσ̂iz

and â† â hσ̂iz i decreases as N1 , it isprecisely the N1 terms that gives rise to the linewidth according to the cumulant theory.
The merit of the 2nd order cumulant theory lies in its simplicity. It is simple to derive the

101

< z>
1.0
0.8
0.6
0.4
0.2
0.0
-0.2
-2
10

+

0.1

+

<a a>/(N )
1.0
0.8
0.6
0.4
0.2
0.0
10 10 -2
0.1-

< 1 2 >
0.12
0.10
0.08
0.06
0.04
0.02
0.00
0.1
1
10 10 -2
w/w opt

1

w/w opt

(2)

G (0)
2.0
1.5
1.0
1

0.5
10 10 -2

0.1

1

10

w/w opt

w/w opt(a) Universal
! /"
1
0.1
-2
10
-3
10
-4
10
-5
10
-2
10

0.1

1

w/w opt(b) Superradiance

! /"
1
0.1
-2
10
-3
10
-4
10
-5
10
10 10 -2

0.1

1

! /"
1
0.1
-2
10
-3
10
-4
10
-5
10
10 10 -2

0.1

1w/w opt

w/w opt(c) Crossover(d) Lasing

10

Figure 5.2: Solutions using the various methods in the superradiance (ξ = 0.1), crossover (ξ = 1),
and lasing (ξ = 10) regions for N=10000 and C = 0.1. For 1/T2 = 0, the analytic Langevinsolution is shown in solid red (solid light gray), the 2nd order cumulant solution is shown in solid
blue (solid dark gray), and the c-number Langevin simulation results are shown by black circles. For
1/T2 = 15 wopt , the analytic Langevin solution is shown in dashed red (dashed light gray), the 2ndorder cumulant solution is shown in dashed blue (dashed dark gray), and the c-number Langevinsimulation results are shown by black diamonds. (a) All observables considered except linewidth∆ν show universal behavior in the superradiance, crossover, and lasing regions, after appropriatescaling. The red (light gray) curves have been artificially widened in order to be seen underneaththe blue curves. (b) ∆ν/κ in the superradiance region (c) ∆ν/κ in the crossover region, (d) ∆ν/κin the lasing region.

102
cumulant equations and it is not computationally intensive to solve them numerically, while thetheory still captures the correct qualitative behavior of the system.

5.4.2

Simulations with large atom numbers

Now that an agreement between the c-number Langevin and exact SU(4) theories has been
established, we study more experimentally realistic systems with N = 10000 using the semi-classicaltheory based on the c-number Langevin equations. We also include the mean field Langevin theory,the phase diffusion model for the linewidth. The results of these simulations are shown in Fig. 5.2.
We consider both the case of vanishing inhomogeneous broadening, 1/T2 = 0, as well as 1/T2 =wopt /5.
As seen in Fig. 5.2 (a), when 1/T2 = 0, the inversion hσ̂ z i, the correlation between atoms
σ̂1+ σ̂2− , the intracavity photon number â† â , and the intensity correlation function G(2) (0) allshow universal behavior in the superradiance, crossover, and lasing regimes after appropriate scaling. For 1/T2 = wopt /5 these observables still show universal behavior after appropriate scaling,
and the values do not change significantly from the 1/T2 = 0 case.
It is worth noting that, even though σ̂1+ σ̂2− has universal behavior throughout the crossover,typical lasers operate just above threshold in terms of the scaled quantities, i.e. w  wopt . Just
above threshold the atom-atom correlations are typically very small,σ̂1+ σ̂2− opt  1 .

(5.53)

Thus atom-atom correlations are usually not important in lasers.
The linewidth ∆ν, however, does not show universal behavior in the superradiance, crossover,and lasing regimes. As seen in Fig. 5.2 (b), in the superradiance region, when 1/T2 = 0, ∆ν/κ isflat in the region of w/wopt < 1. In contrast, the linewidth in the lasing regime, shown in Fig.5.2(d), linearly decreases as w/wopt increases towards unity. This is the typical Schawlow-Townes
behavior, which can be seen by considering κ/ â† â using Eq. 5.30. In the crossover region, shownin Fig. 5.2 (c), we see that for w/wopt  1 , ∆ν/κ is flat, and as w/wopt approaches unity, ∆ν/κ

103

! /"
1
0.1
-2
10 -3
10 -4
10 -5
10 -6
10
-2
10 0.1+<a a>/N
100
10
1
0.1
1

10 10010-2

10-2 0.1

1

10 100w/w opt , c

w/w opt , c

(a)

(b)

Figure 5.3: Comparison of (a) linewidth and (b) Intracavity intensity for a system operated in the
crossover (ξ = 1) shown in blue solid, lasing (ξ = 10) shown in red dashed, and far lasing (ξ = 100)
shown in magenta dotted, regions. wopt,c is the optimum w value in the crossover region. For all
systems, N = 10000.

104
starts to linearly decrease as in the lasing regime. Therefore, a system in the crossover region
displays characteristics of both superradiance and lasing.
When 1/T2 is increased to 1/T2 = wopt /5, Fig. 5.2 (b) shows that ∆ν/κ increases for
w/wopt  1, but it is not significantly affected by 1/T2 as w/wopt approaches unity. Note that thelatter is the pumping region where the emitted light intensity is largest. Near w ∼ wopt we obtain
both the largest intensity and the best frequency stability.
In the crossover region, seen in Fig. 5.2 (c), when 1/T2 = wopt /5, ∆ν/κ increase for w/wopt 
1, where the system is displaying superradiant behavior. On the other hand, when w ∼ wopt where
the system starts to display lasing behavior, ∆ν/κ is insensitive to the inhomogeneous broadeningof the atomic line.
Figure 5.2 (d) shows that ∆ν/κ does not increase in the lasing region when 1/T2 = wopt /5.
Rather the linewidth decrease in the region slightly below w/wopt = 1, when compared to the
1/T2 = 0 case. This reduction has also been observed for smaller atom numbers using the exact
SU(4) code. It can be shown that the value of 1/T2 that maximally reduces the linewidth isw√ .
1/T2 = 1+opt
2

As mentioned previous, a laser normally operates just above threshold, where w/wopt  1.
As seen in Fig. 5.3 (a), for the same pump strength w, and a fixed N , a system in the crossover
region can operate with w/wopt = 1, allowing the crossover system to have a linewidth orders ofmagnitude smaller than the linewidth of the lasing system. Fig. 5.3 (b) shows that this improvement
in linewidth can be achieved without paying the penalty of a reduced output intensity. Anotherway of saying this is that, for the same output intensity, one can achieve greater brightness of thelight source by using atoms with a smaller linewidth.

5.4.3

Effect of Cavity and Atomic Level Instabilities on the Spectrum

The effect of instabilities in the cavity frequency and atomic frequency on the linewidth is
now investigated. Cavity frequency instability can be caused by fluctuations in cavity length dueto the thermal fluctuations of the cavity mirrors, which are unavoidable in an experiment. The

105

dΩ
dΩc
1.0
0.8
0.6
0.4
0.2
0.01

dΩ
dΩa
1.0
0.8
0.6
0.4
0.2
0.1

1
Ξ
(a)

10

100

0.01

0.1

1
Ξ

10

100

(b)

Figure 5.4: Instability in the atom-cavity system frequency ω with respect to (a) the cavity frequency ωc and (b) atomic frequency ωa as a function of crossover parameter ξ106
energy levels in an atom can also be shifted by stray electromagnetic fields. The frequency of the
combined atom-cavity system ω lies somewhere between the bare atom and cavity frequencies, andis given by Eq. 5.28.
If ω is varied, the resultant time averaged linewidth will be an average over these values,
and will hence be much broader than the quantum limited linewidth. Therefore, the derivativeof ω with respect to ωc , and with respect to ωa , gives us an idea as to how much of an effect an
instability in these frequencies will have on broadening the linewidth. These derivatives are plottedin Fig. 5.4.
Fig. 5.4 (a) tells us that on the lasing side of the crossover, ω is shifted proportionally to theshift in ωc , while on the superradiant side, ω is insensitive to shifts in ωc . Conversely, in Fig. 5.4(b) it can be seen ω is shifted proportionally to the shift in ωa on the superradiance side of thecrossover, while on the lasing side, ω is insensitive to shifts in ωa . In the intermediate regime(ξ = 1), ω is not as sensitive to a shift in the ωc as it is in the laser regime, and it is not as sensitive
to a shift in the ωa , as it is on the superradiance side.

5.5

Conclusion
A model that is capable of describing a system operating in any parameter region of thecrossover between steady state superradiance and lasing has been introduced, and three solution
methods to that model were demonstrated. The first method solves a set of c-number Langevinequations that approximate the quantum Langevin equations that are equivalent to our model.
This method was seen to have outstanding agreement with the exact SU(4) method. Also, analytic
expressions for the expectation values of the quantum Langevin equations were derived. Alongwith these solutions, an analytic solution for the linewidth of the power spectrum, which assumes
Gaussian phase fluctuations, and neglects amplitude fluctuations was derived. It is demonstrated
that these analytical results agree with SU(4) and c-number Langevin methods in the region inwhich the value of re-pumping is such that the system has the maximum number of intra-cavity
photons, but disagree far from this region. Finally, a closed set of equations to second order in a

107
cumulant expansion were derived, and shown to agree qualitatively with the SU(4) and c-number
Langevin methods.
Although a system in the lasing parameter region is capable of possessing the smallest
linewidth, typically lasers operate with a re-pumping rate that puts them just above threshold,which does not allow the smallest linewidth for the system to be realized. For the same re-pumping
rate, a system operating in the crossover region will be much farther above threshold, potentiallyallowing for a linewidth orders of magnitude smaller than in the laser region. Also, a crossover
system will have a larger intracavity intensity than a superradiant system, making the former moreexperimentally accessible.
It was also demonstrated that a system in all regions of the crossover is insensitive to atomic
dephasing as long as the re-pumping rate is smaller than the dephasing rate. Also, as the crossoveris traversed into the lasing region, atomic dephasing can actually cause a reduction in the linewidth.
Finally, it was shown that the linewidth on the superradiance side is insensitive to fluctuations
in cavity length, and that these fluctuations become more important as the crossover is traversed
toward the laser side. Conversely, it was shown that the linewidth on the lasing side is insensitive tofluctuations in the atomic energy levels, and that these fluctuations become more important on the
superradiance side. A system in the crossover region will be less sensitive to cavity frequency fluctuations than a lasing system, and it will be less sensitive to atomic frequency instabilities than asuperradiant system. Both of these fluctuations cause a broadening of the linewidth. These conclusions are significant since fluctuations in the cavity frequency are a principal limitation preventing
further reduction of the linewidth in today’s most ultrastable lasers.

Chapter 6

Laser Stabilization using Saturated Absorption in a Cavity QED System

The bulk of this chapter was published in Physical Review A as Phys. Rev. A 92,
013817 (2015) [80].

Today’s ultra-precise and accurate atomic clocks continue to make important contributions
to fundamental physics as well as applied technology. Atomic clocks have imposed significant
constraints on the drift of fundamental constants [58, 45, 35], may have the potential to enhance
the sensitivity of gravitational wave detectors, and have provided ultimate tests of the general
theory of relativity [74]. With the current stability of optical clocks at the 1×10−18 level, there areprospects for applying atomic clocks for the detailed mapping of the Earth’s gravity field [38, 11].
A highly stabilized laser is an integral component of high precision measurements, such as
optical atomic clocks and precision spectroscopy. Current technology for achieving highly phasestable laser sources relies on locking a laser to a high-Q reference ultra low expansion (ULE) glass
cavity [23, 91, 48]. The phase stability of this method is currently limited by the thermal noise
induced in the mirrors, spacers, and coatings of the reference cavity [68], but has been significantly
reduced over the past few years with new engineered materials [49, 62].
As an alternative approach to overcoming the thermal noise problem, it was recently proposed[63] to lock the laser to the saturated resonance feature exhibited by a collection of atoms with anultra-narrow electronic transition trapped in an optical cavity. Here, the atoms were assumed to be
trapped in an optical lattice inside the cavity. Due to the narrow atomic line, such a system would

109
typically operate in the parameter region corresponding to the bad cavity limit of cavity QED.
There, the atomic linewidth is significantly narrower than the cavity linewidth. In contrast to thereference cavity stabilization method described above, the cavity QED method offers a distinct
advantage since no drift compensation is needed.
The cavity QED system exhibits optical bistability in the intracavity intensity [12, 24] where
several solutions exist for the steady-state intracavity field. Working at an input intensity in theregion where bistability is present would in principle allow the greatest degree of stabilization [63].
However, it is not practical to work in the bistable region since quantum and classical fluctuations
between the semi-classical eigenmodes cause the system to be dynamically unstable. Therefore, one
is restricted to working at input intensities above the bistability, where the achievable stabilizationis orders of magnitude worse. Still, it was shown [63] that phase stability corresponding to thesub-mHz level should be achievable.
Recently [86], an experimental effort to demonstrate the cavity QED system was made byprobing the

1S0− 3 P1 intercombination line of 88 Sr atoms (i.e. γa /2π = 7.6 kHz) inside an

optical cavity. There, however, the atoms were not trapped in an optical lattice, but loaded into
the center of the cavity using a MOT, which was then turned off during probing. Typical MOT
temperatures correspond to a few millikelvin which is equivalent to a Doppler width of several MHz.
Considering the narrow 7.6 kHz line of the optical transition, this implies that motional effects will
be important.
In this chapter, we extend the many-atom cavity QED theory of [63] to include atomic
motion, and study its effect on the stabilization precision. In spite of the large Doppler effect, thestanding-wave nature of the cavity field induces sharp saturated absorption and dispersion features
to appear in the considered observables. These features are nestled in the center of the overall
Doppler broadened features [77, 32, 76, 59]. The stabilization that is achievable by utilizing thesesharp features is impeded by multi-photon scattering processes that occur when an atom’s velocity
matches one of its Doppleron resonances [1, 79]. The dependence of the stabilization on the numberof atoms and the temperature due to the Dopplerons is discussed. We demonstrate that the motion

110

FSR

ωc
Probe

Atomic sample
Cavity resonancelocked to probe laser

iDetuning

Laser stabilization on atomic resonance, i = 0

Figure 6.1: Schematic of the cavity QED experiment with a thermal sample of atoms with Doppler
width Γd , each of which have a narrow optical transition of width γa << Γd . The atoms are probedwith a carrier with frequency ωc , and two sidebands located at ωc ±FSR, where FSR is the free
spectral range of the cavity. The carrier frequency ωc , which is close to the atomic frequency ωa , islocked on the cavity mode frequency ω, while the sidebands are assumed far off resonance, typically∼ 104 atomic linewidths. By demodulating the light transmitted through the cavity at the FSR,
detection of the non-linear phase response of the transmitted light is achieved. This phase response
results in a photocurrent that serves as a frequency discriminating error signal that stabilizes thelaser frequency through the requirement ∆i = 0.

111
of the atoms causes the bistability region to disappear, so that no restrictions on input power arenecessary to avoid the dynamic instability that would otherwise result.

6.1

Model
We model our system as a collection of N two-level atoms inside a single mode optical cavity

using the quantum Born-Markov master equation to describe the open quantum system,id1 hĤ, ρ̂ + L̂ [ρ̂] ,
ρ̂ =
dt
i~

(6.1)

where,
Ĥ =

N
N




X
~∆ X zσ̂j + ~η â† + â + ~gj (t) â† σ̂j− + σ̂j+ â ,
2
j=1

(6.2)

j=1

and L̂ [ρ̂] denotes the Liouvillian.
The Hamiltonian H describes the coherent evolution of the coupled atom cavity system in aninteraction picture which rotates at the frequency of the cavity, and ∆ is the atom-cavity detuning.
The Pauli spin matrices for the atoms are σ̂j+ ,σ̂j− and σ̂jz , and â is the annihilation operator ofpthe cavity mode. Furthermore, η = (κPin )/(~ω) is the classical drive amplitude, where κ is the
decay rate of the cavity, Pin is the input power, and ω is the frequency of the cavity mode. The
atom-cavity coupling rate is gj (t) = g0 cos(δj t), where g0 is the maximum coupling amplitude, andδj = kvj is the Doppler shift in terms of the velocity vj of the jth atom, and wave number k of thelight.
The incoherent evolution describes the various forms of dissipation in this system and is
described by the Liouvillian L̂ [ρ̂],oκn †â âρ̂ + ρ̂â† â − 2âρ̂â†2Noγ Xn + −−σ̂j σ̂j ρ̂ + ρ̂σ̂j+ σ̂j− − 2σ̂j− ρ̂σ̂j+
2L̂ [ρ̂] = −

j=1

N

+

1 X z zσ̂j ρ̂σ̂j − ρ̂ ,
2T2
j=1

(6.3)

112
where ρ̂ is the system’s density matrix, γ is the spontaneous emission rate for the atoms, and T12 is
the inhomogeneous dephasing rate.
We derive Langevin equations corresponding to Eq. 6.1. Assuming that the classical drive η
is sufficiently strong, a mean-field description provides an accurate representation [37]. We define
D E
D E
D E
the mean values for the field α = −i hâi, and for the atoms, σj− = σ̂j− , σj+ = σ̂j+ , σjz = σ̂j+ ,which evolve according to the semiclassical evolution,α̇ = −κα + η +NX

gj (t)σj− ,

(6.4)

j=1
− (γp + i∆) σj− + gj (t)ασjz ,

(6.5)




σ̇jz = −γ σjz + 1 − 2gj (t) ασj+ + α∗ σj− ,(6.6)

σ̇j− =

where γp ≡ T12 + γ2 is the total decay rate of the atomic dipole. In the moving frame of reference ofthe jth atom, the cavity field appears as a traveling wave, containing two frequencies shifted aboveand below the cavity frequency by the Doppler shift δj (refer to Fig. 6.1).
It is convenient to approximate Eqns. 6.4 – 6.6 as a function of the continuous variable δ = kv,
Zα̇ = −κα + η + g0dδP (δ) cos(δt)σ − ,(6.7)

σ̇ − = − (γp + i∆) σ − + g0 cos(δt)ασ z ,

(6.8)


σ̇ z = −γ (σ z + 1) − 2g0 cos(δt) ασ + + α∗ σ − ,(6.9)

where P (δ) is the Maxwell velocity distribution of width δ0 , which is related to the temperature by
the equipartition theorem.
To solve this problem that intrinsically contains a bi-chromatic drive, we proceed in two
ways. Our first approach is to numerically integrate Eqns. 6.7 – 6.9, partitioning the integral into
finite-size velocity bins. The velocity partition must be chosen with care, since the system exhibits
Doppleron resonances, which have a strong dependence on the atomic velocity. Specifically, atlower velocity, more resolution in the partition is required.
Our second approach is semi-analytic, and involves a Floquet analysis [31, 2], in which we

113
expand σ − , σ + , and σ z in terms of their Fourier components,
σ− =

X

(l)

eilδt x1 ,

l

σ

+

=

X

=

X

(l)

eilδt x2 ,

l

σ

z

(l)

eilδt x3 ,(6.10)

l
(l)

(l)

(l)

where x1 ,x2 , and x3 are the amplitudes of the lth Fourier component. Upon substitution of
Eqns. 6.10 into Eqns. 6.7 – 6.9, equations for the amplitudes are found:

(l)

ẋ1

(l)

= − (i(∆ + lδ) + γp ) x1 +


g0 α  (l+1)
(l−1)
x3
+ x3
,
2
(6.11)

(l)

ẋ2

(l)

= (i(∆ + lδ) − γp ) x2 +
g0 α∗  (l+1)
(l−1)
x3
+ x3
,
2
(6.12)

(l)

ẋ3

(l)

= −γδl,0 − (ilδ + γ) x3


(l+1)
(l+1)
(l−1)
(l−1)
−g0 αx2+ α ∗ x1
+ αx2
+ α∗ x1.
(6.13)In order to find a steady state solution, we set the time derivatives in Eqns. 6.11 – 6.13 tozero, and substitute Eqns. 6.11 – 6.12 into Eq. 6.13, yielding,

(l)

(l+2)

0 = γδl,0 + al x3 + dl x3(l−2)+ bl x3,(6.14)

where δl,0 is a Kronecker delta, andg 2 |α|2al ≡ ilδ + γ + 0
2



1
1
1
1
+
+
+
Ql+1 Pl+1 Ql−1 Pl−1g02 |α|211
bl = dl−2 ≡+,2
Ql−1 Pl−1

,(6.15)(6.16)

where,Pl = i(lδ + ∆) + γp ,(6.17)

114
Ql = i(lδ − ∆) + γp .

(6.18)For a given α, Eq. 6.14 defines a tridiagonal linear system that can be solved by truncating
l at some finite value, and applying the Thomas algorithm for matrix inversion [70].
Since the atoms have motion, the condition for resonance between an atom and photon isachieved when the atomic frequency and photon frequency are offset by δ. However, higher ordermulti-photon processes, known as Doppleron resonances, involving 2n + 1 photons where n is an
integer, are also possible. The nth order Doppleron resonance corresponds to the terms of orderl/2 in the Floquet theory [1, 79]. This correspondence is why only even values of l can couple into
Eq. 6.14.
In steady state, Eq. 6.7 simplifies to,α=η g02 N
+
κ
2κ

Z



(−1)
(1)
dδP (δ) x1 + x1 .

(6.19)

1
Since x−1
1 and x1 depend on α, the self-consistent field amplitude α that solves Eq. 6.19 is foundnumerically by applying Newton’s method for root finding [70].
We have seen excellent agreement between the two previously described solution methods,
and for the remainder of the chapter, focus our attention on the Floquet solution, which mosttransparently illuminates the underlying physics.

6.2

Discussion of Steady State Solutions
We first consider the lowest order solution to Eq. 6.19, by truncating at l = 0, which means wehave not included higher order Doppleron processes. We have verified that this solution displays thecorrect physics qualitatively by comparing to higher order solutions that are truncated at increasing√√
values of l. We define scaled intracavity and input field amplitudes x ≡ α/ n0 , y ≡ η/(κ n0 ),where n0 = (γγp )/(4g02 ) is the saturation photon number. Eq. 6.19 then becomes,(Z
1 − i(∆ + δ)/γp
N C0
y = x 1+dδP (δ)2
2
4
1 + (∆+δ)+ |x|4 (1 + ξ + )γp2)!
1 − i(∆ − δ)/γp
+
,
2
|x|2
−)
1 + (∆−δ)
(1
+
ξ
+
4
γ2
p

(6.20)

115
g2

γ 2 +(∆±δ)2

where C0 ≡ κγ0p , and ξ ± = γp2 +(∆∓δ)2 .
p

It is interesting to consider the relation corresponding to Eq. 6.20 for a ring cavity system thathas a field traveling in only one direction, assuming equal intracavity power. The corresponding
expression is given by,

N C0
y =x 1+
2

ZdδP (δ)1 − i(∆ + δ)/γp2+ |x|2
1 + (∆+δ)
γ2

2

!
.

(6.21)

p

In the experiment described in [86], the measured observables are the cavity transmitted
power T ≡ |x/y|2 and transmitted phase shift φ ≡ arg(x/y) of the intracavity light relative tothe input light. Fig. 6.2 shows that the presence of ξ ± in Eq. 6.20 results in extra absorption anddispersive features (blue solid) around resonance in the transmission and phase shift, as compared
to a ring cavity field (red dashed) where ξ ± = 0.
These extra absorption and dispersive features are caused by the following: The distributionof atomic velocities results in a different Doppler frequency shift of the light for each atomic velocity
class, so that each velocity class will be resonant at a different detuning ∆. When |∆|  γp |x|,
the resonant velocity class of the atoms is interrogated by both components of the standing wave
field, whereas when |∆|  γp |x|, the resonant velocity class of atoms is interrogated by only one
component of the standing wave field. Thus, there is an increased saturation in atomic absorption,with a corresponding saturated dispersive feature, for |∆|  γp |x|. These sharp features are absent
from the traveling wave cavity situation where there is only one propagating field.
As seen in Fig. 6.2 (a), for low input intensities, there is no atomic saturation at any detuning.
In Fig. 6.2 (b), the atoms in the velocity class around resonance are saturated by both componentsof the field, and the velocity classes away from resonance are saturated by only a single componentof the standing wave field. Therefore, the features caused by the two component saturation andthe features caused by the single component saturation are clearly able to be distinguished. It canbe seen in Fig. 6.2 (c) that as the amount of saturation becomes large, the central feature becomespower broadened, but is still identifiable. Fig. 6.2 (d) shows the intracavity intensity for a given

116(a) 0.6
0.4
0.2
f 0
-0.2
-0.4
-0.6

0.8

T
0.4
0.2
-600 -400 -200 0 200 400 600
Dg p-600 -400 -200 0 200 400 600
Dg p(b)0.6
0.4
0.2
f 0
-0.2
-0.4
-0.6

d0g p

0.6

0.8
0.7

T 0.6
0.5
0.4

-600 -400 -200 0 200 400 600
Dg p-600 -400 -200 0 200 400 600
Dg p(c)0.6
0.4
0.8
0.2
0
T 0.7
-0.2
-0.4
-40 -20 0 20 40

0.4
0.2
f 0
-0.2
-0.4
-0.6

0.85
0.8
0.75
-40 -20 0 20 40

0.6-600 -400 -200 0 200 400 600
Dg p-600 -400 -200 0 200 400 600
Dg p

(d)

103

102
2

ÈxÈ 101
100
10-1
100

101

ÈyÈ2

102

103

Figure 6.2: Development of the extra absorption and dispersion features in the transmission (T)and phase shift (φ), as the input intensity is increased. Blue (dark gray) solid curves are fora standing wave cavity, red (light gray) dashed are for a traveling wave cavity. For all plots,
N C0 = 600 and δ0 /γp = 260, which in the case of 88 Sr, corresponds to a temperature of ∼ 15 mK.
(a) |y|2 = 5 × 10−1 , (b) |y|2 = 6 × 101 , (c) |y|2 = 9 × 102 ; The inset is zoomed in to emphasize thecentral features. (d) Input vs intracavity intensity at resonance. The black dots label the input
and intracavity intensities of (a) (b) and (c), and the dashed line is |y|2 = |x|2 for reference. Notethat there is no bistability.

117
input intensity with the values used in Figs. 6.2 (a-c) labeled by black dots.
In Fig. 6.3 and Fig. 6.4, a comparison between the solutions to Eq. 6.19, with l truncated at
l = 6, and the experimental results of [86] can be seen. Fig. 6.3 shows the transmitted phase shift
φ ≡ arg(x/y) of the intracavity light relative to the input light around resonance for several systemsof varying input probe laser power. Here is demonstrated a linear increase in phase slope as the
input laser power is reduced. Fig. 6.4 shows the transmitted phase shift for several systems withvarying numbers of atoms. Here is demonstrated a linear increase in phase slope as the number ofatoms is increased. In both figures, an excellent agreement between the theory and experimental
results can be seen.

6.3

Shot Noise Limited Laser Stabilization
For laser stabilization, the central part of the phase response close to the atomic resonanceserves as an error signal, see Fig. 6.2 (c), and allows for the generation of a feedback signal to
the laser frequency. Through the photodetector, the error signal is converted to a measurable
photocurrent. Any photocurrent measured by the photo detector should, in principle be zeroed by
an ideal feedback loop of infinite bandwidth to the laser frequency.
To determine the potential phase stability that could be achievable using our system, weconsider the shot noise limited stabilization linewidth assuming a strong local oscillator, as derived
in [63],

∆ν =

~ω


8πεPsig

∂φ
∂∆

2 =

C0
 2 ,∂φ4πεγ|x|2 ∂∆(6.22)∂φ
where Psig is the signal power, ε is the photo-detector efficiency and ( ∂∆) is the dimensionless phaseslope at resonance ∆ = 0. In the case that FM spectroscopy, such as NICEOHMS [86], is used forthe detection of the absorption or cavity transmitted phase, Eq. 6.22 must be modified as follows.
In general, in a configuration where sidebands are applied at the free spectral range of the cavity,
Eq. 6.22 must be multiplied with (1 + Psig /2Psideband ), where Psideband is the sideband power.

118

PHYSICAL REVIEW LETTE

PRL 114, 093002 (2015)
(a)

0.2

0.0

0.0

-0.2

-0.2

-0.4

Pin = 1950 nW-400 -200

0

200

400

Phase (radians)

0.4
0.2

Pin = 900 nW-400 -200

0

200

0.2

(a)

0.0

-0.2

400

-400 -200

0.4

(c)

0.2

0.0

(d)

0.0

-0.2
-0.4

-0.4

(b)

Phase (radians)

0.2

0.4

-0.2

Pin = 700 nW-400 -2000200∆/2π (kHz)FIG. 3 (color online).

400

-0.4

Pin = 650 nW-400 -200

0

200

400

∆/2π (kHz)Phase (radians)

Phase (radians)

0.4

0.2

(c)

0.0

-0.2
-400 -200∆/2π

Measured phase shift of the cavity-

FIG. 4 (color on
Figure 6.3: Measured phase shift of the cavity transmitted field when scanned across the atomic
transmitted
field
when
scannedacrosstheatomicresonance.
resonance. The input probe laser power, Pin , is progressively decreased from 1950 nW (a), 900whennW scanned acr
7
The
input
probe
laser
power
P
is
progressively
decreased
from(b), 700 nW (c) to 650 nW (d). The number ofinatoms is about 2.5 × 10 . Each point is an averagein the cavity is prof three data
points.
are theof 650
the theoreticalmodel.
1950
nW The(a), solid
900 linesnW (b),
700predictionsnW (c) tonW (d). The
number7
7
of atoms is about N cavity ¼ 2.5 × 107 . Each point is an average ofthree data points. The solid lines are theoretical predictions basedon our theoretical model.

laser stabilization. In Fig. 3 we show the phase signal for a
fixed number of atoms as a function of laser detuning fordifferent input powers in the range 650–1950 nW. For highinput powers we strongly saturate the dipole and power
broaden the central saturated absorption peak. As we
gradually lower the input power, the power broadening
is reduced, leaving the central phase feature with a larger

10 (b), 1.7 × 10
plots was 650 nW
The solid lines armodel. The centr

In conclusion
highly saturated
strontium atoms
transmission thr
but, apart from
response of theeffects while dis

119week ending
6 MARCH 2015

Pin = 900 nW
200

Pin = 650 nW
0

200

0.2

400

2π (kHz)

(a)

0.2

0.0

Ncavity = 2.5 x 107-400 -200

0.2

(b)

0.0

-0.2

400

Phase (radians)

0Phase (radians)

YSICAL REVIEW LETTERS

0

200

400

(c)

-400
0.2

0.0

Ncavity = 2.0 x 107

-0.2
-200

0

200

400

(d)

0.0

Ncavity = 1.7 x 107-0.2-400 -2000200∆/2π ( kHz )400

Ncavity = 1.2 x 107

-0.2
-400

-200

0

200

400

∆/2π ( kHz )t of the cavityFIG. 4 (color online). Measured phase shift of the probe lightomic resonance.
Figure 6.4:when
Measuredphaseacrossshift ofthe
theatomic
cavity transmitted
scannedacross the atomic
scannedresonance.field
Thewhennumberof atoms
7
resonance.inThein the cavity,
Ncavity , isfrom
progressively
7 (a), 2.0from
y decreased from
thenumbercavity ofis7atoms
progressively
decreased
2.5 × 10decreased
× 2.5 × 10
77(a), 2.0 × 10 7 (b), 1.7 × 10 (c)to 1.2 × 10 (d).7 The input power used for all plots was 650 nW.
(d). The number
7 (c)
10
(b),
1.7
×
10
to 1.2
× 10
The
input
power
used forall theoretical
Each point is an average of three data
points.
The(d).
solid
lines
are thepredictionsof thet is an averagemodel.
ofplots was 650 nW. Each point is an average of three data points.
predictions based
The solid lines are theoretical predictions based on our theoreticalmodel. The central slope scales linearly with atom number.

ase signal for aer detuning for
0 nW. For highole and power
peak. As wewer broadening
e with a larger

In conclusion, we have constructed a system dominated byhighly saturated multiphoton absorption with laser-cooledstrontium atoms coupled to a low finesse optical cavity. The
transmission through the cavity is altered by thermal effects
but, apart from a small decrease in slope, the central phase
response of the atoms remains relatively immune to theseeffects while displaying a high SNR owing to the cavity anddetection technique. The atomic phase signature was

120
Eq. 6.22 is the smallest, and hence the phase is most stable, when the product of the slope
around resonance and the intracavity intensity is as large as possible. The optimal input intensity,
which allows this product to be as large as possible, is the value used in Fig. 6.2 (c), and is labeledby the black dot at |y|2 ≈ 103 in Fig. 6.2 (d).
To achieve a quantitative agreement between theory and experiment [86], higher orders in l
in Eq. 6.19 must be included. These higher order terms correspond to Doppleron resonances, i.e.
multi-photon scattering processes between the atoms and cavity mode.
To study the importance of these higher order Dopplerons, we calculate the linewidth from
Eq. 6.22 at the optimum input power while varying the order of l at which the truncation occurs.
Fig. 6.5 shows the dependence of the linewidth at the optimum input intensity on the orderof l at which the truncation occurs for 3 different sets of parameters. The first set, shown in blue
circles, converges by l = 12. The linewidth calculated with up to l = 12 included before truncationis around 5 times larger than the linewidth with only l = 0 included. This shows that Doppleroneffects are crucial to include for a correct quantitative analysis of this system.
Shown in red squares, we increase the number of atoms by a factor of 10, and again calculate
the linewidth as a function of the order of l at which the truncation occurs. As a result of N being
increased, the optimum value of |y|2 is also increased. Convergence occurs around l = 16. Now, thedifference between the converged linewidth and the linewidth with only l = 0 included has increased
by a factor of around 2. This demonstrates that as N is increased, higher order Dopplerons play
an increasing role.
We also decrease the temperature by a factor of 10, and again calculate the linewidth at the
new optimum value of |y|2 , as shown by the black diamonds. Even though the optimum |y|2 occurs
at a lower value, there is still an increase in the difference between the converged linewidth andthe l = 0 linewidth. Convergence occurs around l = 12. This shows that as the temperature isdecreased, higher order Dopplerons also play an increasing role.
We next study the effect of optical bistability, and its effects on the optimum input intensity.
As can be seen in Fig. 6.6 (a), when δ0 /γp = 0 there is a bistability in input (|y|2 ) vs intracavity

121

8

ì

ì
à

æ

æ

æ

12

16

24

à

7
à

6

IDΝoptMl 5

ì

4

à
æ

IDΝoptM0

ì
à

ì

æ

ì
à
æ

3
2
1 ì
æ
à
0

2

4

8

l
Figure 6.5: Linewidth at the optimum |y|2 , which gives the smallest linewidth, as a function ofincluded orders of l, normalized by the l = 0 linewidth. Blue circles: N C0 = 600, δ0 /γp = 260,|y|2 = 1270. Red squares: N C0 = 6000, δ0 /γp = 260, |y|2 = 7700. Black diamonds: N C0 = 600,
δ0 /γp = 80, |y|2 = 850.

122(|x|2 ) intensities. As the temperature is increased, the bistability becomes less pronounced, until
eventually it disappears entirely, which can be seen in the δ0 /γp = 30 curve.
The slope at resonance, which can be seen in Fig. 6.6 (b), has very a different behavior forthe zero and non zero temperature cases. For δ0 /γp = 0, the slope is largest on the lower branch,decreases slightly in the bistable region, and then drops exponentially to zero as the intensity is
increased on the upper branch. When the temperature is increased, the bistability disappears, anda dispersion feature with negative slope appears. Then the slope goes to a maximally negativevalue before power broadening eventually causes it to increase to zero.
The disappearance of the bistability is important for achieving the highest possible degree∂φof stabilization. In the zero temperature case, the optimal combination of |x|2 and ∂∆to givethe smallest ∆ν occurs on the far left side of the upper branch, when the input intensity has the
fixed value |y|2 = 4N C0 , labeled with the black dot in Fig. 6.6 (a) and (b). However, this valueof intensity is in the bistable region, so the system is dynamically unstable when the full quantum
dynamics are accounted for. If the tunneling rate between the different branches is not small, weare constrained to work at input intensities that are above the bistable region where the achievablefrequency precision is much less. The grey dot in Fig. 6.6 (a) and (b) is at the far right side of
the bistable region, at the fixed input intensity value |y|2 = (N C0 )2 /4. The stabilization linewidthis orders of magnitude worse here, because the slope is so small. Nonetheless, working above thebistability allows for shot noise limited stabilization line widths of ∼ 1 mHz [63].
Fig. 6.6 (c) shows the stabilization linewidth as a function of N C0 for several temperatures. In
the δ0 /γp = 0 limit of zero temperature, the black dashed curve corresponds to the input intensity
corresponding to the black dot, and the black solid curve corresponds to the input intensity of thegray dot. When the temperature is increased sufficiently, the bistability disappears. Then, there are
no longer any regions of dynamic instability and therefore no restrictions on input intensity. Eachof the δ0 /γp 6= 0 cases in this plot are calculated at their respective optimal input intensity, and stopat a critical value of N C0 , where the bistability reappears. It follows as reciprocal information thatfor a given value of N C0 , there exists a critical temperature at which the bistability disappears.

123
106

(a)

4

10

d0 gp = 0
d0 gp = 30100

d0 gp = 75

ÈxÈ 2

d0 gp = 1501

d0 gp = 300

0.01
10-4
0.01

1

100

ÈyÈ

104

106

2

1.0
0.8
0.6

g p df
2p dD 0.4

g p df
2p dD

0.2

-3

=3 x 10

0.0
-0.2

(b)

0.01

1

100

ÈyÈ

104

106

2

101

10-1

Dn  2 p

10-3

C0 g p 2  gM

10-5

(c)
10-7

102

103

104

105

N C0

Figure 6.6: (a) Scaled input (|x|2 ) vs intracavity (|y|2 ) intensities for several temperatures with
N C0 = 800 (b) Slope at resonance for several temperatures with N C0 = 800 (c) Stabilization
linewidth as a function of N C0 for several temperatures. For δ0 /γp = 0, black dashed is calculated
at the fixed input intensity |y|2 = 4N C0 ; black solid is calculated at the fixed input intensity|y|2 = (N C0 )2 /4. For δ0 /γp 6= 0, the linewidth is calculated with the input intensity fixed to the
optimum value that gives the smallest linewidth. In the case of 88 Sr, δ0 /γp = 150 corresponds to atemperature of ∼ 5 mK.

124
Therefore, by using Fig. 6.6 (c) to note the values of temperature and N C0 at which the bistability
disappears for the different curves, it is seen that the critical temperature increases as N C0 increases.
In general, a lower temperature will yield a smaller linewidth. However, because one must
avoid working in a region of optical bistability, the optimal shot noise limited linewidth for certain
values of N C0 and δ0 /γp in which there is no bistability can actually be smaller than the δ0 /γp = 0
linewidth which is constrained above the bistability.

6.4

Conclusion
We have seen that thermal atoms in a standing wave cavity field exhibit additional phenomenathat were not observed when considering a frozen arrangement of atoms. Specifically, when the
detuning of the laser and atoms is less than the power broadened linewidth, the system interacts
with both components of the standing wave field. This causes new absorptive and dispersive
features in the observables, which are the features that can be used as an error signal for frequency
stabilization. Multi-photon scattering processes due to Doppleron resonances cause the stabilization
linewidth to increase. This effect becomes more dominant as the collective cooperativity N C0 is
increased, and as the temperature is decreased. A system with sufficient N C0 and no atomic motion
will exhibit optical bistability. Atomic motion may cause this bistability to disappear. When
optical bistability occurs, one is generally restricted to using an input intensity that lies outside of
the bistable region in order to prevent hopping between semiclassical solutions. When there is no
optical bistability, there are no dynamically unstable regions, so that no such restrictions on input
intensity are necessary, allowing the optimally smallest stabilization linewidth to be achieved.

Chapter 7

Conclusion

In this thesis, we have examined a variety of open quantum systems in an attempt to discovera scheme in which an extremely spectrally pure coherent light source might be realized. Before
these systems were examined, the necessary preliminary material was provided. To begin, theelectromagnetic field was quantized. This was necessary to correctly describe the Hamiltonian
dynamics of our open quantum systems when only a few photons are present. Also, a quantizedfield is needed in order to describe spontaneous irreversible processes, such spontaneous emission.
The power spectrum, and its relation to the two-time correlation functions through the WienerKhinchin theorem was also discussed. This allowed the spectral properties of a system, such as the
linewidth of the power spectrum, to be calculated.
Next, we derived the master equation, which describes the Schrödinger evolution for a systemthat is weakly coupled to a much larger reservoir. The coupling was treated in time dependent
perturbation theory, and the reservoir was approximated as having infinite bandwidth, and thenaveraged over. Exact solution methods to the master equation, such as the quantum jump andquantum state diffusion methods were presented, and computer simulations of these methods were
demonstrated. The cumulant expansion method, which consists of an approximate set of coupledequations for expectation values of system observables, was also demonstrated.
The quantum Langevin equations, which describe open quantum systems in the Heisenberg
picture of quantum mechanics, were then derived. The same infinite bandwidth approximationthat was made in the master equation derivation was made here. Since these equations prove

126
difficult for computer simulations, a Fokker-Planck equation for a quasi-probability distribution
describing the system was derived. It was shown that a set of c-number Langevin equations areequivalent to the Fokker-Planck equation. These c-number Langevin equations can be simulated on
a computer with relative ease. Once this procedure for deriving the c-number Langevin equations
from the Fokker-Planck equation was demonstrated, an easier, but less mathematically rigorous
way to derive c-number Langevin equations was shown. This new procedure becomes increasingly
advantageous as the system complexity increases. Finally, computer simulations of the c-number
Langevin equations for several systems were demonstrated.
The theory of the crossover between lasing and steady state superradiance was presented. It
was shown how lasing and steady state superradiance can actually be thought of as the extreme
limits of a single and more general phenomenon. In a laser, the phase information that allows forcoherence is stored in the light field. In steady state superradiance, this phase information is stored
in the atomic ensemble. In the crossover region between these two limits, the phase information thatallows for coherence can be stored in the light field and the atomic ensemble simultaneously. The
parameter space between these two extreme limits was explored, and the benefits and drawbacksof operating a system at a given set of parameters, to achieve the most spectrally pure light source,
were discussed.
Specifically, it was shown that for a given rate of re-supplying energy, a crossover system will
be much farther above threshold than a laser, which is typically operated just above threshold.
This can allow the linewidth of a crossover system to be orders of magnitude smaller than the laser,
while maintaining a similar intracavity intensity. A crossover system also exhibits a much larger
intracavity intensity than a superradiant system, so that the crossover system is more experimentally accessible. The linewidth and intracavity intensity were shown to be insensitive to atomic
dephasing in all parameter regions, so long as the rate of atomic dephasing is smaller than the rateof energy being resupplied to the system.
The sensitivity to fluctuations in cavity length and cavity frequency were also discussed.
The superradiance region was shown to be insensitive to fluctuations in cavity frequency, while a127lasing system is insensitive to fluctuations in the atomic frequency. A crossover system was shown
to be less sensitive to fluctuations in cavity frequency than a laser system, and less sensitive tofluctuations in atomic frequency than a superradiant system. This is important, since fluctuations
in cavity length are the bottleneck in further reducing the linewidth in today’s narrowest linewidthlasers.
Finally, we considered the phase stability of a local oscillator locked to a cavity QED systemcomprised of atoms with an ultra-narrow optical transition. The atoms were modeled as being
cooled to millikelvin temperatures and then released into the optical cavity. Even though the atomic
motion introduced Doppler broadening, the standing wave nature of the cavity caused saturated
absorption features to appear, which are much narrower than the Doppler width. These features
can be used to achieve an extremely high degree of phase stabilization, competitive with the current
state-of-the-art. A comparison between the developed theory, and the results of an experiment on
the system that was considered was given. It was demonstrated that the inhomogeneity introducedby finite atomic velocities can cause optical bistability to disappear so that there exist no regionsof dynamic instability that would otherwise restrict operational parameters in the experiment tobe outside the optimum region where the minimum linewidth occurs.

Bibliography[1] G. S. Agarwal and Yifu Zhu. Velocity-tuned resonances and resonance fluorescence in astanding-wave field. Phys. Rev. A, 46:479–484, Jul 1992.
[2] G. S. Agarwal, Yifu Zhu, Daniel J. Gauthier, and T. W. Mossberg. Spectrum of radiation fromtwo-level atoms under intense bichromatic excitation. J. Opt. Soc. Am. B, 8(5):1163–1167,
May 1991.
[3] A. Ashkin, J. M. Dziedzic, Bjorkholm J. E., and S. Chu. Observation of a single-beam gradient
force optical trap for dielectric particles. Optics Letters, 11:288–290, 1986.
[4] Claude Audoin and Bernard Guinot. The measurement of time: time, frequency and theatomic clock. Cambridge University Press, 2001.
[5] Claus Bagger and Flemming O Olsen. Review of laser hybrid welding. Journal of laser
applications, 17(1):2–14, 2005.
[6] Barry C Barish and Rainer Weiss. Ligo and the detection of gravitational waves. Physics
Today, 52:44–50, 1999.
[7] S Blatt, AD Ludlow, GK Campbell, Jan Westenkær Thomsen, T Zelevinsky, MM Boyd, J Ye,
X Baillard, M Fouché, R Le Targat, et al. New limits on coupling of fundamental constants
to gravity using sr 87 optical lattice clocks. Physical Review Letters, 100(14):140801, 2008.
[8] BJ Bloom, TL Nicholson, JR Williams, SL Campbell, M Bishof, X Zhang, W Zhang, SL Bromley, and J Ye. An optical lattice clock with accuracy and stability at the 10-18 level. Nature,506(7486):71–75, 2014.
[9] Justin G Bohnet, Zilong Chen, Joshua M Weiner, Kevin C Cox, and James K Thompson.
Relaxation oscillations, stability, and cavity feedback in a superradiant raman laser. Physicalreview letters, 109(25):253602, 2012.
[10] Justin G. Bohnet, Zilong Chen, Joshua M. Weiner, Dominic Meiser, Murray J. Holland, andJames K. Thompson. A steady-state superradiant laser with less than one intracavity photon.
Nature, 484:78–81, April 2002.
[11] Ruxandra Bondarescu, Mihai Bondarescu, Gyrgy Hetnyi, Lapo Boschi, Philippe Jetzer, andJayashree Balakrishna. Geophysical applicability of atomic clocks: direct continental geoid
mapping. Geophysical Journal International, 191(1):78–82, 2012.

129[12] R. Bonifacio and L. A. Lugiato. Optical bistability and cooperative effects in resonance fluorescence. Phys. Rev. A, 18:1129–1144, Sep 1978.
[13] Jingbiao Chen. Active optical clock. Chinese Science Bulletin, 54(3):348–352, 2009.
[14] JingBiao Chen. Active optical clock. Chinese Science Bulletin, 54:348–352, 2009.
[15] Bjarke TR Christensen, Martin R Henriksen, Stefan A Schäffer, Philip G Westergaard, Jun
Ye, Murray Holland, and Jan W Thomsen. Non-linear spectroscopy of sr atoms in an opticalcavity for laser stabilization. arXiv preprint arXiv:1507.05757, 2015.
[16] G. D. Cole, W. Zhang, M. J. Martin, Jun Ye, and M. Aspelmeyer. Tenfold reduction ofbrownian noise in high-reflectivity optical coatings. Nature Photonics, 7:644, 2013.
[17] Jean Dalibard, Yvan Castin, and Klaus Mølmer. Wave-function approach to dissipative processes in quantum optics. Phys. Rev. Lett., 68:580–583, Feb 1992.
[18] G John Dick. Local oscillator induced instabilities in trapped ion frequency standards. Technical report, DTIC Document, 1987.
[19] R. H. Dicke. Coherence in spontaneous radiation processes. Phys. Rev., 93:99–110, Jan 1954.
[20] SA Diddams, JC Bergquist, SR Jefferts, and CW Oates. Standards of time and frequency atthe outset of the 21st century. Science, 306(5700):1318–1324, 2004.
[21] J.A. Dixon. Surgical application of lasers. 2nd edition. Year Book Medical Publishers,Chicago,
IL, Jan 1987.
[22] Peter Domokos and Helmut Ritsch. Collective cooling and self-organization of atoms in acavity. Phys. Rev. Lett., 89:253003, Dec 2002.
[23] RWP Drever, John L Hall, FV Kowalski, J Hough, GM Ford, AJ Munley, and H Ward. Laser
phase and frequency stabilization using an optical resonator. Applied Physics B, 31(2):97–105,
1983.
[24] P D Drummond and D F Walls. Quantum theory of optical bistability. i. nonlinear polarisability model. Journal of Physics A: Mathematical and General, 13(2):725, 1980.
[25] FJ Duarte. Broadly tunable external-cavity semiconductor lasers. Tunable Laser Applications(CRC, 2009), pages 143–178, 2008.
[26] F.J.E.K.C. Duarte. Two-laser therapy and diagnosis device, September 28 1988. EP PatentApp. EP19,880,302,468.
[27] R. Dum, P. Zoller, and H. Ritsch. Monte carlo simulation of the atomic master equation for
spontaneous emission. Phys. Rev. A, 45:4879–4887, Apr 1992.
[28] L Essen. Atomic time and the definition of the second. Nature, 1956.
[29] L Essen and JVL Parry. An atomic standard of frequency and time interval: a caesium
resonator. Nature, 1955.

130[30] Louis Essen and JVL Parry. The caesium resonator as a standard of frequency and time.
Philosophical Transactions of the Royal Society of London A: Mathematical, Physical andEngineering Sciences, 250(973):45–69, 1957.
[31] Z. Ficek and H. S. Freedhoff. Resonance-fluorescence and absorption spectra of a two-level
atom driven by a strong bichromatic field. Phys. Rev. A, 48:3092–3104, Oct 1993.
[32] Charles Freed and Ali Javan. Standing-wave saturation resonance in the co2 10.6-µ transitions
observed in a low-pressure room-temperature absorber gas. Applied Physics Letters, 17(2):53–
56, 1970.
[33] James G Fujimoto, Mark E Brezinski, Guillermo J Tearney, Stephen A Boppart, Brett Bouma,
Michael R Hee, James F Southern, and Eric A Swanson. Optical biopsy and imaging using
optical coherence tomography. Nature medicine, 1(9):970–972, 1995.
[34] Crispin W Gardiner et al. Handbook of stochastic methods, volume 4. Springer Berlin, 1985.
[35] R. M. Godun, P. B. R. Nisbet-Jones, J. M. Jones, S. A. King, L. A. M. Johnson, H. S. Margolis,
K. Szymaniec, S. N. Lea, K. Bongs, and P. Gill. Frequency ratio of two optical clock transitionsin 171 yb+ and constraints on the time variation of fundamental constants. Phys. Rev. Lett.,
113:210801, Nov 2014.
[36] Leon Goldman. Dye lasers in medicine. In F.J. DuarteLloyd W. Hillman, editor, Dye Laser
Principles, Quantum ElectronicsPrinciples and Applications, pages 419 – 432. Academic Press,
San Diego, 1990.
[37] J. Gripp, S. L. Mielke, and L. A. Orozco. Evolution of the vacuum rabi peaks in a detuned
atom-cavity system. Phys. Rev. A, 56:3262–3273, Oct 1997.
[38] J. Guéna, M. Abgrall, D. Rovera, P. Rosenbusch, M. E. Tobar, Ph. Laurent, A. Clairon, and
S. Bize. Improved tests of local position invariance using 87 Rb and 133 Cs fountains. Phys. Rev.
Lett., 109:080801, Aug 2012.
[39] H. Haken. Theory of intensity and phase fluctuations of a homogeneously broadened laser. Z.
Phys., 190:327, 1966.
[40] Hermann Haken. Laser Theory. Springer Berlin Heidelberg, 1983.
[41] John L Hall. Nobel lecture: Defining and measuring optical frequencies. Reviews of Modern
Physics, 78(4):1279, 2006.
[42] TW Hansch and B Couillaud. Laser frequency stabilization by polarization spectroscopy of areflecting reference cavity. Optics communications, 35(3):441–444, 1980.
[43] Stephan Hartmann. Generalized dicke states. arXiv, 1201.1732, 2012.
[44] N Hinkley, JA Sherman, NB Phillips, M Schioppo, ND Lemke, K Beloy, M Pizzocaro,
CW Oates, and AD Ludlow. An atomic clock with 10–18 instability. Science, 341(6151):1215–
1218, 2013.
[45] N. Huntemann, B. Lipphardt, Chr. Tamm, V. Gerginov, S. Weyers, and E. Peik. Improved
limit on a temporal variation of mp /me from comparisons of yb+ and cs atomic clocks. Phys.
Rev. Lett., 113:210802, Nov 2014.

131[46] Steve Jefferts. Nist-f1 cesium fountain atomic clock, 2009.
div688/grp50/primary-frequency-standards.cfm.

http://www.nist.gov/pml/[47] James Jespersen and Jane Fitz-Randolph. From sundials to atomic clocks: understanding timeand frequency. Courier Corporation, 1999.
[48] YY Jiang, AD Ludlow, ND Lemke, RW Fox, JA Sherman, L-S Ma, and CW Oates. Making optical atomic clocks more stable with 10-16-level laser stabilization. Nature Photonics,
5(3):158–161, 2011.
[49] T. Kessler, C. Hagemann, C. Grebing, T. Legero, U. Sterr, F. Riehle, M. J. Martin, L. Chen,
and J. Ye. A sub-40-mhz-linewidth laser based on a silicon single-crystal optical cavity. Nat.
Photon, 6:687–692, Sep 2012.
[50] T Kessler, C Hagemann, C Grebing, T Legero, U Sterr, F Riehle, MJ Martin, L Chen, andJ Ye. A sub-40-mhz-linewidth laser based on a silicon single-crystal optical cavity. Nature
Photonics, 6(10):687–692, 2012.
[51] H. J. Kimble, Benjamin L. Lev, and Jun Ye. Optical interferometers with reduced sensitivity
to thermal noise. Phys. Rev. Lett., 101:260602, Dec 2008.
[52] P.E. Kloeden and E. Platen. Numerical Solution of Stochastic Differential Equations. Stochastic Modelling and Applied Probability. Springer Berlin Heidelberg, 2011.
[53] R Kodama, PA Norreys, K Mima, AE Dangor, RG Evans, H Fujita, Y Kitagawa, K Krushelnick, T Miyakoshi, N Miyanaga, et al. Fast heating of ultrahigh-density plasma as a step
towards laser fusion ignition. Nature, 412(6849):798–802, 2001.
[54] R Kodama, H Shiraga, K Shigemori, Y Toyama, S Fujioka, H Azechi, H Fujita, H Habara,
T Hall, Y Izawa, et al. Nuclear fusion: Fast heating scalable to laser fusion ignition. Nature,418(6901):933–934, 2002.
[55] M. I. Kolobov, L. Davidovich, E. Giacobino, and C. Fabre. Role of pumping statistics anddynamics of atomic polarization in quantum fluctuations of laser sources. Phys. Rev. A,
47:1431–1446, Feb 1993.
[56] Ryogo Kubo. Generalized cumulant expansion method. Journal of the Physical Society ofJapan, 17(7):1100–1120, 1962.
[57] S. J. M. Kuppens, M. P. van Exter, and J. P. Woerdman. Quantum-limited linewidth of abad-cavity laser. Phys. Rev. Lett., 72:3815–3818, Jun 1994.
[58] R. Le Targat, L. Lorini, Y. Le Coq, M. Zawada, J. Guna, M. Abgrall, M. Gurov, P. Rosenbusch, D. G. Rovera, B. Nagrny, R. Gartman, P. G. Westergaard, M. E. Tobar, M. Lours,
G. Santarelli, A. Clairon, S. Bize, P. Laurent, P. Lemonde, and J. Lodewyck. Experimental
realization of an optical second with strontium lattice clocks. Nat. Commun., 4(2109), Jul
2013.
[59] Marc Levenson. Introduction to Nonlinear Laser Spectroscopy. Elsevier, 2012.
[60] TH MAIMAN. Stimulated optical radiation in ruby. Nature, 1960.

132[61] Eckhard Mandelkow and Eva-Maria Mandelkow. Kinesin motors and disease. Trends in cellbiology, 12(12):585–591, 2002.
[62] M. J. Martin, M. Bishof, M. D. Swallows, X. Zhang, C. Benko, J. von Stecher, A. V. Gorshkov,
A. M. Rey, and J. Ye. A quantum many-body spin system in an optical lattice clock. Science,
341(6146):632–636, 2013.
[63] M. J. Martin, D. Meiser, J. W. Thomsen, Jun Ye, and M. J. Holland. Extreme nonlinear
response of ultranarrow optical transitions in cavity qed for laser stabilization. Phys. Rev. A,
84:063813, Dec 2011.
[64] D. Meiser and M. J. Holland. Intensity fluctuations in steady-state superradiance. Phys. Rev.
A, 81:063827, Jun 2010.
[65] D. Meiser and M. J. Holland. Steady-state superradiance with alkaline-earth-metal atoms.
Phys. Rev. A, 81:033847, Mar 2010.
[66] D. Meiser, Jun Ye, D. R. Carlson, and M. J. Holland. Prospects for a millihertz-linewidthlaser. Phys. Rev. Lett., 102:163601, Apr 2009.
[67] Pierre Meystre and Murray Sargent. Elements of quantum optics. Springer Science & BusinessMedia, 2007.
[68] Kenji Numata, Amy Kemery, and Jordan Camp. Thermal-noise limit in the frequency stabilization of lasers with rigid cavities. Physical review letters, 93(25):250602, 2004.
[69] M. B. Plenio and P. L. Knight. The quantum-jump approach to dissipative dynamics inquantum optics. Rev. Mod. Phys., 70:101–144, Jan 1998.
[70] William H. Press, Saul A. Teukolsky, William T. Vetterling, and Brian P. Flannery. Numerical
Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press, New York,NY, USA, 3 edition, 2007.
[71] Audrey Quessada, Richard P Kovacich, Irène Courtillot, André Clairon, Giorgio Santarelli,and Pierre Lemonde. The dick effect for an optical frequency standard. Journal of Optics B:
Quantum and Semiclassical Optics, 5(2):S150, 2003.
[72] Norman F Ramsey. History of atomic clocks. J. Res. Nat. Bur. Stand, 88:301, 1983.
[73] A. L. Schawlow and C. H. Townes. Infrared and optical masers. Phys. Rev., 112:1940–1949,
Dec 1958.
[74] S. Schiller, G.M. Tino, P. Gill, C. Salomon, U. Sterr, E. Peik, A. Nevsky, A. Grlitz, D. Svehla,
G. Ferrari, N. Poli, L. Lusanna, H. Klein, H. Margolis, P. Lemonde, P. Laurent, G. Santarelli,
A. Clairon, W. Ertmer, E. Rasel, J. Mller, L. Iorio, C. Lmmerzahl, H. Dittus, E. Gill,
M. Rothacher, F. Flechner, U. Schreiber, V. Flambaum, Wei-Tou Ni, Liang Liu, Xuzong
Chen, Jingbiao Chen, Kelin Gao, L. Cacciapuoti, R. Holzwarth, M.P. He, and W. Schfer.
Einstein Gravity Explorera medium-class fundamental physics mission, volume 23. Springer
Netherlands, 2009.
[75] Marlan O. Scully and M. Suhail Zubairy. Quantum Optics. Cambridge University Press, 1
edition, September 1997.

133[76] Stig Stenholm. Foundations of Laser Spectroscopy. Dover Publications, 2005.
[77] Stig Stenholm and Willis E Lamb Jr. Semiclassical theory of a high-intensity laser. Applied
Physics Letters, 181(2):618, 1969.
[78] Karel Svoboda, Christoph F. Schmidt, Bruce J. Schnapp, and Steven M. Block. Direct observation of kinesin stepping by optical trapping interferometry motility. Nature, 365(6448):721–727,
1993.
[79] Andrée Tallet. Raman emission versus dopplerons in a two-level atomic cell. J. Opt. Soc. Am.
B, 11(8):1336–1349, Aug 1994.
[80] D. A. Tieri, J. Cooper, Bjarke T. R. Christensen, J. W. Thomsen, and M. J. Holland. Laserstabilization using saturated absorption in a cavity-qed system. Phys. Rev. A, 92:013817, Jul
2015.
[81] D. A. Tieri, Minghui Xu, D. Meiser, J. Cooper, and M. J. Holland. The crossover between
lasing and steady state superradiance, 2015. In preperation.
[82] Charles H. Townes. How the laser happened : adventures of a scientist / Charles H. Townes.
Oxford University Press New York, 1999.
[83] Th Udem, Ronald Holzwarth, and Theodor W Hänsch. Optical frequency metrology. Nature,416(6877):233–237, 2002.
[84] M. P. van Exter, S. J. M. Kuppens, and J. P. Woerdman. Theory for the linewidth of abad-cavity laser. Phys. Rev. A, 51:809–816, Jan 1995.
[85] Daniel F Walls and Gerard J Milburn. Quantum optics. Springer Science & Business Media,
2007.
[86] Philip G. Westergaard, Bjarke T. R. Christensen, David Tieri, Rastin Matin, John Cooper,
Murray Holland, Jun Ye, and Jan W. Thomsen. Observation of motion-dependent nonlinear
dispersion with narrow-linewidth atoms in an optical cavity. Phys. Rev. Lett., 114:093002,
Mar 2015.
[87] H. M. Wiseman and G. J. Milburn. Quantum theory of field-quadrature measurements. Phys.
Rev. A, 47:642–662, Jan 1993.
[88] Minghui Xu, D. A. Tieri, and M. J. Holland. Simulating open quantum systems by applyingsu(4) to quantum master equations. Phys. Rev. A, 87:062101, Jun 2013.
[89] Minghui Xu, DA Tieri, EC Fine, James K Thompson, and MJ Holland. Synchronization of
two ensembles of atoms. Physical review letters, 113(15):154101, 2014.
[90] Minghui Xu, DA Tieri, and MJ Holland. Simulating open quantum systems by applying su(4) to quantum master equations. Physical Review A, 87(6):062101, 2013.
[91] B. C. Young, F. C. Cruz, W. M. Itano, and J. C. Bergquist. Visible lasers with subhertz
linewidths. Phys. Rev. Lett., 82:3799–3802, May 1999.