Stochastic Models of Atom-Photon Dynamics with
Applications to Cooling Quantum Gases
byJosh W. Dunn
B.A., University of Colorado, Boulder, 2002

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
2007

This thesis entitled:
Stochastic Models of Atom-Photon Dynamics with Applications to Cooling Quantum
Gases
written by Josh W. Dunn
has been approved for the Department of Physics

Chris H. Greene

John L. Bohn

Date

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

iii
Dunn, Josh W. (Ph.D., Physics)
Stochastic Models of Atom-Photon Dynamics with Applications to Cooling Quantum
Gases
Thesis directed by Professor Chris H. Greene

Through the years, stochastic physics has provided important insight into natural phenomena that possess an inherently random nature. From its foundations in the
study of Brownian motion up through its myriad present applications, a stochastic description of nature has yielded an elegant theoretical understanding, as well as providing
practical and efficient simulation techniques. One particularly important application isto the field of quantum optics, in which the interaction of light and matter is treatedin a fundamentally quantum-mechanical manner. The work presented here utilizes themethods of stochastic physics to understand a variety of quantum-optical phenomena
involving the dynamics of atoms interacting with photons. A solid theoretical understanding of such phenomena is often necessary to describe laser cooling of atoms, andmany such applications are discussed here. A detailed model of atoms with complex
internal structure interacting with three-dimensional laser fields is presented, as wellas the rich dynamics of three-level atoms interacting with two lasers. Applications tocavity cooling of atoms and molecules are discussed, and a method for describing nonMarkovian dynamics using relaxation techniques is presented. A novel cooling scheme
utilizing Feshbach resonances in the scattering of two atoms is also treated.

Acknowledgements

I would like to thank Chris Greene for his many years of guidance and support
during my graduate career. His passion for physics is contagious and it inspires those
who interact with him strive for a deeper understanding of nature. His curiosity extends
over a very broad range, and he is often willing to tackle projects outside of his main
expertise in order to expand his understanding of physics. I appreciate his willingness
to allow me to pursue a study of quantum optics and stochastic physics even given his
limited experience in these fields. I believe that our collaboration has been very fruitful.
Combining my day-to-day learning about quantum optics with Chris’s comprehensive
understanding of traditional atomic physics has allowed us to develop models of atomphoton dynamics with unprecedented attention to an accurate treatment of atomic
internal degrees of freedom. Our collaboration has often led to the realization of parallels
between quantum-optical concepts and the many theories developed by Fano, which
Chris knows so well. For me, this conceptual bridge has been very useful, illuminating
tough subjects and often revealing an elaborate and scattered collection of methodologies
to be unified under a simple physical picture. I am proud to be a part of this scientificlineage.
I am grateful to the many physicists at JILA who helped me during my time here.
Jun Ye has been kind enough to provide me with many opportunities to put theory
to use in describing various aspects of his vast array of cold-atom and cold-molecule
experiments. Over the years, Jinx Cooper has been a helpful and reliable fountain of

v
physics knowledge. He was always willing to utilize his lifetime of experience to help
me understand the tricky and subtle aspects of quantum optics. Robin Santra, during
his time as a post doc in Chris’s group, was patient and willing to answer my manyquestions involving atomic structure.
I would also like to thank various other collaborators that made my work possible.
Jan Thomsen from Denmark and Flavio Cruz from Brazil provided the experimental
background and the impetus for our joint project that led to a better understanding ofmultilevel atomic dynamics as well as a cooling scheme that they are currently implementing. Ben Lev’s strong knowledge of cavity QED led to many interesting discussions
between us about the dynamics of atoms and molecules in an optical cavity.
Finally, I would like to thank the National Science Foundation for funding myresearch as a graduate student.

viContents

Chapter
1 Introduction

1

2 Feshbach-Resonance Cooling8

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2

Collision-theory basics . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

2.2.1

Potential scattering theory . . . . . . . . . . . . . . . . . . . . .

9

2.2.2

Cold collisions . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

2.2.3

Feshbach resonances . . . . . . . . . . . . . . . . . . . . . . . . .

13

Zero-range-potential model for two particles in a harmonic trap . . . . .

14

2.3.1

One channel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

2.3.2

Two channels . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.4

Square-well-potential model for two particles in free space . . . . . . . .

19

2.5

Digression: application to molecular dissociation . . . . . . . . . . . . .

20

2.5.1

Zero-range potential . . . . . . . . . . . . . . . . . . . . . . . . .

21

2.5.2

Square-well potential . . . . . . . . . . . . . . . . . . . . . . . . .

28

2.6

The cooling scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30

2.7

Results and experimental possibilities . . . . . . . . . . . . . . . . . . .

37

2.8

Conclusion

46

2.3

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vii
3 Dynamics of Multilevel Atoms in Three-Dimensional Polarization-Gradient Fields 47
3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

47

3.2

Master equation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

50

3.2.1

General form of the master equation . . . . . . . . . . . . . . . .

52

3.2.2

Master equation in the low-intensity limit: adiabatic eliminationof excited states . . . . . . . . . . . . . . . . . . . . . . . . . . .

56

3.3

Direct solutions of the semiclassical master equation . . . . . . . . . . .

57

3.4

Stochastic wave-function solutions of the fully quantum master equation

59

3.5

Calculations for 25 Mg and 87 Sr . . . . . . . . . . . . . . . . . . . . . . .

64

3.6

Conclusion

71. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4 Bichromatic Cooling of Three-Level Atoms

73

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

73

4.2

Three-level atomic systems . . . . . . . . . . . . . . . . . . . . . . . . .

75

4.3

Fully quantum model and its advantages . . . . . . . . . . . . . . . . . .

78

4.4

Numerical solution . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

78

4.5

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

79

4.6

EIT explanation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

80

4.7

Atom-photon dynamics for asymmetric lineshapes . . . . . . . . . . . .

83

4.8

Conclusion

86. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5 Cavity Cooling of Atoms and Molecules

87

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

87

5.2

Physical picture of the cooling mechanism . . . . . . . . . . . . . . . . .

89

5.3

Theoretical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

91

5.3.1

Single atom . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

92

5.3.2

Multiple atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

5.3.3

Accounting for Raman losses and other experimental effects . . .

96

viii
5.4

Comparison of fully quantum to semiclassical model . . . . . . . . . . .

5.5Conclusion

97. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101

6 Nonlinear Optical Spectroscopy

103

6.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

6.2

Basics of nonlinear optical spectroscopy . . . . . . . . . . . . . . . . . . 105

6.3

Phenomenological damping model . . . . . . . . . . . . . . . . . . . . . 107

6.4

Transient multiwave mixing using relaxation theory . . . . . . . . . . . . 109

6.5

6.4.1

Relaxation theory of static pressure broadening . . . . . . . . . . 109

6.4.2

Extension to transient nonlinear-optics problems . . . . . . . . . 111

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

Bibliography

113

Appendix
A Liouville Space121

B Spontaneous-Emission Relaxation Operator for Atoms with Complex Internal
Structure
C Semiclassical Master Equation for Atoms with Complex Internal Structure

124
138

C.1 Equations without spontaneous decay . . . . . . . . . . . . . . . . . . . 138
C.1.1 Atom-field interaction . . . . . . . . . . . . . . . . . . . . . . . . 138
C.1.2 Transformation to moving, rotating frame . . . . . . . . . . . . . 143
C.2 Optical Bloch equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 147
C.2.1 Class I: optical coherences . . . . . . . . . . . . . . . . . . . . . . 148
C.2.2 Class II: ground-state populations and ground-state coherences . 149
C.2.3 Class III: excited-state populations and excited-state coherences

150

ix

TablesTable
4.1

Transition linewidths, Γi = γi /2π for i =1, 2, of the 1 S0 -1 P1 -1 S0 threelevel Ξ systems for common alkaline-earth atoms 24 Mg, 40 Ca, and 88 Sr,
as well as the alkaline-earth-like atom 174 Yb. The ratio of Γ1 /Γ2 is also
shown, indicating the factor below the lower-transition Doppler limit
which might be expected using a three-level bichromatic scheme, neglecting the more detailed lineshape analysis provided in this chapter. . . . .

76

xFiguresFigure2.1

Eigenspectrum in oscillator units ~ω for two particles in a harmonic
trap, as a function of the two-particle scattering length a in units of
the reduced-mass oscillator length Losc . . . . . . . . . . . . . . . . . . .

2.2

17

An example of the type of magnetic field pulse used in molecular-dissociation
experiments. This pulsed magnetic field is used here as a perturbation
in a in a two-channel model of dissociation. The amplitude of the pulse
is Bini , the width of the envelope is tpert , and the frequency of oscillationis νpert . In this plot, νpert = 10/tpert . . . . . . . . . . . . . . . . . . . . .

2.3

24

An example of a typical molecular photodissociation spectrum obtainedusing the two-channel square-well model. Plotted is the photodissociation
rate as a function of energy, both in arbitrary units. . . . . . . . . . . .

2.4

29

Eigenenergies of an interacting two-atom system in a harmonic trap neara Feshbach resonance, as a function of the magnetic field. The shift
in the energies as the resonance is traversed can be seen, as well as theenergy dependence of the location of the shift in magnetic field. The inset
illustrates an ideal series of population transfers induced by magnetic fieldramps for the Feshbach-resonance cooling scheme. This figure is adaptedfrom Ref. [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

xi
2.5

Upper figure: Wave function in a box before and after rapidly (diabatically) expanding the box. Lower plot: Wave function in a box before andafter slowly (adiabatically) expanding the box. . . . . . . . . . . . . . .

2.6

36

Illustration of the effect of a single Feshbach resonance cooling cycle for
an atom pair taken from a thermal distribution with T = 1 mK andtrapped with frequency ν = 1 MHz. The black (red) line represents the
population distribution before (after) application of one slow and one fastmagnetic field ramp. The state Q (here Q = 15) is indicated. This figure
is adapted from Ref. [1]. . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.7

40

Diagram illustrating the effect of the Feshbach-resonance cooling schemeon an s-wave distribution of atoms. Shown is the population in arbitrary
units as a function of the relative-coordinate quantum number n. Population from the state n = Q is increased in energy, and assumed to be
removed from the trap. All population with n > Q is decreased in energy
by one relative-coordinate oscillator unit of energy. . . . . . . . . . . . .

2.8

42

Probability that a pair of atoms remains trapped vs. the average total
kinetic energy of the two atoms in oscillator units (note that kB T =hEtot i/6 for two harmonically trapped atoms). Three different coolingparameters are used: 2~ωQ = 5τ (solid line), 9τ (dashed line), and 12τ
(dot-dashed line). It is assumed that rethermalization occurs between
cooling cycles (see text), although this scheme does not necessarily requireit. Inset: probability to remain trapped vs. the number of cooling cycles
for the same three cooling parameters. This figure is adapted from Ref. [1]. 45

xii
3.1

Energy level diagram of an atom with multiple hyperfine manifolds. If theenergy spacing of the excited-state manifolds are of the order or smaller
than the natural linewidth of the transition, the usual sub-Doppler cooling transition (Fg ↔ Fe = Fg + 1) is not isolated and the other manifolds
must be taken into account. This figure is adapted from Ref. [2]. . . . .

3.2

51

An example of a characteristic MCWF stochastic trajectory. Shown isthe average of the atomic center-of-mass energy over 500 independent
stochastic wave functions, as a function of time, for a two-level atom in
a 1D standing-wave laser field. The energy is given in units of the recoil
energy Erec = ~2 k 2 /2m, and time is given in units of the inverse recoil
frequency ωr−1 = ~/Erec . All wave functions are initialized in the groundstate of the atom and localized in momentum space with zero momentum.
The steady state, wherein the system fluctuates around an average value,
is seen to be achieved after a transient relaxation period. Error bars
indicate the variance in the data at each given time for the ensemble of
500 stochastic wave functions. An estimate of the steady-state atomic
center-of-mass energy is obtained by performing a time-average over all
wave functions for all times after the relaxation regime. The error bar ofsuch an average will be smaller than the error bars in the figure, which
apply only to the data for a given time. This figure is adapted from Ref. [2]. 63

xiii
3.3

Simulation of a two-level atom with δ = −Γ/2 and Γ = 400 ωrec interacting with a one-dimensional laser field, as a test of the Monte Carlowave-function technique. Plotted is kB multiplied by the steady-state
temperature, in units of the recoil energy ~ωrec . The dots are the numerical results of the Monte Carlo simulation, and the line indicates aquadratic least-squares fit the data. The calculation yields a low-intensity
limit for the one-dimensional Doppler temperature in agreement with the7
known value TD = 40
~Γ/kB = 70 ~ωrec /kB . . . . . . . . . . . . . . . . .

3.4

65

Results for calculated ensemble-average energies (rms momentum squared)for 25 Mg and 87 Sr, as a function of the light-shift parameter ~|δ3 |s3 /(2Erec ).
For comparison, also shown are the calculated energies for atoms with isolated transitions, Je = Jg + 1, with Jg =1, 2, 3, and 4, with detuningδ = −5Γ. This figure is adapted from Ref. [2]. . . . . . . . . . . . . . . .

3.5

69

Semiclassical force curves illustrating the effect of increasing degeneracy
on the laser cooling of idealized atoms with isolated cooling transitionand no overlap of hyperfine manifold. As the internal atomic degeneracy
is increased and all other parameters held fixed, the appearance of a
sub-Doppler force grows. See text for details. . . . . . . . . . . . . . . .

3.6

70

Three-dimensional renderings illustrating the time evolution of an atomic
cloud in momentum space, based on the numerical simulations of cooling
87 Sr.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

72

xiv
4.1

Left side: Atomic configuration for generic three-level Ξ system. Rightside: The dressed eigenenergies for the 1 S0 -1 P1 -1 S0 Ξ system in 24 Mg
with s1 (δ1 ) = 0.001 and s2 (δ2 ) = 1. Real (top) and imaginary (bottom)
parts of the eigenvalues of Eq. (4.3), with dressed atomic states labeled.
The real parts are the energies and the imaginary parts are the effective
linewidths of the dressed atomic system. Both are plotted as functionsof δ2 , with fixed δ1 = 0. This figure is adapted from Ref. [3]. . . . . . . .

4.2

77

Steady-state temperatures for bichromatic three-level laser cooling, as afunction of both detunings δ1 and δ2 , with both atomic transitions perturbatively probed, s1 ¿ 1 and s2 ¿ 1. The bare two-photon resonanceis indicated with a dashed line. . . . . . . . . . . . . . . . . . . . . . . .

4.3

81

Steady-state temperatures for bichromatic three-level laser cooling, as afunction of both detunings δ1 and δ2 , with the lower atomic transition
perturbatively probed, s1 ¿ 1, for both plots, and the upper atomic
transition dressed to varying degrees in the two plots, with s2 = 1 in theupper plot and s2 = 5 in the lower plot. The bare two-photon resonanceis indicated with a dashed line. This figure is adapted from Ref. [3]. . .

4.4

82

Upper plot: Comparison of the true absorption spectrum of the dressed
three-level Ξ system (solid line) with a simplistic absorption spectrumwith Lorentzian lineshapes (dotted line). Optimum laser-cooling detuning is indicated with arrows for each type of spectrum. Lower plot: Comparison of the numerical three-level-cooling temperature results (datapoints) with the ratio of the maximum slope of a Lorentzian lineshapewith width Γ1 to the slope of the asymmetric lineshape, as a function of
δ1 with δ2 = −Γ1 /2 (solid line). This ratio provides an indication of the
expected cooling for the dressed system relative to the Doppler limit for
the lower transition. This figure is adapted from Ref. [3].

. . . . . . . .

85

xv
5.1

Time evolution of the momentum probability distribution for an ensembleof particles undergoing cavity laser cooling. Time is indicated, measured
in units of the inverse recoil frequency ωr−1 . . . . . . . . . . . . . . . . .

5.2

99

Time evolution of the average kinetic energy of an ensemble of particles
undergoing cavity laser cooling. Energy is given in units of the recoil
energy ~ωr and time is measured in units of the inverse recoil frequency
ωr−1 . The black line was obtained using the fully quantum model while
the red line was obtained using the semiclassical model. . . . . . . . . . 100

5.3Time evolution of the cavity-mode population. The black line was obtained using the fully quantum model, calculating the quantum-mechanicalaverage of the cavity photon number operator n = a† a. The red line wasobtained using the semiclassical model, calculating the average of thesquare of the cavity-mode amplitude α2 . Time is measured in units of
the inverse recoil frequency ωr−1 . . . . . . . . . . . . . . . . . . . . . . . 102

B.1 The partitioning of the density operator for an atom with multiple coupled excited-state manifolds, each potentially having multiple substates.

128

Chapter 1
Introduction

The beginnings of the stochastic description of nature can be traced back to Einstein’s 1905 work [4] concerning the observation of Brownian motion, in which small
particles immersed in water follow erratic, seemingly random trajectories. Einstein’s
analysis, performed concurrently and independently by Smoluchowski [5], seems quite
simple, yet it contains many of the elements of modern stochastic physics such as an
idea later called the Markov approximation and a precursor to the Fokker-Planck equation [6]. Einstein’s idea, that a physical theory could contain aspects that are random
at a fundamental level, was revolutionary at the time, going beyond the traditional deterministic equations with which physicists had been primarily concerned prior to that.
His equations employed diffusion physics, the nature of which could be determined by
the microscopic interaction properties of the observed particle with atoms in the liquid.
Inspired by Einstein’s theory, Langevin [7] later developed an equivalent theory which,instead of a diffusion term, utilized the idea of a random fluctuating force. Langevin
was able to work out the basic properties of such forces, and equations of this type bear
his name today. Through the years, many physicists have built on these foundations to
describe a large spectrum of stochastic problems. The work presented in this thesis continues humbly along these lines, modeling various physical phenomena using stochastic
equations of motion and hopefully adding useful knowledge to the canon.
Quantum optics, another venerable branch of physics, explores the quantum-

2mechanical nature of the interaction of light and matter. This encompasses a diverse
and rich set of phenomena, from the nature of spontaneous emission to the scatteringof intense, coherent electromagnetic waves provided by a laser to the use of an opticalcavity to modify the properties of the photon vacuum as observed by an atom. The
field of quantum optics contains its own set of unique concepts, techniques, and tools.
Quantum opticians like to reduce the complex internal structure of an atom down to asimple two-level system whenever possible. This system can be mapped onto a spin-1/2
system for which powerful conceptual tools like the Bloch sphere can be employed. Laserinteractions with the system can then be thought of as rotating a vector in a certain
way around the Bloch sphere. Quantum optics also often intersects with the field of
stochastic physics. Most of the work in this thesis involves light-matter interactions andthus can be classified under the heading of quantum optics.
Contained in the example of Brownian motion are the elements of the theoryof system-reservoir interactions [8]. The particle’s erratic motion in the liquid can bequalitatively separated into two unique influences. Over a long time scale, one type offorce acts to gradually change the course of the particle’s motion. On a much shortertime scale, another type of force induces the sharp motions that give the trajectory
its random appearance. Treating the particle’s interaction with individual atoms in
the liquid as interaction with a reservoir, a detailed equation of motion can be written
down. Formal solutions of this equation can be manipulated by viewing progressionof time to be in discrete steps of length greater than the short time scale but shorter
than the longer time scale, a procedure known as coarse-graining. In this manner, theeffects of the frequently occurring collisions with atoms in the liquid can be wrapped
up into an average effect. Such an equation, termed a master equation, is useful because it removes the difficulty of explicitly treating the overwhelmingly large number ofmicroscopic processes that contribute to the particle’s overall motion.
Generally speaking, a master equation is obtained from the differential Chapman-

3
Kolmogorov equation by setting the drift and diffusion terms to zero. The master equation describes evolution of the probability density with stochastic jumps that can occur
at random intervals, in addition to any deterministic evolution. In the context of quantum optics, the master equation is useful for describing a system of one or more atoms,
perhaps coupled to a laser, interacting with the ever-present vacuum photon reservoir.
The system undergoes a deterministic evolution due to its own internal interactions, aswell as a stochastic evolution manifested by jumps that occur due to interactions withthe vacuum photon field. The quantum-optical master equation is useful because its
solutions yield the density operator (often called, less generally, the density matrix) for
the system at any time. The density operator can then be used to obtain average system properties by taking its product with quantum-mechanical operators and tracing
over the system degrees of freedom. In this thesis, quantum-optical systems are almost
exclusively modeled using master equations, as in Chapters 3 through 6.
Another aspect of the utility of the quantum-optical master equation is its amenability to stochastic Monte Carlo techniques [9, 10, 11, 12, 13, 14, 15]. The system-reservoirinteraction portion of the master equation can be unraveled and interpreted as a sum
over an infinite number of individual trajectories for wave functions rather than density matrices. These trajectories each execute a random walk in phase space withrandomly spaced jumps occurring due to the interaction with vacuum photon states.
Besides offering an elegant fundamental interpretation of system-reservoir interactions,this unraveling can be put to use as a simulation technique. A set of wave functions
can be independently propagated with jump processes occurring according to a numerical pseudo-random number generator. Although the exact master-equation evolution is
only obtained in the limit of an infinite number of such trajectories, a rather small number of trajectories can provide a very good approximation of the dynamics. Moreover,the error involved in using a finite number instead of an infinite number of trajectories
is implicitly determined by the technique. For systems with a large number of degrees

4of freedom N , the propagation of wave functions of size N rather than a density matrixof size N 2 offers extensive numerical advantages, both in terms of computational speed
and storage. This actually understates the problem, since the master equation for adensity operator cannot be phrased as a linear N × N matrix equation. Only by goingto a space of N 2 × N 2 can the problem be expressed in a form for which a numericallinear matrix solver can be used. In this thesis, stochastic-trajectory techniques are
utilized in Chapters 3 and 5. This subject is formalized by the concept of Liouvillespace, which is discussed in detail in Chapters 3 and 4, as well as Appendix A.
A good theoretical understanding of cooling processes is of practical use as well
as being fundamentally intriguing. Cooling is basically a means of decreasing a system’s
footprint in phase space, and as such, it makes use of many ideas in thermodynamic
theory. Ideas for cooling schemes are continually surfacing, and it can be an enlightening
exercise to test them against the rigors of the second law of thermodynamics. Once suchidea is presented here in Chapter 2. From a practical standpoint, cooling is tremendously
important in the today’s field of atomic, molecular, and optical physics. The variety ofquantum phenomena observed and employed in current experiments is made possible
primarily when the system, often a gas of atoms, is cooled to very low temperatures.
Such temperature regimes reduce the noise associated with higher temperature systems,
as in experiments involving atomic clocks and precision standards. In addition, they
can allow for quantum degeneracy to be obtained, as in the Bose-Einstein condensatesand degenerate Fermi gases common to experiments now.
Of all the various cooling methods, one particular type — laser cooling — hasbeen perhaps the most influential. The idea of scattering detuned laser light off anensemble of atoms to reduce their kinetic energy has contributed in some way to justabout every experiment involving cold atoms. A thorough understanding of the physics
of laser cooling [16, 17, 18], developed throughout the 1980’s and 1990’s, paved the way
to a Nobel Prize for its major contributors in 1997. Laser cooling comes in a variety of

5forms and complexities, and Chapters 3, 4, and 5 discuss the theoretical modeling of afew of these.
The main themes of the work presented here are stochastic quantum mechanics
and cooling. Although a physical model that describes a cooling process often involves
some sort of stochastic element, this will not be the case in Chapter 2, in which a
cooling idea for two trapped interacting atoms is explored. This process will be modeledusing deterministic (in the sense that quantum mechanics is deterministic for the wave
function) scattering models. Thus, this chapter is connected to the others only by the
subject of cooling. At the other end of the spectrum, Chapter 6 discusses the applicationof relaxation techniques to the field of nonlinear optical spectroscopy for understanding
non-Markovian scattering phenomena in recent experiments. This chapter utilizes thestochastic physics developed earlier in the thesis, but the subject has really nothing to dowith cooling. Nevertheless, the progression of chapters is intended to offer a somewhat
continuous train of thought.
Chapter 2 presents a novel cooling technique that makes use of Feshbach resonances in the scattering of two atoms. A few types of basic scattering models, using one
or more channels and different descriptions of the interaction, are first developed. Aftera brief example of using these models to describe molecular photodissociation, details of
the cooling scheme are presented. The scattering models are employed to simulate theworkings of the cooling scheme, and the results yield an understanding of how effective
the scheme is. The work discussed here was published in Ref. [1].
Chapter 3 develops the master equation for an atom with complicated internal
structure interacting with a three-dimensional polarization-gradient laser field. Detailsof the derivation of the relaxation operator for this system are relegated to Appendix B,
and a thorough derivation of the semiclassical optical Bloch equations for the system arepresented in Appendix C. The method of stochastic trajectories is presented, and thenapplied to the system at hand. This yields numerical results for the cooling dynamics

6of the fermionic alkaline-earth atoms 25 Mg and 87 Sr. These results are compared to thecooling behavior of atoms with simpler internal structure. The work presented here has
been published in Ref. [2].
The laser cooling of three-level systems using two lasers is presented in Chapter 4. The equations of motion are derived in one dimension for this system, and asparse-matrix method is used to obtain direct and exact solutions, bypassing the use
Monte Carlo techniques for this problem. A rich variety of phenomena is observed, theinterpretation of which leads to an understanding of the cooling technique in terms ofelectromagnetically induced transparency. The dynamical effects of scattering light off
atoms with asymmetric Fano lineshapes is discussed, and it is shown that such scattering can lead to lower temperatures than would be expected using a more simplistic
analysis with Lorentzian lineshapes. The work presented here has been submitted for
publication and a preprint version can be found in Ref. [3].
The master-equation models from the previous chapters are extended to include
an optical cavity in Chapter 5, with the goal of calculating cavity-cooling dynamics formolecules. The full master equation is then manipulated to yield a master equation
for just the molecule system, coupled to a Langevin equation for the cavity mode. The
cavity mode operator is then approximated as a classical field, which will be valid in thelimit of large photon numbers in the cavity. The set of coupled equations is then solved
to illustrate the effects of cooling on the population distribution, and the validity of thesemiclassical approximation is tested.
Lastly, Chapter 6 illustrates a possible idea for modeling non-Markovian dynamics
in the interparticle scattering of an atomic gas. A relaxation-type theory of atomic
collisions, originally applied to understanding static pressure broadening in gases, isextended for use the framework of nonlinear optical spectroscopy. The use of relaxation
theory should provide detailed and accurate scattering information and may lead toa new method for modeling the cutting-edge transient four-wave mixing experiments

7being performed today.

Chapter 2
Feshbach-Resonance Cooling

2.1

Introduction
The achievement in 1995 [19, 20] of quantum degeneracy in an atomic gas, known

as Bose-Einstein condensation (BEC) for bosonic constituents, was the result of yearsof development of various cooling techniques. In the end, this feat required not asingle cooling method, but a series of methods, each particularly suited to traversing
a certain temperature range of the gas. In the first condensate, a relatively hot gasof rubidium atoms was first cooled by basic Doppler laser cooling, after which moresophisticated polarization-gradient laser-cooling schemes produced further cooling, andfinally a creative series of evaporative-cooling trajectories did the trick.
Although many well-developed cooling schemes exist and have been successful fora variety of atoms, there is always a desire for new cooling ideas. An atom’s internal
structure might make laser cooling ineffective, or the details of the interatomic scattering
properties of an atom might not yield a favorable ratio of elastic to inelastic collisions,something that evaporative cooling requires. The scheme presented in this chapter is in
the spirit of the search for novel cooling techniques that may potentially find success insystems where other techniques have failed. It also tells a somewhat entertaining story
about those who attempt to defy the second law of thermodynamics.
For completeness, this chapter begins with a brief summary of the elements ofcollision theory necessary to understand the following development. The next section

9develops the quantum-defect theory of two particles in a harmonic trap, which will be
used to model the Feshbach-resonance cooling scheme. This section explores the singlechannel as well as the two-channel scattering models, utilizing a zero-range-potential
to describe the interparticle interactions. In order to provide more general insight,
the following section develops a two-channel square-well-potential version of the twoparticle theory. An example of using the previously developed theories for a perturbative
treatment of molecular photodissociation is briefly illustrated. Next, the Feshbachresonance cooling scheme is described, the results of calculations are presented, and theefficacy of the scheme is assessed.
Parts of this chapter are developed from research by this author which has been
published in Ref. [1], and some of the text and figures have been adapted from the work
therein.

2.2

Collision-theory basics
The purpose of this section is to provide a summary overview of the basic concepts

and results of two-particle collision theory. Once introduced, these concepts will be
useful in understanding the terminology used in developing the scattering models of
various sophistication later in the chapter. The general theory of potential scattering
theory is first presented, followed by a treatment of the so-called cold-collision regime,in which a great deal of simplifications can be made for the low-energy collisions which
occur in the quantum-gas experiments found in many of today’s laboratories.

2.2.1

Potential scattering theory

In understanding the physics of the collision of two particles, it is equivalent andsimpler to study the problem of a single particle colliding with a potential V (r). The

10time-dependent Schrödinger equation for this situation is
¸
·
∂
~2 2
−
∇ + V (r) Ψ(r, t) = i~ Ψ(r, t)
2m
∂t

(2.1)

where m is the mass of the particle, ~ is a constant related to Planck’s constant hhby ~ = 2π, and Ψ(r, t) is the wave function representing the particle. For V (r) realand independent of time, an assumption that is typically valid, stationary state wave
functions ψ(r) can be found which are solutions to the time-independent Schrödinger
equation,

·
−

¸
~2 2
∇ + V (r) ψ(r) = Eψ(r).
2m

(2.2)

Here, the particle’s energy E is defined as
E=~2 k 2
p2
=
,
2m
2m

(2.3)

where k is the wave vector, defined in terms of the particle’s momentum p as
k=

p
.
~

(2.4)

If V (r) goes to zero faster than 1/r for large r, then outside of a certain range
Eq. (2.2) describes a free particle. In this region, the wave function ψki for a given
initial wave number ki , can be written as the sum of an incoming plane wave and an
outgoing scattered spherical wave,·
iki ·rψki (r) ∼ A e
r→∞

¸
eikr
+ f (k, θ, φ)
.
r

(2.5)

where f (k, θ, φ), called the scattering amplitude and which is a function of the incoming
wave vector k and the spherical coordinates θ and φ, describes the amplitude of theoutgoing wave as a function of the direction and momentum of the incoming wave. It
can be shown that the differential cross-section is related to the scattering amplitude
by the expressiondσ= |f (k, θ, φ)|2 ,
dΩ

(2.6)

11
with the total cross-section then given by integrating the differential cross section overall solid angles,Z
σtot =

dσ
dΩ =
dΩ

Z
|f (k, θ, φ)|2 dΩ.

(2.7)

If V (r) is spherically symmetric,
V (r) → V (r),(2.8)

then ψki can be expanded in a series of products of Legendre polynomials and radialfunctions,ψki (k, r, θ) =∞X

Rl (k, r)Pl (cos θ),

(2.9)

l=0

where Rl (k, r) is a radial function, independent of angle, and Pl (cos θ) is the lth Legendre
polynomial. This series is referred to as a partial wave expansion, with each term of
the expansion representing an individual partial wave, parameterized by a single valueof angular momentum l. It will be useful to redefine the radial functions asul (k, r) ≡ rRl (k, r).

(2.10)

Then the separated Schrödinger equation for the ul radial functions, which is called theradial equation, is·¸
d2l(l + 1)
2−
− U (r) + k ul (k, r) = 0,dr2
r2

(2.11)

where
U (r) =

2µ
V (r)
~2

(2.12)

is the rescaled interaction potential. The asymptotic solutions to the radial equationare
Al (k)
sin [kr − lπ/2 + δl (k)],r→∞kul (k, r) ∼(2.13)

where Al and δl are constants (with respect to r), and δl is referred to as the phase shift
for the lth partial wave. In the partial-wave formalism, various scattering properties

12
can be expressed as sums of partial-wave contributions. For example, the scattering
amplitude is
∞
n
o
1 X
f (k, θ) =(2l + 1) e2iδl (k) − 1 Pl (cos θ),
2ik(2.14)

l=0

and the total cross-section is equal to the sum of all partial wave cross-sections,
σtotal (k) =∞Xσl (k),(2.15)

l=0

where the partial wave cross-sections areσl (k) =2.2.2

4π(2l + 1) sin2 δl (k).
k2

(2.16)

Cold collisions

In the radial equation for ul (k, r), Eq. (2.11), the l-dependent term can be added
to the bare potential U (r) and the two treated as an effective potential,
Ueff (r) = U (r) +l(l + 1).
r2

(2.17)The second term, called the centrifugal barrier, will cause the effective potential toincrease with increasing values of l. For small collision energies, this centrifugal barrier
will become insurmountable for partial waves beyond a certain value of l. For collisionsof low enough energy, only the l = 0 partial waves will contribute to the scattering; theothers will be prevented from colliding. Collisions occurring in this low energy range
are termed cold collisions. Since only s-wave (l = 0) collisions occur in this regime,
the number of partial waves that need to be considered is reduced, and many of thescattering properties simplify. This reduction of the partial-wave series given in Eq. (2.9)
to include only the first term is very useful as it can allow scattering properties to often
be determined analytically when working in the cold-collisions regime.
If a quantity called the scattering length is defined as the low energy limit of thetangent of the phase shift for an s-wave collision,tan δ0 (k),k→0ka = − lim

(2.18)

13then the scattering properties reduce to simple expressions involving a. In particular,
the scattering amplitude becomesf → −a,(2.19)dσ→ a2 ,dΩ k→0(2.20)

σtot → 4πa2 .

(2.21)

k→0

the differential cross section becomes

and the total cross section reduces to

k→0

A positive scattering length indicates that the interactions between the particles are repulsive, while a negative scattering length indicates that the interactions are attractive.

2.2.3

Feshbach resonances

Feshbach resonances have had an enormous impact on the study of cold quantum
gases in recent years, as have been extensively utilized experimentally to control theinteraction strength between atoms in these gases [?, 21, 22, 23, 24]. Controlled byadjusting the magnetic field applied to the gas, Feshbach resonances amount to a knob
by which experimentalists can tune the interparticle interactions. Such control, whenapplied to, for example, optical lattices containing one or more atoms at each site, canallow unprecedented insight into the physics of condensed matter systems.
A Feshbach resonance [25, 26, 27], as opposed to a shape resonance, is fundamentally a multichannel phenomenon, as it occurs for two atoms when their collision energy
becomes degenerate with a bound state in a closed collision channel, producing brieftransitions into and out of this state. It is primarily in order to describe such a resonancethat motivates the use of various multichannel models in the following sections.

14

2.3

Zero-range-potential model for two particles in a harmonic trap
To describe the dynamics of the Feshbach-resonance cooling scheme, it is neces-sary to utilize a quantum-mechanical model for two particles interacting in a harmonic
trap. Drawing from the basic collision theory presented in the preceding section, atheory with both a single channel as well as a two-channel theory will be presented.
The single-channel theory is simpler, and provides a good starting point, while the introduction of a second channel allows for a more accurate description of the scatteringby allowing inelastic phenomena to take place. Many papers have addressed the issueof the applicability of the single-channel model for various types of resonances, fromnarrow to broad. Such details are beyond the scope of the present discussion, and a
two-channel theory will be the primary tool utilized to model the Feshbach-resonancecooling scheme.
Although a variety of methods can be used to solve the zero-range-potential problem in a harmonic trap, this section will utilize a quantum-defect-theory approach.

2.3.1

One channel

In the past, a single-channel quantum-defect-theory model of two interacting
trapped particles has been successfully used to describe the coherent atom-molecule
quantum beats produced in an atomic gas with a Feshbach resonance, which has been
subjected to magnetic-field ramps [28]. The development of this model is outlined here
to provide a basis for understanding the subsequent development of a two-channel theory.
The interaction between the two harmonically trapped particles is described byan energy-independent regularized zero-range potential [28, 29, 30],
V (r) =∂
4π~2 a (3)δ (r) r,m
∂r

(2.22)

where r is the (vector) interparticle coordinate and a is the two-particle scattering

15
length. The Hamiltonian for the system is separable into center-of-mass and relative(interparticle) parts. This is due to the fact that the interaction potential given above
depends only on the relative coordinate, while the free-space kinetic energy operator fora pair of particles can always be written as the sum of a relative and a center-of-mass
term. The center-of-mass degree of freedom exhibits a trivial simple-harmonic-oscillator
spectrum, so only the relative motion is considered in the following. Also, since thezero-range interaction in Eq. (2.22 affects only s-wave interactions, the equation canbe restricted to consider only l = 0. Defining the reduced mass, µ = m/2, the timeindependent Schrödinger equation for the relative-coordinate wave function ψ(r) of two
atoms interacting in a harmonic-oscillator trap is given by
¶
µ 2 2
~ d
1 2 2
−
+ V (r) + µω r ψ(r) = Eψ(r).
2µ dr2
2
Defining the dimensionless units of length x = r/Losc , where Losc =(2.23)p~/(ωm/2) is the

reduced-mass oscillator length, and of energy ² = E/~ω, rescaling the s-wave eigenfunction asψ²,l=0 =u(x) 1√ ,x
4π

and noting that Eq. (2.22) is nonzero only at the origin, Eq. (2.23) becomes
¶
µ
1 d2
1 2
−
+ x u(x) = ²u(x),2 dx2 2(2.24)(2.25)

with the zero-range potential interaction imposing the boundary condition on u(x) at
the origin ofu0 (0)
Losc
=−
.
u(0)
a

(2.26)The solutions to Eq. (2.25), regular and irregular at the origin, can then be
found and cast in terms of the effective quantum number ν = ²/2 − 3/4. Applying thetechniques of quantum defect theory, the quantization condition for the scaled energy ²
in terms of the quantum defect β is obtained,3² = 2(n − β) + ,2(2.27)

16
where n = 0, 1, 2, 3, . . . . Enforcing the boundary condition on the logarithmic derivativeof u(x) at the origin, given in Eq. (2.26), yields the quantum-defect equation¡¢a 2Γ 2² + 43¡¢.
tan πβ = −Losc Γ 2² + 14

(2.28)

Equation (2.28) along with Eq. (2.27) can be solved self consistently to obtain theenergy eigenvalues ² for a specified scattering length a. The eigenenergies ² obtained
are illustrated in Fig. 2.1 as a function of the scattering length in oscillator units a/Losc .
Note that the energy levels undergo a shift as the inverse of the scattering length goes
through zero, or equivalently, the scattering length goes through infinity. This effect isdue to the introduction of a new bound state, and the characteristics of these level shifts
will provide the basic motivation for the Feshbach resonance cooling scheme presented
later in this chapter.
Note that a somewhat different approach to solving this problem is given in
Ref. [29], yielding the transcendental equation
¢
¡
+ 34
2Γ −²
L
2
¡ −² 1 ¢ = osc ,aΓ 2 +4(2.29)

which can be shown to be equivalent to the system of equations specified by Eq. (2.27)
and Eq. (2.28).

2.3.2

Two channels

In order to later model the effects of the magnetic field ramps for the Feshbachresonance cooling scheme, a two-channel Feshbach resonance model is now developed,
based on the single-channel model described in Ref. [28]. Both of these models describe
a two-atom Feshbach resonance for a harmonic trap, and utilize a zero-range potential
to describe the interaction between the two atoms. The two-channel model has theadvantage of allowing for a field-dependent resonance state. In the two-channel problem,the s-wave radial solutions for the relative coordinate r of the atom pair satisfy the

17

10
8
6
4

E 2
0
-2
-4
-8

-4

0

4

8

a
Figure 2.1: Eigenspectrum in oscillator units ~ω for two particles in a harmonic trap, asa function of the two-particle scattering length a in units of the reduced-mass oscillator
length Losc .

18matrix equationHψ(r) = Eψ(r),(2.30)

where the Hamiltonian matrix H is defined as
H=−

~2 d2
1
+ V (r) + E th ,2µ dr2(2.31)and ψ(r) is a two-component wave-function vector with components ψi (r) corresponding
to the ith channel. In Eq. (2.30), the first term on the right hand side is the two-channel
kinetic energy operator, where 1 is the identity matrix. The matrix V represents thetwo-particle interactions, with diagonal elements describing single-channel scatteringand the off-diagonal elements describing channel coupling. The matrix E th is a diagonal
matrix that describes the bare (zero interaction) splitting of the thresholds of the twochannels. As such, it will be a diagonal matrix, and in most cases the zero of energy
will be chosen so that only one element is nonzero. Here it will be written as0 0 
E th = ,
0 ∆²(2.32)

where ∆² is the energy shift of the second channel from the first channel.
Along the lines of the single-channel model developed above, the interactions
within each channel will be described by independent zero-range potentials, and a parametric channel coupling will be introduced. Using the same scaled length and energy
parameters from the single channel model, this results in decoupled time-independent
Schrödinger equations for the two channels,
¶
µ
1 2
1 d2
+ x u1 (x) = ²u1 (x)−2 dx2 2
µ
¶
1 d2
1 2
−
+ x u2 (x) = (² − ∆²) u2 (x),2 dx2 2(2.33)(2.34)

where ui for i = 1, 2 is the wave function for the ith channel, scaled as in Eq. (2.24).
The zero-range potential imposes a boundary condition at the origin, which is written

19here asu(x)d  1 dx
u2 (x)x=0

1
 − a1=ζζ  u1 (x)
u2 (x)− a12.

(2.35)

x=0

Equation (2.35) contains the in-channel scattering physics in the diagonal terms, withai being the single-channel scattering length for the ith channel, and with the channel
coupling represented in the off-diagonal elements by the parameter ζ.
Again applying a quantum-defect-theory treatment, along the lines of Ref. [28]for a single-channel model, a transcendental equation is obtained [31],
Ã ¡
!Ã ¡
!
¢
¢
3
2Γ − 2² + 34
2Γ − ²−∆²
+
1
1
2
4
¡
¢ −
¡
¢ −
− ζ 2 = 0,1
a1a
Γ − 2² + 14
Γ − ²−∆²+
2
2
4(2.36)

the solutions of which, for specified values of a1 , a2 , ζ, and ∆², are the eigenenergies of
the interacting system.
The aggregate scattering length for the full two-channel problem predicted bythis model, in the low-collision-energy limit when ω → 0, is
!−1
Ã
|ζ|2
1
+p
a(E, B) =.
a1
2µ∆²(B)/~2 − 2µE/~2 − 1/a2

(2.37)

This expression can be compared to the known scattering length, determined fromeither experiment or more detailed scattering calculations, to determine the values ofthe parameters a1 , a2 , and ζ. The parameters also affect the magnetic-field dependenceof the adiabatic energy states for the interacting system, and can also be fitted to.
Although it is often found that a given set of parameters will not provide accurate
scattering predictions over all regimes, they can be chosen to suitably describe thescattering over a certain range of interest.

2.4

Square-well-potential model for two particles in free space
In this section, the theory will be generalized to describe interactions between

the two particles described by square-well potentials, which, although not usually as20
analytically appealing as zero-range potentials, can provide a better description of shortrange interaction physics when such a thing is important. The theory presented here isbased on a model developed by Chris Greene [32].
The two-channel Hamiltonian will be written as Eq. (2.31) in the previous section,with the constant potential in the inner region given by
−V1 V12 
V (r) = 
 θ(r0 − r),
V12 −V2(2.38)and the upper threshold specified by0 0 
E th = 
.
th
0 E2

(2.39)The solution matrix, given by u(r), obeys u(0) = 0.
Define the matrix



µ
µ
¶³
¶ ²+V
´
−V
1
12
2m
2m 

W = − 2V (r) + E th − ²1 =

.
2
~
~
−V12 ² − E2th + V2

(2.40)The eigenvector and eigenvalue matrices corresponding to W are X and w2 , respectively:
W X = Xw2 . Two uncoupled single-channel equations for y(r) = X T u(r) are thenobtained. The regular solutions at the origin that satisfy these uncoupled equations are
simply0 
sin w1 r
y(r) = .
0sin w2 r

2.5(2.41)

Digression: application to molecular dissociation
As an example of the utility of the theory developed earlier in this chapter, thissection applies two-channel scattering models to the problem of molecular dissociationinduced by an oscillating magnetic field of perturbative amplitude. For the present
purpose, only processes in free-space will be considered, and the scattering equations
will be used neglecting the harmonic-oscillator trapping potential.

21
2.5.1

Zero-range potential

The model consists of two scattering channels, each with a zero-range potential at
the origin. The radial wave function satisfies the following free-space coupled equations,~2 d2
u1 = Eu1 ,2m dx2(2.42)~2 d2
u2 = (E − ∆²(B)) u2 ,2m dx2(2.43)−−

where E is the two-channel energy, ∆²(B) is the channel energy separation and is afunction of the magnetic field, and u1 and u2 are the lower and upper channel radial
wave functions, respectively. The effects of the true potentials in the channels are
enforced by a boundary condition at the origin that fixes the log derivative of the radial
wave function at that point. With the exception of the origin, the wave functions
everywhere satisfy the free-particle Schrödinger equations given above. The boundary
condition at the origin is given in Eq. (2.35).
Molecular dissociation consists of transitions from states below the lower-channel
threshold (both channels closed) to states in which one or both of the channels are open.
Thus, it is necessary to consider two separate cases: (a) the lower channel open and theupper channel closed, and (b) both channels closed. Dissociation processes consist oftransitions from the bound case (b) to the continuum case (a). The two cases will beconsidered individually.
(a) For ∆² > E > 0 (lower channel open, upper channel closed), the energynormalized solutions are





r
sin(k
x
+
δ)
u
(x)
1a
2m 
 1a 



,
=
2
πk1a ~
Da e−κ2a x
u2a (x)(2.44)

where Da is a constant. From the differential equations,Ea =2~2 κ22a~2 k1a= ∆² −.
2m
2m(2.45)

22
The boundary condition at the origin gives

 


ζ  sin δ k1a cos δ  −1/a1

=

,
∗
−κ2a Da
ζ
−1/a2
Da
or


k1a
−1/a1 − tan δ

ζ∗



ζ
 sin δ  = 0.
−1/a2 + κ2aDaThis equation has a non-trivial solution if
¯
¯
¯
¯
k
1a
¯−1/a1 −
¯
ζ
tan δ
¯
¯
¯
¯ = 0.
¯¯¯ζ∗−1/a2 + κ2a ¯(2.46)(2.47)(2.48)

Using the solution to this equation, and defining the energy-dependent scattering length
to be a(Ea ) = − tan δ(k1a )/k1a ,
1
1
|ζ|2
=
+q
.
a(Ea )
a12m(∆² − E ) − 1/a~2a(2.49)

2

Going back to Eq. (2.47) closed-channel amplitude is determined to beµ¶
k1asin δ 1+
Da = −
ζ
a1 tan δ
µ
¶
sin δ
1
1
=
−
.
ζ
a(Ea ) a1(b) For E − ∆² < 0 (both channels closed), the solutions are




r
−κ
x
1b
2m  eu1b (x)


=

.
πκ1b ~2−κx2bu2b (x)Db e(2.50)(2.51)

From the differential equations,
Eb = −~2 κ21b~2 κ22b= ∆² −.
2m
2m

The boundary condition at the origin gives

 
 
ζ  1 −κ1b cos δ  −1/a1

=
 ,
∗
κ2b Dbζ−1/a2
Db(2.52)(2.53)

23
or



 
−1/a
+
κ
ζ
1
1b

 1    = 0.
ζ∗
−1/a2 + κ2b
Db

This equation has a non-trivial solution if¯¯¯¯¯−1/a1 + κ1b
¯
ζ
¯
¯
¯
¯ = 0.
¯¯∗¯ζ
−1/a2 + κ2b ¯(2.54)(2.55)

Using Eq. (2.54), the upper-channel amplitude is found to be
1
Db =
ζ

¶
µ
1
,
κ1b −
a1

(2.56)Db =ζ∗.
κ2b − 1/a2

(2.57)

or equivalently,

2.5.1.1

Dissociation rate

At this point, time-dependent perturbation theory will be used to determine therate of transitions from case (b), where a true bound state exists, to case (a), where thelower channel is open, and the states of the coupled system are all continuum.
The perturbation being applied is a magnetic field, oscillating with frequency
νpert , with amplitude Bpert , and with an envelope of length tpert , as used, for example,
in experiments at JILA [33]. This type of oscillating field pulse is shown in Fig. 2.2.
This oscillation can be represented by¶µtsin (2πνpert t) ,B(t) = Bini + Bpert sin π
tpert(2.58)

where Bini is the magnetic field at which the pulse is centered, and Bpert is the amplitudeof the oscillations.
The first-order transition amplitude in time-dependent perturbation theory is−ic(1)n (t) =
~

Z t
t0

0

dt0 eiωni t Vni (t0 ),(2.59)

24

BBini
1
0.5
0.2

0.4

0.6

0.8

1

ttpert-0.5-1
Figure 2.2: An example of the type of magnetic field pulse used in molecular-dissociation
experiments. This pulsed magnetic field is used here as a perturbation in a in a twochannel model of dissociation. The amplitude of the pulse is Bini , the width of theenvelope is tpert , and the frequency of oscillation is νpert . In this plot, νpert = 10/tpert .

25
where Vni is the transition matrix element between the state |ii and the state |ni. The
perturbation that considered here is of the form


0
0

V (t) = ,
0 ∆²(B(t)) − ∆²(Bini )


0
0

=
,
0 µ (B(t) − Bini )(2.60)

where µ is the slope of ∆² with respect to B.
The transition matrix element for each of the two cases will now be considered.

2.5.1.2

Bound-bound transitions

Transitions that occur between case (b) and case (b) are bound-bound transitions.
The transition matrix element is


¶
µ
¶
µ
t
u1b0 (x)sin (2πνpert t) u2b0 (x).
 = u2b (x)Bpert sin π
u1b (x) u2b (x) V (t) tpert
u2b0 (x)(2.61)The transition amplitude is then given by
Z ∞2mµBpert(1)
cEE 0 (t) = √Db Db0dxe−(κ2b +κ2b0 )x
π κ1b κ1b0 ~2
0
¶
µ
Z t
t
0 i(E−E 0 )t/~
×
dt e
sin π
sin (2πνpert t)
tpert
t0

(2.62)

2mµBpert
1ζ ∗2= √f (E − E 0 , νpert , t),2
π κ1b κ1b0 ~ (κ2b − 1/a2 ) (κ2b0 − 1/a2 ) κ2b + κ2b0
where
f (E − E 0 , νpert , t) =¶µt0sin (2πνpert t) .
dt0 ei(E−E )t/~ sin π
tpert
t0

Z t

(2.63)

26
The bound-bound transition rate is then
¯
¯2
¯ (1)
¯
REE 0 (t) = ¯cEE 0 (t)¯
µ
¶
µ
¶2 µ
¶2
2mµBpert 2
1
1
|ζ|2
=
π~2
κ1b κ1b0 κ2b + κ2b0
(κ2b − 1/a2 ) (κ2b0 − 1/a2 )× f 2 (E − E 0 , νpert , t)¶µ¶2µµBpert 2 1
1
√
√
√
=
π
∆² − E + ∆² − E 0
EE 0×  ³q|ζ|22m(∆² − E) − 1/a2
~2

´ ³q

2m
(∆² − E 0 ) − 1/a2
~2

2

´

× f 2 (E − E 0 , νpert , t).
(2.64)
2.5.1.3

Bound-free transitions

Transitions that occur between case (b) and case (a) are bound-free transitions.
The transition matrix element is


µ
¶
¶
µ
t
u1b (x)
sin (2πνpert t) u2b (x).
 = u2a (x)Bpert sin π
u1a (x) u2a (x) V (t) tpert
u2b (x)(2.65)The transition amplitude is then given by
Z ∞2mµBpert(1)dxe−(κ2a +κ2b )x
cba (t) = √Da Db
π k1a κ1b ~2
0
µ
¶
Z t
t
0 i(Ea −Eb )t/~
×
dt e
sin π
sin (2πνpert t)
tpert
t0
µ
¶
2mµBpert ζ ∗ sin δ
1
1
1
1
= √
−
f (Eb − Ea , νpert , t).
2
ζa(Ea ) a1 κ2b − 1/a2 κ2a + κ2b
π k1a κ1b ~(2.66)

27
The transition rate is then calculated to be¯¯¯ (1) ¯2
Rba (t) = ¯cba (t)¯
¶
µ
µ
¶ µ
¶2 µ
¶2
2mµBpert 2 2
1
1 2
1
1
√
=
sin δ
−
a(Ea ) a1
κ2b − 1/a2
κ2a + κ2b
π k1a κ1b ~2
× f 2 (Eb − Ea , νpert , t)µ=2mµBpert
π~2¶2sin2 δ
k1a κ1bµ1κ2a + κ2b

¶2


q

|ζ|2
2m
(∆² − Ea ) − 1/a2
~2

2


µ

1
κ2b − 1/a2

× f 2 (Eb − Ea , νpert , t)µ¶µ¶2µBpert 2 sin2 δ
1
√
√
√
=
π
−Ea Eb∆² − Ea + ∆² − Eb

2 
2
|ζ|
|ζ|
 q

× q
2m
2m
(∆² − Ea ) − 1/a2(∆² − Eb ) − 1/a2
~2
~2
× f 2 (Eb − Ea , νpert , t).
(2.67)In the above equation, the expression for the energy-dependent scattering length given
in Eq. (2.49) has been inserted. From Eq. (2.48),
tan δ =

k1a
,
|ζ|2
−
1/a
1
κ2a −1/a2(2.68)

2x
so, using the fact that sin2 (arctan x) = 1+x2,

sin2 δ = ³

2
k1a
|ζ|2
κ2a −1/a2 − 1/a1

´2.

(2.69)

2+ k1a

For t > tpert , the function f 2 will be a sharply peaked function of νpert centeredon the transition energy.
At this point, it would also be possible to calculate the dissociation properties
taking into account interference between the two possible pathways of bound-bound-freeindirect and bound-free direct dissociation. This would yield behavior characterized by
Fano lineshapes [34], originally developed in the context of autoionization theory.

¶2

28
2.5.2

Square-well potential

Since both channels are closed, the solutions for r > r0 must decay exponentially,


−q1 r
N
e
1

~ phys (r) = ψ(2.70)
.
N2 e−q2 r
Smoothly matching the solution and its derivative at r0 requires~ phys (r0 ) = Xy(r0 )~z,
ψ

(2.71)~ 0phys (r0 ) = Xy 0 (r0 )~z,
ψ

(2.72)Xy(r0 )~z = D(r0 )~s,

(2.73)

Xy 0 (r0 )~z = D0 (r0 )~s,

(2.74)

or

with

and



−q
r
1
0 
e
D(r) = 
,
−q
r
2
0
e
 
N1 
~s =   .
N2(2.75)(2.76)

Eliminating z and defining the R-matrix gives¡¢
RD(r0 ) − D0 (r0 ) ~s = 0,(2.77)

which has a non-trivial solution when¡¢
det RD(r0 ) − D0 (r0 ) = 0,or¯¯¯(R11 − q1 ) e−q1 r0¯¯¯¯
R21 e−q1 r0

¯
¯
¯
R12
¯
¯ = 0.
¯−qr
20¯(R22 − q2 ) ee−q2 r0(2.78)(2.79)

29

Rate

EnergyFigure 2.3: An example of a typical molecular photodissociation spectrum obtainedusing the two-channel square-well model. Plotted is the photodissociation rate as afunction of energy, both in arbitrary units.

This gives the quantization conditione−(q1 +q2 )r0 [(R11 − q1 ) (R22 − q2 ) − R12 R21 ] = 0.

(2.80)The normalization condition requires
µ



¶



Z ∞ ¯
¯
 sin w1 r ¯ ~ phys ¯2
dr ¯ψ
(r)¯
1=
dr sin w1 r
X X
+
r0
0
B sin w2 r
Z r0
Z r0
Z ∞
Z ∞
=
dr sin2 w1 r + |B|2
dr sin2 w2 r + |N1 |2
dre−2q1 r + |N2 |2
dre−2q2 r .
Z r0

B † sin w2 r

0

0

T

r0

r0

(2.81)
Square-well calculations of the dissociation rate utilizing this model and to findbound-bound and bound-free transition rates using Mathematica yield spectra of theform shown in Fig. 2.3. The shape of this spectrum agrees with the predictions of moresophisticated scattering models, and comparison with such models or with experimentally obtained spectra can be used to determine the correct model parameters.

30

2.6

The cooling scheme
In this section, the workings of the Feshbach-resonance cooling scheme are pre-sented. This scheme makes use of the unique characteristics of a Feshbach resonanceto reduce the energy of pairs of externally confined atoms. It will be shown how thismethod can be used to cool pairs of atoms taken from a thermal distribution. The
basis of the cooling scheme originates from observing, as shown in Fig. 2.1 that thequantum-mechanical energy levels of two atoms in a harmonic trap shift by an energy
of approximately two trap quanta as the scattering length is adjusted through the resonance. If the scattering length is a function of some control parameter, then as thecontrol parameter is swept in one direction across the resonance, it will induce thisenergy shift. Throughout this section, the particular control parameter considered willbe the magnetic field B, which here is used to manipulate the atom-atom scattering
length a in the vicinity of a pole. However, it is important to note that, in other contexts, the shift of the energy levels could be induced by varying the detuning of an
off-resonant dressing laser, or by varying an electric field strength. The cooling schemethat is presented here in terms of the control parameter B can be extended to thoseother contexts in a straightforward manner. In the following, the basic mechanism of
the Feshbach resonance cooling process is first developed. After that, the method of
simulating the cooling, utilizing the two-body scattering models developed earlier in thischaper, will be presented. Next, the feasibility and effectiveness of the proposed coolingscheme are then illustrated through an application to a system of two atoms in a trap,
as would be encountered in a realistic experimental setup. Finally, the possibility ofapplying the scheme to atoms in optical traps or lattices will be discussed.
As mentioned earlier in this chapter, the Schrödinger equation for two identical,
interacting mass m atoms under spherical harmonic trapping potential of frequency ν
can be decoupled into two equations: one involving the three relative coordinates of

31
the pair of atoms, and another involving the three center-of-mass coordinates [29, 30].
In the Schrödinger equation for the relative coordinate of two trapped atoms, a central
potential is used. It is assumed for the time being that the center-of-mass coordinate
can be neglected, which will be the case if the system is translationally cold. An appliedexternal magnetic field B is introduced and its effect is described in a QDT context asa B-dependent quantum defect βEl (B). Then, the energy levels Enl (B) for the relative
motion of an atom pair are given by [30]
Enl (B) = (2n − 2βEl (B) + l + 3/2) ~ω,(2.82)

where ω = 2πν. Note that the quantum defect βEl (B) depends strongly on the relative
orbital angular momentum l of the pair but it varies only weakly as the radial oscillator
quantum number n is changed. The defect βEl (B) has only a weak dependence on theenergy when considered on the scale of an oscillator quantum. Explicitly, this can be
written as |dβEl (B)/dEnl | ¿ 1/~ω.
It will be shown later that the quantum defect for one relative partial wave l for an
atom pair (for example, the s wave, p wave, or d wave) must change by approximately
unity across the energy range kB ∆T of interest, and across the accessible range of thecontrol parameter, ∆B. The Feshbach resonance, which causes the scattering phaseshift
to change by a value of π is the result of this unit variation of βEl (B). A simple closedform expression can be derived for the quantum defect βEl , as is shown in Refs. [30]and [35]. This expression can be simplified when considering energies that are higher
than the energy of first few trap states to becomeµβEl (B) ≈ arctan¶
a (Enl , B) ~ω√,
2Losc Enl ewhere a(Enl , B) is the energy- and field-dependent scattering length and Losc =(2.83)p~/ (µω)is the characteristic oscillator length for the system, with the reduced mass µ = m/2
defined.

32
For simplicity, the present development is constrained to consider only an s-wave
resonance. However, the ideas described could be easily extended to higher partialwave resonances. The limiting low-energy scattering phaseshift is proportional to thewavenumber k = (2µE/~2 )1/2 for an s-wave Feshbach resonance. The E- and Bdependent scattering length, along the lines of Section 1.3, (dropping the subscript
l for notational efficiency, and because only the s wave is being considered) is thengiven byp
ΓE ~2 /(8µEn )a(En , B) = abg +,0 (B)En + (B − Bres )Eres(2.84)

where abg is the background scattering length. A zero-energy bound state is created
when the the magnetic field strength is equal to Bres , the resonance magnetic field . The
width ΓE in energy of the resonance is related to the width ∆ in the control parameter
0 (B)∆, where E 0
βEl by the relation ΓE = 2kabg Eres
res signifies the rate at which theresonance energy Eres changes as a function of the control parameter [36]. Figure 2.4
provides a plot of the characteristic s-wave energy levels En for the relative coordinateof two atoms in a spherical harmonic-oscillator trapping potential, as functions of theapplied magnetic field B for the particular case of a magnetic Feshbach resonance inthe scattering of 85 Rb(2, −2)+85 Rb(2, −2). A somewhat large and perhaps unrealistic trapping frequency of ν = 1 MHz is used in this in order to accentuate the fielddependence of the energy levels for visualization purposes. The idea for using such a
setup for cooling involves ramping the magnetic field B slowly from B1 to B2 and thenquickly back to B1 . A magnetic-field ramp for a more realistic situation encounteredin an experiment would likely cross more level shifts (that is, it would cover a larger
range of magnetic field). The quantum-mechanical state that undergoes an energy-level
shift for for the field value of B = B2 (which will be defined by n = Q later in thischapter) is indicated by a dashed line. In the example illustrated in Fig. 2.4, the model
0 = −3.5 MHz/G, Γ = 10 G,
resonance parameters considered are Bres = 155.2 G, EresB

33
and abg = −380 a0 , where a0 is the Bohr radius.
The particular choice of parameters affects the magnetic-field dependence of theadiabatic energy states. These parameters can thus be adjusted to provide satisfactory
agreement with experimental data in the regions presently of interest. For example, for
85 Rb, empirically determined parameters are found to be a = −435 a , a = 1.49 a ,10
20

and ζ = 0.00116 a−1
0 . Simulations of the field ramps can then be performed by specifying
an initial state of the system and numerically solving the Schrödinger equation.
The basic idea of cooling with a Feshbach resonance involves ramping the magnetic field through the region in which the energy levels shift by a value of approximately
2~ω. Figure 2.4 denotes the internal energy eigenvalues for the pair of atoms as a function of magnetic field B. If the pair of atoms is initially in an eigenstate at the fieldstrength of B = B1 , then a sufficiently slow field ramp from B1 to B2 will decrease theinternal energy of the pair by approximately 2~ω if the energy level undergoes a shift
in that field range (see inset of Fig. 2.4). More precisely, a slow ramp here is meant
to mean an adiabatic ramp. A fast ramp from B2 back to B1 would seem, ideally, toproject the wave function of the atom pair onto an eigenstate |n(B2 )i with the sameenergy as |n(B1 )i. This type of ramp is more precisely defined as a nonadiabatic ramp.
The result of the projection is to induce a net decrease in energy of approximately 2~ω.
Once this set of field ramps is finished, additional ramps can then be performed. Ideal
cooling is described in a diagrammatic manner by arrows in the inset of the figure,which has the same axes as the main figure. This diagram illustrates a process by which
population is transferred from point a to point b during the adiabatic (slow) field ramp
and from point b to back to point c during the nonadiabatic (fast) field ramp.
Figure 2.5 illustrates the application of diabatic vs adiabatic changes in terms ofa quantum particle in a box. For a particle initially in a one-dimensional box, the wave
function of the particle is depicted on the left side of the top and bottom of the figure.
In the upper portion of the figure, the box is rapidly expanded to a larger size. In this

34

a

60

_
E/hω

c

b

40

20

0
140

B1

150

B2 160170

B (G)Figure 2.4: Eigenenergies of an interacting two-atom system in a harmonic trap near a
Feshbach resonance, as a function of the magnetic field. The shift in the energies as theresonance is traversed can be seen, as well as the energy dependence of the location of
the shift in magnetic field. The inset illustrates an ideal series of population transfers
induced by magnetic field ramps for the Feshbach-resonance cooling scheme. This figure
is adapted from Ref. [1].

35
case, assuming that the change in box size is fast enough to be in the diabatic limit,
the initial wave function will project onto the eigenstates of the new system, and thewave function will be unchanged. In the bottom portion of the figure, the box is slowly
expanded to a larger size. If this change is slow enough to be in the adiabatic limit, thewave function will slowly track the corresponding eigenfunction of the box, resulting in
the displayed wave function. This difference of projection for fast changes and adiabatic
following for slow ramps is an important concept for understanding Feshbach-resonancecooling.
The dynamics of the cooling scheme can be modeled by specifying an initial stateof the system and then integrating the time-dependent Schrödinger equation,i~∂|ψ(B, t)i = H(B) |ψ(B, t)i .
∂t(2.85)The system state vector can be expanded in a set of eigenkets of the the system Hamiltonian for a given magnetic field,|ψ(Bi , t)i =Xa(i)νi (t) |νi i ,(2.86)

νi(i)where the aνi (t) is the expansion coefficient of the |νi i energy eigenstate for the systemwith a magnetic field of Bi , which has eigenvalue ~ωνi ,H(Bi ) |νi i = ~ωνi |νi i .

(2.87)

Consider propagating the system from an initial magnetic field value of B1 at time t1
to a final value of B2 at time t2 . Using the expansion in Eq. (2.86) for the value of theinitial magnetic field, Eq. (2.85) implies thatd (1)a (t) = −iων1 a(1)
ν1 (t).
dt ν1(2.88)

It has already been illustrated in Fig. 2.4 how the eigenvalues ων1 change as a functionof magnetic field. This time dependence, ωνi = ωνi (t), by way of the time-varying

36

Figure 2.5: Upper figure: Wave function in a box before and after rapidly (diabatically) expanding the box. Lower plot: Wave function in a box before and after slowly(adiabatically) expanding the box.

37
magnetic field, yields a solution to Eq. (2.88) of
R t2−i t dt ων1 (t)
(1)
1
a(1)
.
ν1 (t2 ) = aν1 (t1 )e(2.89)

Putting things in terms of the eigenstates of the system for the new magnetic field B2 ,a(2)
ν2 (t2 ) =

X

R t2

−i t dt ων1 (t)1hνf |νi i a(1),
ν1 (t1 )e(2.90)

νiwhere the overlap matrix hν2 |ν1 i is defined by
aνj =

X

hνj |νi i aνi ,(2.91)

νiand can be determined in a straightforward manner from the QDT development of twochannel model presented in Section 2.3.2. If the system is initially diagonalized over arange of magnetic-field values of interest for the problem at hand, these eigenenergies
can be used to perform numerical integrations in Eq. (2.90) and, supplying initial valuesof the expansion coefficients and the appropriate overlap matrices, the time evolution
can be determined.

2.7

Results and experimental possibilities
Upon performing simulations of the cooling using the numerical propagation tech-nique presented in the preceding section, the effect of the adiabatic and nonadiabatic
field ramps can be determined. Initially it is assumed that the atom pair is in a purequantum state at a magnetic field of B = B1 . As predicted by the above qualitative
discussion, the adiabatic field ramp (B1 to B2 and shown in Fig. 2.4) decreases theenergy of the atom pair regardless of the initial eigenstate chosen. A nonadiabatic ramp
(B2 to B1 and also shown in Fig. 2.4) for a state which is initially at a field value of B2
and assuming that this initial state is not degenerate with the resonance state, causes
the system to project onto a state at B = B1 having roughly the same energy as theinitial state at B2 . If, however, the pair is initially in a state that is degenerate with the38resonance state for the field value of B = B2 , the fast ramp is found to result in a strong
projection of the initial state onto the resonance state. This is to be expected since thecharacter of the dressed state of the interacting system having the same energy as thebare resonance state is very similar to the character of this bare resonance state. The
projection between these two states is thus very effective. In this case, the atom pairmust gain energy due to the fact that the resonance state at B = B1 has a energy thatis larger than that of the initial state at B = B2 .
It is now possible to generalize the above argument to allow an initial quantum
state that is in a mixed state. From this, it can be shown how cooling can be performed
for a more realistic ensemble of atom pairs instead of a pure state. If an atom pair istaken from a thermal distribution of pairs, the probability of occupying a level withrelative-motion energy En is given in terms of a Boltzmann factor by e−En /τ /Z(τ ) withP
τ = kB T , where kB is Boltzmann’s constant, Z(τ ) = i e−Ei /τ is the partition function,noting that the sum runs over all states of the system, and T is the temperature of thesource distribution from which the two atoms were taken. Based on the conclusions forfield ramps applied to pure states above, the application of a cooling cycle, defined asa slow field ramp from B = B1 to B = B2 , followed by a fast ramp back to B = B1 ,will exhibit two different types of behavior when applied to a mixed state. The first
behavior is to decrease by approximately 2~ω the energy of the population in states
which happen to undergo a full shift in the energy level between the field strengths of
B = B1 and B = B2 . The second behavior is to increase the energy of the populationthat happens to be the state that is degenerate with the resonance state, in other words,
the state which is undergoing an energy shift at B = B2 . This happens because thepopulation is moved to the state or states that are degenerate with the resonance at
B = B1 . The state from which the population is heated originates will be denoted fromthis point by n = Q.
Figure 2.6 illustrates the results of a numerical simulation of a single cooling cycle

39applied to a mixed, thermal distribution of s-wave states initial state of the system. The
same model Feshbach-resonance parameters were used as the example shown in Fig. 2.4.
The black line in Fig. 2.6 represents an initial s-wave relative-coordinate probability
distribution for the the atom pair in a harmonic trap, again with frequency ν = 1 MHz,and for a temperature T = 1 mK of the source. The red line denotes the same s-wave
probability distribution after the application of a slow and a fast magnetic field ramp,
for a value of the cooling parameter Q = 15. It is clear that the application of a single
cooling cycle transfers the population of the n = Q state to states with much larger
energy. In the example in the figure, this state is n ≈ 88, as is evidenced by the spikeseen in Fig. 2.6. As can also been seen, the field ramps transfer the population of eachstate with n > Q, states higher in energy than the Q state, each to the state next-lowest
in energy, which has approximately 2~ω less energy than the initial state. The numericalsimulations reveal that the net result of these two types of processes is to increase theoverall average energy of the system. That is, the energy increase from the few cases
where the pair gains a large amount of energy overwhelms the total energy decrease
due to the many cases where the pair loses a small amount of energy. It turns out
though, that this behavior is not entirely unexpected, as can be seen in Ref. [37] which
discusses the general nature of cooling processes based on time-dependent Hamiltonians.
However, it is important to note that, since the atoms gaining energy can be displaced
to an arbitrarily high energy state, as determined by the values of the extremum of the
field ramp (B = B1 in Fig. 2.4), it should be possible to remove these atoms by utilizing
methods similar to those used for atom removal in evaporative cooling.
At this point, it is interesting to observe the comparison of the Feshbach-resonancecooling scheme with the well-known evaporative cooling technique by which the mostenergetic atoms in a gas are allowed to escape, resulting in a colder remaining gas.
The process described here is similar to a conceptual “fractional” evaporative coolingexperiment, whereby the parameter Q is equivalent to the evaporative-cooling cutoff

40

Occupation probability

0.1
0.04

0.08

0.03

0.06

0.02

0.04

0.01
0.02
0

10
0

Q

Q

20

25
50
75
Pair s-wave state n

25
100

Figure 2.6: Illustration of the effect of a single Feshbach resonance cooling cycle for an
atom pair taken from a thermal distribution with T = 1 mK and trapped with frequency
ν = 1 MHz. The black (red) line represents the population distribution before (after)
application of one slow and one fast magnetic field ramp. The state Q (here Q = 15) isindicated. This figure is adapted from Ref. [1].

41
parameter and, instead of all the population with n > Q being removed, only a fractionof this population is skimmed off and removed from the trap. The similarity between
the two methods can be clearly seen in Fig. 2.6. The nature of the cooling scheme
are more explicitly shown in Fig. 2.7. Note, however, that cooling with a Feshbachresonance scheme is differentiated from evaporative cooling by the fact that it is notfundamentally of a statistical nature. The more that is known about the atom pair’s
initial state, the more effectively the pair can be cooled. In fact, if the initial state of
the pair is known with exact precision, not pulled from a thermal reservoir (as is thecase for a pure state), it is possible to design a sequence of field-ramp that will allow thetransition of an atom pair to the ground state of the trap 100% of the time in principle.
As with evaporative cooling, the requirement of atom removal leads to considerations of efficiency. An estimate of the efficiency can be made by assuming that the
population of the state Q is removed from the trap, while the population of all stateswith n > Q are moved to the next-lowest state, that is, to states with n − 1. This
assumes that the range of the field ramps is such that the heated fraction (n = Q) ends
up at an energy corresponding to negligible thermal population (as in Fig. 2.6), andthat all population above a specified energy can be removed. If the level energies at
B = B1 are approximated by En (B1 ) ≈ 2n~ω, the probability to remove an atom pairduring a cycle isPrem (Q, τ ) =e−2~ωQ/τ
.
Z(τ )(2.92)

The average energy decrease in a cooling cycle is due to the energy of the n ≈ Qpopulation removed from the trap, plus the energy loss for states n > Q:
∆E(Q, τ ) = (2~ωQ + hECM i)∞X
e−2~ωQ/τ
e−2~ωn/τ
+
2~ω
.
Z(τ )
Z(τ )

(2.93)

n=Q+1

Noting that hECM i = 3τ (since hEtot i = 3τ for a single atom in a harmonic trap), andP−2~ωn/τ ≈ e−2~ωQ/τ τ /2~ω, Eq. (2.93) becomeswith ∞
n=Q+1 e∆E(Q, τ ) =e−2~ωQ/τ
(2~ωQ + 4τ ) .
Z(τ )(2.94)

42

Figure 2.7: Diagram illustrating the effect of the Feshbach-resonance cooling schemeon an s-wave distribution of atoms. Shown is the population in arbitrary units asa function of the relative-coordinate quantum number n. Population from the staten = Q is increased in energy, and assumed to be removed from the trap. All population
with n > Q is decreased in energy by one relative-coordinate oscillator unit of energy.

43
The energy efficiency Eeff , defined as the amount of energy removed per atom removed,is then given byEeff (Q, τ ) = 2~ωQ + 4τ.

(2.95)

Since Q determines the efficiency of the cooling process, it is referred to as the coolingparameter. Results from the numerical model indicate that Eq. (2.95) provides a good
estimate of the efficiency.
The time scale for one cooling cycle is determined by the speed of the adiabatic
field ramp. This speed in turn is determined by the strength of the coupling between theresonance state and the trap states. The smaller the coupling for an avoided crossing,
the slower is the field ramp required to maintain adiabaticity. The coupling between
the resonance state and the trap states is related to the resonance width parameter ΓE ,which can be used in a Landau-Zener estimate of the transition probability [36],µ∼Ptr = exp −2
ωΓE|dB/dt| |dE/dB|¶.

(2.96)

Motivated by the possibility of experimentally trapping a small, deterministic
number of atoms [38, 39], the experimental feasibility of this cooling scheme is now explored. A Feshbach-resonance cooling experiment involves a sequence of cooling cycles.
As discussed above, a single experiment could result in a heated atom pair, which inturn would be lost from the trap. To see the effect of multiple field ramps, Eqs. (2.92)
and (2.93) can be iterated. For a variety of cooling efficiency parameters Q, Fig. 2.8
shows the probability for an atom pair to remain trapped vs. the average total kinetic
energy (the energy of both the relative and the center-of-mass degrees of freedom) of
the two atoms in oscillator units. Included in this calculation is the probability that theatom pair is in an s-wave state to begin with, because the field ramp has no effect onother partial waves. It is assumed that rethermalization occurs between cooling cycles,
which could be ensured by, for example, introducing a slight anharmonicity into thetrapping potential. (In the absence of anharmonicity, the relative and center-of-mass44
degrees of freedom would remain uncoupled.)
To be more specific, a crossed-beam optical dipole trap [40] can be considered,
which offers a good blend of large trap frequency (for a large s-wave fraction), isotropy,
and anharmonicity (for rethermalization between the relative and center-of-mass degrees
of freedom). Assuming the dipole trap has an average frequency of ν = 10 kHz andcontains two atoms taken from a source with temperature T = 8 µK (hEtot i/~ω = 100),
it is seen from Fig. 2.8, solid line, that a temperature of 0.16 µK (hEtot i/~ω = 2, both
atoms in the ground state) could be reached 10% of the time by performing less than2000 cooling cycles. For a range in magnetic field for the ramps of ∆B ≈ 1 G, and using
Eq. (2.96) with Ptr = 0.1, it can be seen that such a series of field ramps could take placein under 1 s. A pertubative calculation accounting for the trap anharmonicity indicatesthat rethermalization between the relative and center-of-mass degrees of freedom should
occur on a time scale comparable to a single ramp time for a crossed-beam dipole
trap. This will ensure that the relative s-wave distribution will rethermalize with eachramp and that the cooling of the relative coordinate will also cool the center-of-mass
coordinate (both of which have been assumed up to this point).
Feshbach resonance cooling could also be applied to atom pairs in an opticallattice. In this case, field ramps could be performed on the lattice ensemble of atom
pairs, with a certain percentage of sites resulting in cooled pairs, while other sites will
have either zero atoms or one (uncooled) atom. It may also be possible to prepare
the optical lattice by some other means to have a high probability of exactly double
occupancy at each lattice cite (see for example Ref. [41]). From such an initial state,a Feshbach resonance cooling scheme could be used to efficiently cool atom pairs to
low-lying trap states.

45

Probability to remain trapped

1
0.8
0.6

1
0.8
0.6
0.4 0.4
0.2
0
0.2
2
3
4
5
6
10 10 10 10 10
Number of cycles
0
100
80
60
_40
<Etot> / h ω

200Figure 2.8: Probability that a pair of atoms remains trapped vs. the average total
kinetic energy of the two atoms in oscillator units (note that kB T = hEtot i/6 for two
harmonically trapped atoms). Three different cooling parameters are used: 2~ωQ = 5τ(solid line), 9τ (dashed line), and 12τ (dot-dashed line). It is assumed that rethermalization occurs between cooling cycles (see text), although this scheme does not necessarily
require it. Inset: probability to remain trapped vs. the number of cooling cycles for thesame three cooling parameters. This figure is adapted from Ref. [1].

46

2.8

Conclusion
In summary, this chapter has developed a two-body theory that describes how

the energy of an atom pair can be reduced. The resulting novel cooling scheme, whichmakes use of Feshbach resonances, offers a viable new means to manipulating small,deterministic numbers of trapped atoms [38, 39] with present-day technology. Since it
is not clear at present how efficiently other cooling methods such as evaporation can beapplied to small atom samples, this proposal may prove quite useful. Extension of thisscheme to atom clouds is possible.

Chapter 3
Dynamics of Multilevel Atoms in Three-Dimensional
Polarization-Gradient Fields

3.1

Introduction
This chapter continues the discussion of theoretical descriptions of cooling byfocusing on the technique of using lasers to cool atoms. In particular the dynamics of
atom-laser interactions is analyzed for atoms having multiple, closely spaced, excitedstate hyperfine manifolds. The atom’s motion is described in a polarization gradient
field created by a three-dimensional laser configuration. .
The system is treated first in a semiclassical approximation, in which the centerof-mass of the atom is approximated to be a classical variable rather than a quantummechanical operator. This semiclassical problem allows much easier solutions than afully quantum problem, and yields quantities with intuitive appeal because they areclassical in nature. For example, the friction that an atom feels while moving in the
optical field can be obtained, as well as the diffusion coefficient for the atom’s motion.
These two quantities work against each other to obtain a steady-state thermal distribution of a given temperature. In the end, such semiclassical approximations can be quite
problematic and limited, as is the case when the atomic momentum approaches therecoil limit. Other problems can occur also: it is not unusual for semiclassical methods,
when applied naively, to yield unphysical negative temperatures.
After this, the system is treated fully quantum mechanically, including the atom’s

48
center-of-mass degree of freedom. The master equation describing this system is developed, and then specialized to the low-intensity limit by adiabatically eliminating theexcited states. It is shown how this master equation can be simulated using the Monte
Carlo wave function technique, and details are provided on the implementation of thisprocedure.
Monte Carlo calculations of steady state atomic momentum distributions for two
fermionic alkaline earth isotopes, 25 Mg and 87 Sr, interacting with a three-dimensionallin-⊥-lin laser configuration are presented, providing estimates of experimentally achievable laser-cooling temperatures.
The complex behavior that occurs when a multilevel atom interacts with polarizationgradient fields has been of interest for some time now. Sub-Doppler cooling [42] occurs
because of elaborate optical-pumping processes produced by laser light in atoms with
sublevel structure, as seen, for example, in the lin-⊥-lin and the σ+ -σ− laser configurations. The semiclassical understanding of these interactions [16, 17, 43, 44, 45, 46] in one
or more dimensions has led to a reasonably good qualitative understanding of the underlying mechanisms. Semiclassical analysis has even in some cases provided quantitative
predictions of sub-Doppler laser cooling temperatures measured in experiments [46].
However, the most direct route to a quantitative understanding of atom-laserinteractions is via a fully quantized master equation for the atom, in which the centerof-mass motion of the atom is taken into account quantum mechanically. This allows
behavior at low laser intensities and low atomic velocities, the regime laser cooling
strives to reach, to be described correctly. The drawback of solving such a masterequation, however, is the large number of basis states required for the calculation, due
to the additional momentum states. This problem becomes especially pronounced whenattempting to model three-dimensional systems, where the state space grows as the cubeof the number of one-dimensional momentum states needed.
The Monte Carlo wave-function (MCWF) technique, introduced in the early49
1990’s has allowed significant progress to be made on the subject of atom-photon interactions in three dimensions as well as lower-dimensional calculations. The MCWF
technique is a simulation procedure for the master equation that involves propagation of
single stochastic wave functions, rather than density operators, with random processes
occurring at random intervals due to interactions with the photon field that cause spontaneous emission. It has been shown that this method is equivalent to the masterequation in the limit of a large number of independent stochastic wave functions [10].
The MCWF technique has been successfully utilized to calculate three-dimensional subDoppler laser cooling temperatures for atoms with Zeeman degeneracy in the groundand excited states [11].
The majority of the research done on laser cooling has involved essentially twolevel systems, consisting of a ground state and an excited state, which may or may
not contain degenerate sublevels. However, some investigations have explored atomic
systems in which multiple distinct excited states come into play. In particular, the useof bichromatic laser fields [47, 48] to cool three-level Λ systems have been extensively
studied (see Refs. [49, 50, 51] for example).
The focus in this chapter is primarily on monochromatic laser cooling for atoms
with multiple closely spaced hyperfine excited-state manifolds. Figure 3.1 provides agraphical illustration of this type of atomic configuration. This situation is of importance, for example, in alkaline-earth atoms with nonzero nuclear magnetic moment. If
the excited state manifolds are spaced in energy on the order or smaller than the excitedstate linewidth Γ, coherences between these manifolds become nonnegligible, and canhave a significant effect on the optical pumping processes required for sub-Doppler cooling and on the dynamics of the atom-photon interaction. Sub-Doppler laser cooling was
experimentally identified in fermionic 87 Sr [52], despite significant spectral overlap in the
excited state. At the time, it was hypothesized that the large ground-state degeneracy
in 87 Sr (due to the large nuclear spin I = 9/2) was somehow able to overcome the de-50crease in cooling due to the spectral overlap. Other systems with spectral overlap in the
excited state are 39 K [53], 7 Li [54], and the fermionic isotopes of Yb [55]. In 87 Rb, theeffects of excited-state spectral overlap on the effectiveness of velocity-selective coherent population trapping have been explored, both experimentally and theoretically [56].
The goal in this chapter is to provide a detailed discussion of the theoretical techniques
required to model such systems realistically.
The structure of this chapter is as follows. Section 3.2 develops the master equation for a laser-driven atom with multiple excited-state manifolds. Within this section,
the a general fully quantized master equation is first derived, pulling from the detailed
derivation of the spontaneous-emission relaxation operator in Appendix B. Next, thismaster equation is specialized to the low-intensity limit, which reduces the number ofinternal atomic states necessary to describe the problem, and which is justified for the
low-power laser configurations which result in optimum cooling of the gas. Section 3.3
then goes on to demonstrate the semiclassical method of solving this problem, withresults provided which will aid in a qualitative understanding of the later, more quantitative results obtained from the Monte Carlo. Section 3.4 introduces the MCWF
technique and applies it to the low-intensity master equation. Section 3.5, full Monte
Carlo master-equation simulations are performed for 25 Mg and 87 Sr atoms in a threedimensional lin-⊥-lin laser configuration as an example of using this technique determine
expected temperatures for these atoms in a laser cooling experiment.
Parts of this chapter are developed from research by this author which has been
published in Ref. [2], and some of the text and figures have been adapted from the work
therein.

3.2

Master equation
In this section we develop the master equation describing a multilevel atom in-teracting with a coherent laser field and coupled to a vacuum photon field. It is this

51

d3

m=-2 -1

0 +1 +2

d2
d1

m=-1

0 +1

m=0

Fe3=Fg+1

Fe2=Fg

Fe1=Fg-1

wL

m=-1

0 +1

Fg

Figure 3.1: Energy level diagram of an atom with multiple hyperfine manifolds. If theenergy spacing of the excited-state manifolds are of the order or smaller than the naturallinewidth of the transition, the usual sub-Doppler cooling transition (Fg ↔ Fe = Fg + 1)
is not isolated and the other manifolds must be taken into account. This figure is adaptedfrom Ref. [2].

52equation, with quantized atomic center-of-mass, that will provide an accurate description of atom-photon dynamics, and this master equation will provide the basis for theMonte Carlo simulations that will be discussed later.

3.2.1

General form of the master equation

The full Hamiltonian for the atom-laser system plus the radiation field is
H = HA + HR + VA−L + VA−R ,where HA =

(3.1)

P

P2
i ~ωi Pi + 2m is the bare atomic Hamiltonian, HR is the vacuum radiation

field Hamiltonian, and VA−L and VA−R are the atom-laser and atom-radiation field
coupling terms, respectively. In the atomic Hamiltonian, Pi is a projection operator
onto the i-th internal excited-state manifold, ~ωi is the energy of the i-th excitedstate manifold relative to the ground-state manifold, P is the atomic center-of-mass
momentum operator, m is the atomic mass, and the sum runs over all excited-state
manifolds. Equation (3.1) assumes that the effects of atom-laser and atom-radiationfield coupling are independent [8].
Equation (3.1) can be viewed in terms of system-reservoir interactions. The system consists of the atom, the laser, and their interaction. The system Hamiltonianis
HS = HA + VA−L .

(3.2)

The reservoir is the vacuum radiation field, having many more modes than the system.
With the Markov approximation, along with a few other approximations, the masterequation is then given byσ̇ =i[σ, HS ] + Lsp [σ].
~

(3.3)

The operator σ is the system reduced density operator element, i.e., the reservoir degrees
of freedom have been traced over, σ = TrR ρ. The remaining term, Lsp [σ], encompasses

53
the interaction between the atom and the vacuum photon field, and provides for thephenomenon of spontaneous emission.
The relaxation operator due to spontaneous emission, which is derived in detailin the Appendix B, is given by3Γ
Lsp [σ] =
8π

Z
d2 Ω

XX

†

e−ik·R (²∗ · A(i) )σ(² · A(j) )eik·R²⊥k i,j

−

3Γ
16π

Z

d2 Ω

XXh
†
(² · A(i) )eik·R e−ik·R (²∗ · A(j) )σ²⊥k i,ji†+ σ(² · A(i) )eik·R e−ik·R (²∗ · A(j) ) , (3.4)†

where A(i) and A(i) are vector raising and lowering operators, respectively, between
the ground state and the ith excited state, R is the atomic center-of-mass position, kis the direction of the photon emitted in the relaxation process, and Γ is the decay rateof the exited states. The integral is performed over solid angle in the vector k and thesum over ² ⊥ k refers to the two polarization directions perpendicular to k. Note that
here and throughout this chapter, it is assumed that each of the excited-state hyperfine
manifolds has the same lifetime τ = Γ−1 . Expanding these vector operators in a basis√of spherical unit vectors, ²̂±1 = ∓(x̂ ± iŷ)/ 2 and ²̂0 = ẑ, we haveA(i) =

X(−1)q ²̂−q A(i)
q .

(3.5)

q=0,±1The spherical components of the vector operators are
Aq(i) =X

αFg ,Fei ,Mg ,Mei ,Jg ,Je ,I |Jg IFg Mg i hJe IFei Mei |(3.6)

αFg ,Fei ,Mg ,Mei ,Jg ,Je ,I |Je IFei Mei i hJg IFg Mg | ,(3.7)

Mg ,Mei†Aq(i) =X
Mg ,Meiwhereq(2Fg + 1)(2Fei + 1)(2Je + 1)






1
Fei  Jg Fg I  Fg. (3.8)
×





−Mg Mg − Mei Mei
Fei Je 1

Fg +Fei +Mg +Je +I

αFg ,Fei ,Mg ,Mei ,Jg ,Je ,I = (−1)

54
Eq. (3.4) is written in a way that makes explicit that it is in Lindblad form [57, 58,59]. As will be shown later, it is important for the relaxation operator to be of this formin order to make use of the MCWF technique. Because the complex exponentials in thesecond line cancel each other, the remaining integral over solid angle can be evaluated,
whereby Eq. (3.4) can be equivalently written as [10]3ΓLsp [σ] =
8π

Z
d2 Ω

XX

†

e−ik·R (²∗ · A(i) )σ(² · A(j) )eik·R −²⊥k i,j

ΓX[Pei σ + σPei ] . (3.9)
2i

The atom-laser interaction term, which is given in the electric-dipole approximation by
VA−L (R, t) = −D · EL (R, t),(3.10)

where EL (R, t) is the electric field of the laser and D is the electric dipole operator. As
is typical, the laser is treated as a classical field, since it is a densely populated modeof the electric field. The laser electric field can be written in terms of its positive and(+)

negative frequency components, EL (R, t) = EL (R)e−iωt + c.c., and then expanded
into spherical components,(+)

EL (R) =E0 X(−1)q aq (R)²̂−q ,2(3.11)

q=0,±1

where E0 is the electric-field amplitude and aq (R) are the expansion coefficients. Making
the rotating-wave approximation, so that
(+)

(−)

VA−L (R, t) = −D(+) · EL (R)e−iωt − D(−) · EL (R)eiωt ,where D(+) =

P

(−) =
i Pei DPg and D

VA−L = −(3.12)

Pi Pg DPei , we findΩX
Di (R)e−iωt + H.c..
2(3.13)i

The previous equation introduced the atom-laser raising operator,Di† (R) =X
q=0,±1†aq (R)A(i)q ,(3.14)

55
and lowering operator,X

Di (R) =a∗q (R)A(i)q ,(3.15)

q=0,±1and also introduced the so-called “invariant” Rabi frequency,Ω=E0 hJe ||D||Jg i√,
2Je + 1

(3.16)

where hJe ||D||Jg i is the reduced dipole matrix element between the ground and excitedstates. This form of a Rabi frequency, defined in terms of the reduced matrix element
between the J = Jg ground state and the J = Je excited state, is convenient because,in general, Rabi frequencies for transitions to different excited-state manifolds will not
be the same.
Next, observe that the second term in Eq. (3.9) is comprised of excited-state
projection operators both pre- and post-multiplying the system density operator. Thus,
it is clear that this term can be absorbed into the free-evolution commutator term in
Eq. (3.3), allowing the master equation to be equivalently described by Hamiltonian
evolution determined by an effective Hamiltonian Heff , plus a term which is commonly
called a jump term, and which cannot be written in the form of a commutator with thesystem density operator. This results in
´ 3Γ
i³
†
Heff σ − σHeff+
σ̇ = −
~
8π

Z
d2 Ω

XX

†

e−ik·R (²∗ · A(i) )σ(² · A(j) )eik·R , (3.17)²⊥k i,jwhere the effective Hamiltonian Heff is given by
µ
¶
P2 X
Γ
−
~ δi + i
Pei + VA−L ,Heff =
2m
2(3.18)iwhere VA−L is as given in Eq. (3.13). In obtaining Eqs. (3.17) and (3.18), the usual
rotating-frame transformation, which removes the free-evolution atomic Bohr frequencies from the problem, has been made. The more relevant frequencies are instead thelaser detunings δi = ω − ωi from the ith excited-state hyperfine manifold. The masterequation given in Eq. (3.17) is fully general, but has been written in a form that will

56
facilitate setting up a stochastic wave function simulation using the MCWF techniquedescribed later.

3.2.2

Master equation in the low-intensity limit: adiabatic eliminationof excited states

The master equation just discussed will now be specialized to the limit of low
laser intensity. Specifically, this limit is valid when the saturation parameter for theatom in the ith excited-state hyperfine manifold,
si =

Ω2 /2
,
δi2 + (Γ/2)2(3.19)

is small, which occurs when the laser intensity is small or the laser detuning from theatomic transition is large. In this limit, the excited states are said to adiabaticallyfollow the ground states. The excited states can then be eliminated from the equationsof motion, resulting in a master equation in terms of only the ground-state sub-densitymatrix,
σgg = Pg σPg .

(3.20)In this limit, the master equation becomes (see section 8.3.3 of Ref. [60])σ̇gg = −´i³heff σgg − σgg h†eff +
~

Z
d2 Ω

XX
†
(²∗ ·B(i) (R, k))σgg (²·B(i) (R, k)). (3.21)²⊥k i,jThe new effective Hamiltonian is given byµ¶
P 2 X siΓ†+~ δi − iD(i) (R)D(i) (R).
heff =
2m
2
2(3.22)i

The new decay raising and lowering operators are given byr†
Bq(i) (R, k) =andr
Bq(i) (R, k) =3si Γ (i) † ik·R (i)A e
D (R),
8π q(3.23)

3si Γ (i) −ik·R (i) †A e
D (R).
8π q(3.24)

57
Note that this new lowering (raising) operator contains two components: a raising (low†ering) operator D(i) (R) (D(i) (R)) between the ground state and the ith excited-state(i) †manifold due to the atom-laser interaction, and a lowering (raising) operator Aq eik·R(i)(Aq e−ik·R ) of type q corresponding to coupling with the reservoir photon field via a
photon with polarization q. Thus, the jump operator in the low-intensity equations
describes a transition cycle of the atom involving coupling to both the laser and thereservoir photon field. Note also that this new operator and the effective-Hamiltonian
term in the equation of motion are both proportional to the saturation parameter si ,the perturbation parameter.

3.3

Direct solutions of the semiclassical master equation
In Appendix C, the semiclassical master equation is derived in detail. This section

presents an overview of solution techniques of the resulting equations, called the optical
Bloch equations because of their similarity to the Bloch equations that describe theprecession of a nuclear magnetic moment in a magnetic field.
The optical Bloch equations for the multilevel atomic system considered here are,
h
i
³
´
h
i ³
´ i
(m)
(m)
ρ̃˙ ei , gj = i δω (m) − ∆ − kL v(i − j) ρ̃ ej , gi − eiωL t Ṽ , ρ̃  (m)  (3.25)~
ei ,gj
³
´
h
i ³
´ i
h
i
(m)
(m)
ρ̃˙ gi , ej
= i ∆ − δω (m) − kL v(i − j) ρ̃ gi , ej
− e−iωL t Ṽ , ρ̃  (m)  (3.26)~
gi ,ejhiiρ̃˙ (gi , gj ) = −ikL v(i − j)ρ̃ (gi , gj ) −Ṽ , ρ̃(3.27)~(gi ,gj )
³
´
h
i ³
´ ih
i
(m) (n)(m) (n)ρ̃˙ ei , ej
= i δω (m) − δω (n) − kL v(i − j) ρ̃ ei , ej−Ṽ , ρ̃  (m) (n)  .
~ei ,ej(3.28)
The above equations describe the evolution of the various matrix elements of the density
matrix, with the shorthand notation ρ(i, j) = hi|ρ|ji. The subscripts on the commutators likewise indicate a matrix element, [, ]i,j = hi|[, ]ji. The tildes over the density

58
matrix elements imply the following redefinition corresponding to optical coherences as
³
´
³
´
(m)
(m)
ρ̃ ei , gj → ρ̃ ei , gj eiωL t ,
³
´
³
´
(m)
(m)
ρ̃ gi , ej
→ ρ̃ gi , ej
e−iωL t .

(3.29)(3.30)For elements that are not optical coherences, the redefinition has no effect.
This system of equations can be solved by putting the various matrix elements ofρ̃ into a vector, corresponding to all of the possible combinations of N internal atomic
states and yielding a N 2 elements. The optical Bloch equations given above then provide
an N 2 × N 2 matrix, labeled here as Msc , for the time evolution of the ρ vector,ρ̃˙ = Msc ρ̃.

(3.31)

This equation can be solved by specifying an initial state for ρ̃ and numerically integrating the above differential matrix equation. For atomic systems with on the order of
10 to 100 internal states, this is usually numerically tractable.
If only the steady state of the system is desired, as is the case here in which
final laser-cooling temperature is the goal, it is faster to simply set ρ̃˙ = 0 in the above
equation and then solve for ρ̃. This involves a numerical inversion of the matrix Msc .
However, this matrix is singular and cannot be inverted. It turns out that the system
is overspecified, and this can be remedied by supplying the additional condition thatthe sum of the populations is always equal to unity, a fundamental requirement for thedensity matrix of a system with conserved norm. This condition can be enforced bymodifying one of the rows of Msc corresponding to a population element to have zeros inall locations except population elements. Then, the steady state solution to the matrix
equation is given byρ̃ss = M−1
sc v,(3.32)

where v is a vector with all zeros except for a 1 in the position corresponding to themodified row of Msc .

59
Once the steady-state density matrix of the system is obtained, any steady-state
properties of the system can be determined by performing traces over this density matrix
multiplied by the appropriate quantum mechanical operator Ôhi
O = hÔi = Tr Ôρ̃ss .

3.4(3.33)

Stochastic wave-function solutions of the fully quantum masterequation
The MCWF [9, 10, 11, 12, 13, 14, 15] technique is a means of interpreting asystem-reservoir master equation — which describes the evolution of a density operator
for a system interacting with a large external reservoir — as the evolution of an ensemble of individual wave functions, each undergoing random quantum jumps. The freeevolution of the stochastic wave functions is determined by the effective Hamiltonianthat was derived in section 3.2. The nature of the quantum jumps is determined by theleftover term in the master equation, which cannot be absorbed into the free-evolution
commutator. The components of this leftover term are often called quantum-jumpoperators.
The following discussion will deal primarily with the master equation in the lowintensity limit, as developed in the previous section, although the methods could justas easily be applied to the arbitrary-intensity master equation. The low-intensity limit,
however, provides a reduction in the number of internal atomic states required in thecalculation, and this will be beneficial for performing calculations later. Furthermore,since the lowest temperatures are achieved for low laser intensities, such a specialization
does not hinder the ability to calculate lower bounds of temperature.
Having already expressed the master equation in a form involving an effective
Hamiltonian and a jump term, the application of the MCWF technique is rather straightforward along the lines developed in the literature (see, in particular, Ref. [10]). For asingle stochastic wave function, the procedure is as follows. First, set the wave function

60
to an initial value. Then, numerically propagate the wave function for a time step δtaccording to the effective Hamiltonian Heff only, from an initial value |ψ(t)i to a final
value |ψ (1) (t + δt)i,
|ψ

(1)

µ
¶
iHeff δt(t + δt)i = 1 −|ψ(t)i .
~(3.34)

Restrictions on the size of δt are given such that the first-order truncation of the timeevolution operator in Eq. (3.34) is approximately valid. We note that Heff is nonHermitian by construction, as a result of absorbing parts of the relaxation operator into
the original (Hermitian) bare system Hamiltonian. Because of this, propagation with
Heff will not conserve the norm of the wave function when propagated to |ψ (1) (t + δt)i.
The time step δt of the propagation must be chosen so that δp ¿ 1 in the inner product,¯DE¯ψ (1) (t + δt) ¯ ψ (1) (t + δt) = 1 − δp.

(3.35)The quantity δp is the loss of norm resulting from propagating with Heff for a time step
δt, and is found to beδp = δt hψ(t)|X†

B(i) (R, k) · B(i) (R, k)|ψ(t)ii= δt hψ(t)|

X X
i

=

X X
i

†

Bq(i) (R, k)Bq(i) (R, k)|ψ(t)i(3.36)

q=0,±1

δpi,q .

q=0,±1

The total loss of norm has been decomposed into individual elements each corresponding
to a particular type of interaction with the reservoir (i.e., the q-value of the interaction,or the excited state i involved). These individual contributions are given by†δpi,q = δt hψ(t)|Bq(i) (R, k)Bq(i) (R, k)|ψ(t)i .

(3.37)

It is seen that the loss of norm due to a given type of interaction with the reservoir
is determined by the quantum-mechanical expectation value of the product of jump
operators of this type of interaction. The loss of norm δp can also be interpreted as theprobability for a quantum jump to occur.

61
After the wave function has been propagated as described above, and the values of
δpi,q calculated, it must then be determined whether or not a quantum jump occurred.
This is achieved by generating a pseudo-random number on a computer and comparing
it to the value of the total jump probability δp. If the random number is less than δp, aquantum jump occurred, and if it is greater, no quantum jump occurred. If a quantum
jump does occur, the type of quantum jump must also be calculated by comparing therandom number with the individual sub-probabilities δpi,q in the same manner.
If a quantum jump of type q, i occurs, the quantum jump lowering operator(i)

Bq (R, k) must be applied to the wave function from the beginning of the time step,sδt (i)|ψ(t + δt)i =B (R, k) |ψ(t)i .
(3.38)
δpi,q qThe square-root factor in front of the lowering operator is necessary for renormalization.
If no quantum jump occurs, then the stochastic wave function is simply renormalized.
The resulting wave function is then used as the starting point for propagation
over the next time step, and the procedure is repeated.
A good approximation of the true system density matrix is achieved by combining
the trajectories of a number of independently propagated stochastic wave functions, eachtrajectory having a unique sequence of pseudo-random numbers. (A thorough discussionof the statistical issues involved with the MCWF technique can be found in Ref. [10].)Once a suitable ensemble of stochastic wave function trajectories has been obtained, anestimate of the true expectation value of an operator is found by taking the ensemble
average of the expectation value of that operator with respect to the stochastic wave
functions. For example, an estimate of the average kinetic energy at a time t for a system
for which N independent stochastic wave functions have been calculated is given by
¯ 2 ¯
À
N ¿
¯P ¯
1 X
¯
¯
ψi (t) ¯ψi (t) ,hEi (t) =
N
2m ¯
i=1

where ψi (t) is the ith stochastic wave function, given at time t.

(3.39)

62
Figure 3.2 demonstrates a simple example of the application of the MCWF technique, wherein the average kinetic energy is calculated for a two-level atom interactingwith a one-dimensional standing-wave field. For this calculation, we have used a Rabi
frequency of Ω = Γ/2 and a detuning of δ = −Γ/2, where Γ is the decay rate of theupper to the lower atomic state, and we have set Γ = 400Erec , where Erec = ~2 k 2 /2m
is the recoil energy. The atomic kinetic energy, averaged over 500 stochastic wave functions each initialized to zero momentum, is plotted as a function of time, with error bars
indicating the error in the ensemble average for a given time. The separation of thetransient relaxation period from the steady-state is clear, the steady state regime being
characterized by fluctuations in the average energy about a mean. This noise is due tothe finite number of stochastic wave functions being used, and if a greater number ofwave functions were used, the amplitude of the fluctuations would be decreased. In thelimit of an infinite number of wave functions, the true density-matrix solution of themaster equation would be obtained. An estimate of the steady-state kinetic energy isfound by time-averaging the calculated data over the entire steady-state regime. Since
this is a larger ensemble than the set of wave functions for a single time, the error ofsuch an average will be smaller than the error bars shown in the figure.
As a test of the model, it is useful to perform a calculation for which the analytical solution is known. For a two-level atom moving in a one-dimensional field of
two opposing, linearly polarized lasers, it is expected that the steady-state temperature
will coincide with the Doppler limit when the detuning is set to it’s optimal value ofδ = −Γ/2, and when the intensity of the laser is vanishing, Ω → 0. This limiting temperature is achieved because there is no internal structure to the atom which would allow
sub-Doppler dynamics. Generally, in the literature, the Doppler limit will be quoted as
TD = ~Γ/2kB . This value of the Doppler temperature is useful as a rough estimate, and,indeed, most experiments are unable to sufficiently eliminate heating factors — such as
background-gas collisions and imperfections in the laser — to warrant such precision in

63<E>/Er

150

100

50

0

0

1000

2000

t ωr

3000

4000

5000

Figure 3.2: An example of a characteristic MCWF stochastic trajectory. Shown isthe average of the atomic center-of-mass energy over 500 independent stochastic wave
functions, as a function of time, for a two-level atom in a 1D standing-wave laser field.
The energy is given in units of the recoil energy Erec = ~2 k 2 /2m, and time is given inunits of the inverse recoil frequency ωr−1 = ~/Erec . All wave functions are initialized in
the ground state of the atom and localized in momentum space with zero momentum.
The steady state, wherein the system fluctuates around an average value, is seen to beachieved after a transient relaxation period. Error bars indicate the variance in the dataat each given time for the ensemble of 500 stochastic wave functions. An estimate of
the steady-state atomic center-of-mass energy is obtained by performing a time-averageover all wave functions for all times after the relaxation regime. The error bar of suchan average will be smaller than the error bars in the figure, which apply only to thedata for a given time. This figure is adapted from Ref. [2].

64
the first place. However, to check the Monte Carlo calculations, a more accurate value
will be needed, as is obtained in the 1989 paper by Castin, Wallis, and Dalibard [61], inwhich they survey the analytical understanding of Doppler atom-laser dynamics. Using
an innovative technique called the Family Method, they find that the one-dimensional
Doppler limit is in fact TD = 7~Γ/40kB . Employing a Monte Carlo wave-function calculation for this rather trivial system with δ = −Γ/2 and Γ = 400 ωrec , and performing
the calculation for a variety of laser intensities, the limiting behavior shown in Fig. 3.3
is found, with the dots indicating the numerical results of the Monte Carlo calculation,
and a quadratic fit to the data denoted by a solid line. In the low-intensity limit, the7temperature goes to TD = 40
~Γ/kB = 70 ~ωrec /kB , in agreement with the result ofCastin et al. [61].

3.5Calculations for 25 Mg and 87 Sr
The purpose of this section is to illustrate the application of the theory developed

up to this point to a complicated system. We wish to quantitatively study the dynamics
of particular atoms interacting with three-dimensional polarization-gradient laser fields.
The balance of the frictional cooling forces along with the diffusion experienced by theatom due to spontaneous emission and its interaction with the laser leads to a steadystate momentum distribution that determines the temperature of a gas of such atoms. In
particular, we will study here the cooling of the fermionic isotopes of two alkaline-earth
atoms, 25 Mg (nuclear spin I = 5/2, 1 S0 -1 P1 width Γ/2π = 81 MHz, hyperfine splittings ∆ω13 /2π = 46 MHz and ∆ω23 /2π = 27 MHz, where we have assumed a hyperfine
quadrupole parameter B = 0 [62]) and 87 Sr (I = 9/2, 1 S0 -1 P1 width Γ/2π = 32 MHz,hyperfine splittings ∆ω13 /2π = 43 MHz and ∆ω23 /2π = -17 MHz). These atoms, having nonzero nuclear magnetic moment, have degenerate (assuming zero magnetic field)
Zeeman sublevels. These sublevels allow for the mechanism of sub-Doppler cooling in
an appropriate laser configuration. Both 25 Mg and 87 Sr exhibit significant excited-state65

kB TÑΩrec
400
300
200
100
0.2

0.4

0.6

0.8

1

WWsatFigure 3.3: Simulation of a two-level atom with δ = −Γ/2 and Γ = 400 ωrec interactingwith a one-dimensional laser field, as a test of the Monte Carlo wave-function technique.
Plotted is kB multiplied by the steady-state temperature, in units of the recoil energy~ωrec . The dots are the numerical results of the Monte Carlo simulation, and the line
indicates a quadratic least-squares fit the data. The calculation yields a low-intensity
limit for the one-dimensional Doppler temperature in agreement with the known value
7~Γ/kB = 70 ~ωrec /kB .
TD = 40

66
spectral overlap, with ∆ω13 /Γ = 0.57, ∆ω23 /Γ = 0.33, and ∆ω13 /Γ = 1.3, ∆ω23 /Γ =-0.53, respectively. We consider the three-dimensional lin-⊥-lin laser configuration, consisting of a pair of opposing beams along each cartesian axis, in which each beam islinearly polarized orthogonal to its opposing beam. Furthermore, for this calculation,
we set to zero the relative phases of the three sets of laser pairs.
Having a nuclear spin of I = 5/2, the 1 S0 state of 25 Mg results in a hyperfine ground state with 6 sublevels. Use of the low-intensity master equation given in
Eq. (3.17) allows us to consider only these 6 internal states of the atom, since the excited states have been adiabatically eliminated in this regime. However, as noted inRef. [11], a momentum grid extending to 20~k in each direction with a spacing of ~kwould yield a density matrix with (6 × 413 )2 ≈ 2 × 1011 elements. A direct solution ofthis master equation is not numerically feasible, even without considering the further
increases in matrix size necessary to describe the master equation relaxation operator
in Liouville space [63]. On the other hand, the MCWF method only requires numericalpropagation of individual wave functions, which would be represented by vectors with
6 × 413 ≈ 4 × 105 elements. If the number of independent stochastic wave functions required to achieve satisfactory convergence for the calculation of a particular property ofthe system is not unreasonably large, the MCWF method provides a distinct advantage
over a direct master-equation solution.
We follow the procedure outlined in Section 3.4, working in the low-intensity limit
in order to reduce the number of internal atomic states in the calculation, which increases
the efficiency of calculation. Since laser cooling is most effective at low laser intensities,
this turns out to be a useful regime in which to work, with the additional benefit thatlower temperatures require a smaller number of atomic center-of-mass momentum statesin the calculation. We must determine the effective Hamiltonian as given in Eq. (3.22)
and the jump operators as given in Eqs. (3.23) and (3.24) for each atom, and for theparticular laser field being considered.

67
We consider here the lin-⊥-lin laser configuration in three dimensions, with therelative phases of the beams set to zero. The positive-frequency component of theelectric field isi
E0 h ikXŷe+ ẑe−ikX + ẑeikY + x̂e−ikY + x̂eikZ + ŷe−ikZ2
E0 X(−1)q aq (R)²̂−q ,=2

EL (R, t) =(3.40)

q=0,±1

with spherical coefficients´
1 ³
a+1 (R) = − √ e−ikY + eikZ + ieikX + ie−ikZ ,2´
1 ³
a−1 (R) = + √ e−ikY + eikZ − ieikX − ie−ikZ ,2
a0 (R) = e−ikX + eikY .

(3.41)(3.42)(3.43)

With these coefficients, along with parameters appropriate to the particular atom under consideration, the atom-laser raising and lowering operators given in Eqs. (3.14)
and (3.15) can be constructed. With knowledge of the effective Hamiltonian and theraising and lowering operators, we can then proceed with the MCWF procedure as
outlined.
Our example entails propagating 20 stochastic wave functions each for three different values of the light-shift parameter, ~|δ3 |s3 /(2Erec ) =10, 20, and 30, for both
25 Mg (I = 5/2) and 87 Sr (I = 9/2). We consider only δ

3 = −5Γ. As in Figure 3.2, we

calculate the stochastic trajectories of the ensemble average (i.e., averaged over the 20
wave functions) kinetic energy for each atom as a function of time. We continue thispropagation until the transient regime has been passed for some time, and use the time
average over the steady-state ensemble-average kinetic energy to provide an estimate of
the total average kinetic energy and the final error. The results are shown in Fig. 3.4,
along with the energies for atoms with an isolated cooling transition for comparison,
Je = Jg + 1 with Jg =1, 2, 3, and 4, with detuning δ = −5Γ, as first calculated by
Castin and Mølmer in Ref. [11]. From this cursory analysis, we can see that 25 Mg

68should exhibit a sharp rise in temperature with increasing laser intensity, while 87 Sr
will cool to sub-Doppler temperatures even for higher intensities, as has been noted
experimentally [52].
Detailed calculations of this sort, for realistic atoms, are quite computationallyexpensive. For example, a single data point for the Mg and Sr calculations presented
here required on the order of 200 hours wall time for a 20 processor parallel code,running on a cluster of 2.4 GHz Intel Zeon processors.
It is insightful to make a qualitative illustration of the effects of increasing theatomic degeneracy for an isolated transition as would be the case for a system having nohyperfine overlap. Using the semiclassical calculation method discussed in Section 3.3,
the force on such an atom can be determined as a function of the atom’s velocity.
Figure 3.5 shows plots of this force, in arbitrary units, as a function of the atom’s
velocity in units of the recoil velocity, in for increasing internal atomic degeneracy. The
top-left plot is for a system with a ground hyperfine quantum number of Fg = 0 to an
excited state of Fe = 1. For these parameters, no sub-Doppler cooling is expected, andthe force curve displayed is simply that of normal Doppler cooling. In the top-right plot,
with Fg = 1 and Fe = 2, the appearance of a small sub-Doppler force occurs. This force
usually appears in a small range of velocities around zero. Nevertheless, the force curveis steepened dramatically, and even over this small range can have a significant impact
on the effectiveness of cooling. The bottom-left plot, with Fg = 2 and Fe = 3 shows
an even more dramatic sub-Doppler force, and the bottom-right curve, with Fg = 3
and Fe = 4 even more so. This qualitative observation supports the more quantitative
results shown in Fig. 3.4 for isolated transitions.
As a final note, one of the advantages of performing detailed quantum calculations
is the ability to create illustrative pictures. Figure 3.6 provides one such example, using
the three-dimensional data obtained from the above Monte Carlo simulations in the raytracing rendering program Pov-Ray, and illustrating the reduction in three-momentum

69

200
Jg = 1
Jg = 2
Jg = 3150_ 2<(p / h k) >Jg = 4
25

Mg

87

Sr

100

50

0

0

5

10

15 _ 2025
30
h |δ|s / (2 Erec)

35

40

45

Figure 3.4: Results for calculated ensemble-average energies (rms momentum squared)for 25 Mg and 87 Sr, as a function of the light-shift parameter ~|δ3 |s3 /(2Erec ). For comparison, also shown are the calculated energies for atoms with isolated transitions,
Je = Jg + 1, with Jg =1, 2, 3, and 4, with detuning δ = −5Γ. This figure is adaptedfrom Ref. [2].

Force

70

0.03

0.03

0.02

0.02

0.01

0.01

0

0

-0.01

-0.01

-0.02

-0.02

Force

-0.03
-2
0.02

-1

0

1

-0.03
2
-2
0.02

0.01

0.01

0

0

-0.01

-0.01

-0.02
-2

-1

0
Velocity

1

-0.02
2
-2

-1

0

1

2

-1

0
Velocity

1

2

Figure 3.5: Semiclassical force curves illustrating the effect of increasing degeneracy
on the laser cooling of idealized atoms with isolated cooling transition and no overlapof hyperfine manifold. As the internal atomic degeneracy is increased and all otherparameters held fixed, the appearance of a sub-Doppler force grows. See text for details.

71
space of a cloud of 87 Sr atoms undergoing sub-Doppler polarization-gradient cooling.

3.6

Conclusion
In conclusion, this chapter has provided a detailed description of the fully quantum-

mechanical master equation that describes an atom with multiple internal internal structure interacting with a three-dimensional polarization-gradient laser field. It has beenshown how the spontaneous-emission relaxation operator is generalized for atoms of thistype. The MCWF technique has been applied to these equations of motion, providing
a more efficient means of performing calculations for these systems compared to a full
solution of the master equation. A few example calculations have been presented toillustrate the application of this theory to atomic systems interacting with laser configurations commonly used in experiments.

72

Figure 3.6: Three-dimensional renderings illustrating the time evolution of an atomic
cloud in momentum space, based on the numerical simulations of cooling 87 Sr.

Chapter 4
Bichromatic Cooling of Three-Level Atoms

4.1

Introduction
In this chapter, the subject of laser-cooling dynamics is further explored, thistime applied to a three-level atomic system. This type of atomic configuration is much
simpler than the elaborate multi-manifold systems discussed in Chapter 3. Moreover,instead of applying complicated polarization-gradient electromagnetic fields in three
dimensions, only two standing-wave laser beams in one dimension will be utilized here.
Nonetheless, this basic configuration yields a rich variety of novel phenomena with some
unexpected conclusions.
As in Chapter 3, the types of atoms that will be focused on will be group-II
atoms, the alkaline earths. Alkaline-earth atoms are difficult to laser cool to the µK
or nK regime, but this technique exhibits encouraging potential to circumvent current
roadblocks.
The work here also introduces a sparse-matrix technique that permits efficientsolution of the master equation for the stationary density matrix, including the quantized atomic momentum. This overcomes longstanding inefficiencies of exact solution
methods, and it sidesteps inaccuracies of frequently-implemented semiclassical approximations. The realistic theoretical limiting temperatures are optimized over the fullparameter space of detunings and intensities. A qualitative interpretation based on thephenomenon of electromagnetically-induced transparency reveals dynamical effects due

74
to photon-atom dressing interactions that generate non-Lorentzian lineshapes. Throughcoherent engineering of an asymmetric Fano-type profile, the temperature can be lowered down to the recoil limit range.
The inability to cool alkaline-earth atoms down to the 100 nK range has provento be a formidable bottleneck in the study of ultracold dilute gases. Novel coolingtechniques have limited applicability and so far quantum degeneracy has been obtained
only for ytterbium, a lanthanide atom with similar level structure, despite extensive
efforts by numerous groups [64, 65, 52, 2]. The theoretical description of laser cooling took strides in the 1990s, primarily utilizing the semiclassical approximation which
resulted in complicated calculations [45]. Few realistic calculations have treated theatomic motion quantum mechanically and predict realistic cooling temperatures. Here,a direct method is introduced, based on sparse-matrix techniques, which efficiently
yields the fully quantized stationary density matrix. In contrast to Monte Carlo approaches [10, 12, 2], the temperature and other observables emerge without statisticalmonitoring; the efficiency enables a complete mapping of the large parameter spaceof laser intensities and detunings. Here, this mapping is computed for a three-level
cooling scheme and develop a qualitative interpretation based on the physics of electromagnetically induced transparency (EIT). The calculations predict that promisingregimes exist that should produce record-low temperatures and overcome the difficultyof alkaline-earth atom cooling. This is orders of magnitude cooler than could be obtained using the main resonance line alone. Moreover, the cooling can be controlled andadjusted in real time, which should make this method competitive with or superior toother proposed ways of cooling alkaline-earth atoms.
Although leading to many advances [19, 66, 67], internal atomic structure greatly
restricts the species for which laser cooling can be applied successfully, and determines
whether ultracold temperatures — on the order of a few µK — can be achieved. Basic
Doppler cooling [18, 68] has a minimum temperature proportional to the width of the75
atomic transition. Lower temperatures can be achieved by exploiting the multilevel
hyperfine structure [16, 17, 69, 70], by modifying the atomic scattering [71], or by using
narrow optical transitions as a second-stage for pre-cooled atoms [64, 65]. Here, subDoppler cooling with three-level systems under two-color excitation is explored. The
published literature discusses some similar schemes, for both atoms [72, 73, 74] andtrapped ions [75, 76, 77]. In particular, the work of Tan et al. [73] appears initiallyto be the most exciting, since it claims to permit cooling to below the recoil temperature, believed to be a rigorous quantum lower bound for this type of laser cooling.
However, this is an artifact of the errors associated with a semiclassical treatment of
atomic momentum in Ref. [73], since that approximation is known to produce incorrect,
arbitrarily low (and even negative) temperatures. In order to realistically evaluate thepromise of three-level cooling based on dressed-state ideas, it is thus critical to carryout fully quantum solutions with great efficiency, as in the sparse solution techniqueintroduced here.
Parts of this chapter are developed from research by this author which has been
submitted for publication, and some of the text and figures have been adapted from thework therein. A preprint form of this submitted publication can be found in Ref. [3].

4.2

Three-level atomic systems
There are three basic types of three-level systems — the Λ, V, and Ξ systems —each classified according to the ordering of the bare quantum states in energy, and thepossible decay pathways. The Λ-configuration is commonly used for studying EIT [78]and related phenomena. The focus here is on a Ξ system, shown in Fig. 4.1, becauseof its relevance to current experimental work [79, 80]. Note, however, that the general
conclusions can be applied to any type of three-level system, and can be extended tosystems with more than three levels. The level structure and internal parameters for theΞ system are depicted in the left-hand side of the figure. This cooling technique seems

76
well suited to alkaline-earth atoms, which are good candidates for the next generationof optical atomic clocks, studies of ultracold collisions, optical Feshbach resonances [81],and quantum degeneracy [82].
The linewidths, Γi = γi /2π for i =1, 2, of both the upper and lower transitions
for some of the most relevant alkaline-earth (24 Mg, 40 Ca, and 88 Sr) and alkaline-earthlike (174 Yb)systems are given in Table. 4.1. From an initial analysis, it might seemthat mixing of the bare atomic states could allow cooling to a temperature set by thenarrowest of the two transitions. This seems reasonable because traditional Doppler
cooling with a single laser and two-level system predicts temperatures proportional to
the transition linewidth. Since the upper transition in all cases considered in Table. 4.1
is narrower than the lower transition, this simplistic analysis predicts that sub-Doppler
cooling should occur by a factor of Γ1 /Γ2 below the Doppler limit. This ratio is also
given in Table. 4.1, and it can be seen that, for example, magnesium offers the possibilityof beating the Doppler limit by the rather large factor of 36. Later in this chapter, thedetailed workings of this type of cooling will be explained, and it will be shown that
this limit is in fact an underestimate of effectiveness of the cooling, and the Doppler
limit can be beat by significantly more than a factor of 36.
Table 4.1: Transition linewidths, Γi = γi /2π for i =1, 2, of the 1 S0 -1 P1 -1 S0 three-level Ξsystems for common alkaline-earth atoms 24 Mg, 40 Ca, and 88 Sr, as well as the alkalineearth-like atom 174 Yb. The ratio of Γ1 /Γ2 is also shown, indicating the factor below thelower-transition Doppler limit which might be expected using a three-level bichromatic
scheme, neglecting the more detailed lineshape analysis provided in this chapter.
24 Mg

40 Ca

88 Sr

174 YbΓ1 (MHz)

79

34.63

32

28.01

Γ2 (MHz)

2.2

4.77

2.96

4.81

Γ1 /Γ2

36

7.3

11

5.8

77

|2 >
δ2

Γ1 ω2
δ1
ω1

Ω2

Ω1

|0 >

2

|+, n1 − 1 >1 |0, n1 >0−1|−, n1 − 1 >

−2
1

|1 >
Im(E)/h̄Γ1Γ2

Re(E)/h̄Γ13

0.8|−, n1 − 1 >

0.6
0.4

|+, n1 − 1 >0.2|0, n1 >
0
−2

−1

0

1

2

δ2 /Γ1Figure 4.1: Left side: Atomic configuration for generic three-level Ξ system. Right side:
The dressed eigenenergies for the 1 S0 -1 P1 -1 S0 Ξ system in 24 Mg with s1 (δ1 ) = 0.001 ands2 (δ2 ) = 1. Real (top) and imaginary (bottom) parts of the eigenvalues of Eq. (4.3),
with dressed atomic states labeled. The real parts are the energies and the imaginaryparts are the effective linewidths of the dressed atomic system. Both are plotted asfunctions of δ2 , with fixed δ1 = 0. This figure is adapted from Ref. [3].

78

4.3

Fully quantum model and its advantages
Fig. 4.1 shows the internal atomic states in order of increasing energy as |0i, |1i,

and |2i. The transition energy of the lower transition, |0i → |1i, and of the upper(1)transition, |1i → |2i, are ~ω0(2)

and ~ω0 , respectively. Two lasers are included, of

frequency ω1 and ω2 , and their detunings are defined from the appropriate atomic(i)transitions as δi = ωi −ω0 , for i = 1,2. The intensities of lasers 1 and 2 are characterized
by their Rabi frequencies Ω1 = − h0|d|1i·E1 (x) and Ω2 = − h1|d|2i·E2 (x), respectively,
where d is the electric-dipole operator of the atom and Ei is the electric-field amplitude
for the ith laser. The spontaneous-emission linewidths of states |1i and |2i are Γ1 /2π
and Γ2 /2π, (79 MHz and 2.2 MHz, respectively, for 24 Mg). The time evolution for theatom in the laser field, with mass m and center-of-mass momentum operator p, is givenby the master equation,ρ̇(t) =i[ρ, H] + L [ρ] ,
~

(4.1)

where ρ is the reduced density matrix of the atom system, the vacuum photon field
(1)

(2)

degrees of freedom having been traced over, H = p2 /2m + ~ω0 |1i h1| + ~ω0 |2i h2| +
(1)

(2)

Vlaser (x) + Vlaser (x) is the Hamiltonian of the atom-laser system, and the superoperator
Liouvillian L describes effects due to coupling of the atom to the vacuum photon field,
resulting in spontaneous emission [8].

4.4Numerical solution
Eq. (4.1) treats all the atomic degrees of freedom quantum mechanically, so its

solutions generally provide an accurate description of the atom’s dynamics. Thus, muchof the difficulty associated with semiclassical approximations of the system is avoided.
For reasonable temperatures, however, the number of numerical basis states required todirectly solve the problem, even in one dimension, is impractical for most computers.
Here, this problem is resolved without resorting to Monte Carlo methods by noting

79
that the matrix for the linear system equivalent to Eq. (4.1) is very sparse: only asmall fraction of its elements are nonzero. For a typical calculation, the Hilbert space
involves 3 internal degrees of freedom for the atom, and 150 momentum states in one
dimension. The resulting number of complex density matrix elements is thus 450×451/2
independent elements for this case owing to symmetry, and L would take over 160GB tostore if it was a full matrix. But the particular structure shared by Liouvillian operators
L that describe relaxation processes simplifies the numerics, and requires less than 2 GB
to solve the sparse system; specifically, the microscopic properties of atomic operators in
L permit construction of the matrix L with the zero elements eliminated. The steadystate solution of Eq. (4.1) is then found using a standard sparse-matrix inversion, giving
an exact direct solution of a fully-quantized master equation in just 1-2 minutes of CPU
time on a current workstation.
The steady-state density matrix has been calculated this way in one dimension.
­®
The average kinetic energy is p2 /2m = Tr(ρ p2 /2m), which in thermodynamic equilibrium can be equated with 12 kB T , where T is the temperature. The large parameter
space of the problem has been explored, e.g., varying the detunings, δ1 and δ2 , as well as
the strengths Ω1 and Ω2 of the two lasers. However, the lasers are better characterizedby how strongly they dress the atom. The saturation parameters s1 and s2 for therespective transitions are
si =

4.5

Ω2i
1
.
2 δi2 + (Γi /2)2

(4.2)

Results
The parameter space has three distinct regimes. In the first, with s1 & 1 and s2

arbitrary, only heating occurs, as expected from Doppler-cooling theory since the lower
transition is driven strongly. Numerical results for the second regime, with s1 ¿ 1
and s2 ¿ 1, are shown in Fig. 4.2, with the temperature normalized to the lower(1D)level Doppler limit TD

(1D)

, TD

= 7~Γ1 /40kB [16], which is the optimum temperature

80
expected for cooling with just one laser. In this plot, the location of the bare two-photon
resonance (δ1 + δ2 = 0) is denoted by the dashed line. It is clear from this plot (note thenormalization of the temperature axis) that cooling occurs only to the Doppler limit for(1D)

the lower transition TD= ~Γ1 /2kB , and only in the range near δ1 = −Γ1 /2. In thiscase, laser 2 has no effect, as this amounts to simple two-level Doppler cooling on thelower transition with laser 1.
In the final regime, when s1 ¿ 1 and s2 & 1, cooling occurs down to substantially(1D)

below TD

. Figure 4.3 illustrates cooling in this regime, also for 24 Mg, with thetemperature again normalized to the lower-level Doppler limit and the bare two-photon
resonance indicated. The steady-state temperature is plotted as a function of δ1 andδ2 for s1 (δ1 ) = 0.001 in both the upper and lower plots, and with s2 (δ2 ) = 1 in theupper plot and s2 (δ2 ) = 5 in the lower plot. Note that the saturation parameters arebeing held fixed as the detunings are varied, so that the Rabi frequencies are being(1D)

continuously adjusted. The lowest temperatures, on the order of 10−2 TD, are foundin the quadrant with δ1 > 0 and δ2 < 0, as well as less extreme cooling in other regions.
Observe that the lowest temperatures are obtained for frequencies detuned to the blue ofthe two-photon resonance. This seems counterintuitive, since a red detuning is usuallyrequired in order to have a net decrease of atomic kinetic energy in a photon-scattering
event.

4.6

EIT explanation
Qualitative understanding of the cooling mechanism emerges from analysis of thesimpler Hamiltonian,




Ω1
0

0

~


H = Ω∗1 −2δ1 − iΓ1
.
Ω2

2


0
Ω∗2
−2(δ1 + δ2 )iΓ2

(4.3)

81

Figure 4.2: Steady-state temperatures for bichromatic three-level laser cooling, as afunction of both detunings δ1 and δ2 , with both atomic transitions perturbatively
probed, s1 ¿ 1 and s2 ¿ 1. The bare two-photon resonance is indicated with adashed line.

82

Figure 4.3: Steady-state temperatures for bichromatic three-level laser cooling, as afunction of both detunings δ1 and δ2 , with the lower atomic transition perturbatively
probed, s1 ¿ 1, for both plots, and the upper atomic transition dressed to varying
degrees in the two plots, with s2 = 1 in the upper plot and s2 = 5 in the lower plot.
The bare two-photon resonance is indicated with a dashed line. This figure is adaptedfrom Ref. [3].

83
Its complex eigenvalues have real dressed energies, and imaginary parts giving dressedstate linewidths. These dressed energies and widths are plotted on the right side ofFig. 4.1 as functions of δ2 , for the same parameters used in Fig. 4.3. The cancelationof one of the widths can be viewed as an EIT effect: since laser 1 is perturbative whilelaser 2 strongly dresses the upper transition, the new eigenstates of the system, denoted|+i and |−i, are well approximated as linear combinations of the bare states |1i and |2i.
These states have the modified energies and widths shown in Fig. 4.1. The linewidth
modifications can be viewed as a Fano interference [34], in which the dressing-lasertransitions caused by the probe laser enable multiple coherent pathways among thebare states. Constructive or destructive interference respectively increases or decreases
the atomic linewidth.
The cooling mechanism is thus qualitatively explained as ordinary two-level Doppler
cooling. But instead of using a transition between two bare states of an atom, the transition occurs between a (mostly) unmodified ground state, and a dressed excited state,with a shifted energy and a new linewidth that can be narrower than the bare linewidthof the lower transition. As the probe laser is scanned, the detuning relative to thedressed energy levels is varied, but since these levels are shifted from their bare energies, resonance occurs for different detunings than are encountered in the bare system.
In fact, the shifts of the eigen-energies in Fig. 4.1 from the bare energies explain theapparent observation of blue two-photon cooling in Fig. 4.3. In the dressed system, thebare two-photon resonance is no longer meaningful, and the cooling region is in fact tothe red of a dressed resonance.

4.7

Atom-photon dynamics for asymmetric lineshapes
When mapping this system onto Doppler cooling theory, note that the lineshapes

are not Lorentzian, but are asymmetric Fano lineshapes for the dressed system, as
illustrated by the solid line in the upper part of Fig. 4.4 as a function of δ1 for a fixed84
of δ2 = −Γ1 /2 [74, 3]. This changes the optimum-detuning condition to be that the
maximum cooling for a given transition occurs when the probe laser is detuned fromthe dressed excited state precisely to the inflection point of the absorption spectrum.
This can be understood by noting that the applied force f due to the laser beam is
proportional to the absorption rate, for a given δ1 . As in semiclassical cooling theories,
the friction coefficient α for the atom in the laser field isα=−df (v),
dv(4.4)

where v is the atomic velocity. Since the detuning of the laser and the resonant atomic
velocity are linearly related, the derivative of the absorption spectrum with respectto δ1 also yields a maximum in the cooling force. This is evident in normal Doppler
cooling because the optimum detuning occurs when δ = −Γ/2, the inflection point of
the Lorentzian. In general then, for asymmetric lineshapes, the optimum detuning does
not obey such a simple relation, but depends on the degree of asymmetry.
This point is illustrated in the upper portion of Fig. 4.4, which compares the actual
asymmetric absorption spectrum (solid line) of the dressed system as a function of δ1 andfor fixed δ2 = −Γ1 , with the more simplistic predictions for the absorption spectrum
using symmetric Lorentzian lineshapes (dotted line). The peak of each lineshape islocated at a dressed eigenenergy, and lineshapes at the same energies and widths. As an
example, optimum-cooling detunings relative to the leftmost resonance are illustrated
with a black arrow for the true asymmetric lineshape and with a gray arrow for the
hypothetical symmetric-lineshape case. The value of the optimum detuning, as well as
the slope of the lineshape is seen to be different for these two cases.
From this complete picture of the cooling mechanism, the minimum temperatures
can now be determined, allowing for the detuning modification due to asymmetric lineshapes. The lower part of Fig. 4.4 shows the ratio of the maximum slope of a Lorentzianlineshape with width Γ1 to the slope of the asymmetric lineshape, as a function of δ1

Absorption Rate

85

−2

−1

0

1

2

0

1

2

δ1 /Γ1

0

T /TD

10

−1

10

Γ2 /Γ1

−2

10

−2

−1

δ1 /Γ1Figure 4.4: Upper plot: Comparison of the true absorption spectrum of the dressed
three-level Ξ system (solid line) with a simplistic absorption spectrum with Lorentzianlineshapes (dotted line). Optimum laser-cooling detuning is indicated with arrows foreach type of spectrum. Lower plot: Comparison of the numerical three-level-coolingtemperature results (data points) with the ratio of the maximum slope of a Lorentzianlineshape with width Γ1 to the slope of the asymmetric lineshape, as a function of δ1with δ2 = −Γ1 /2 (solid line). This ratio provides an indication of the expected coolingfor the dressed system relative to the Doppler limit for the lower transition. This figure
is adapted from Ref. [3].

86with δ2 = −Γ1 /2. This ratio provides an indication of the expected cooling for the
dressed system relative to the Doppler limit for the lower transition. For comparison,
fully quantum numerical results are indicated by data points. Note that the expected
temperature, due to the asymmetric lineshape, is predicted to be lower than the uppertransition Doppler limit, indicated by the dotted line in the lower plot of Fig. 4.4, aprediction supported by the numerical data.

4.8

Conclusion
In conclusion, coherent engineering of a three-level system can optimize the ef-

fectiveness of two-level Doppler cooling.

Dressed states are created with modified

linewidths in the range between the smallest and the largest of the two bare linewidths,
and the additional effect of asymmetric lineshapes can lead to temperatures below the
Doppler limit of either bare transition. The ability to tailor the degree of cooling lends
this technique additional utility, and may be particularly useful when applied the notoriously difficult-to-cool alkaline-earth atoms. A dressing scheme can be suited to thecharacteristics of a particular atom, and real-time adjustment of the cooling properties
can allow narrowing of the velocity-capture range as an atomic gas is cooled. Utilizing
such coherent effects should lead to simple schemes for cooling far below the typical
Doppler limit.

Chapter 5
Cavity Cooling of Atoms and Molecules

5.1

Introduction
In this chapter, a cooling scheme is explored which, in addition to an atom (ormolecule) and a laser, utilizes an optical cavity. The idea of using a cavity to help cool
atoms was first proposed in the early 1990’s [83], and early experimental explorations [84]confirmed the feasibility of the method. Further investigation [71, 85] found situations
in which the cooling was much more effective than expected, reaching temperatures far
below the predictions of the basic theory. Further theoretical investigations [86, 87] wereable to attribute the excess cooling to the presence of so-called collective effects: basically, phenomena occurring due to the presence of multiple atoms, all interacting withthe cavity at the same time. This was not really a many-body effect, as the collective
enhancement could be modeled without including any atomic interaction physics.
The previous work in this field has led to a good understanding of cavity cooling applied to atoms. However, the cooling scheme, which stipulates driving an atomwith a laser having a frequency very far detuned from the atomic transition, holds thepossibility of being effective also for molecules, which are traditionally not amenable tolaser-cooling schemes. Almost all types of typical laser cooling require a closed internaltransition in the object being cooled — the excited state must decay back to the original
ground state in order for the process to be repeated — and molecules with their dense
forest of energy levels usually do not provide this. Cavity cooling, being far detuned so

88that coherent-scattering processes are dominant, offers a means of avoiding this problem. Additionally, the collective effects that were found to enhance cavity cooling ofatoms may also bolster the cavity cooling of molecules, leading to even more effective
cooling than might be expected.
The study of the cavity cooling of molecules offers another opportunity to apply the quantum-optics master-equation techniques developed throughout this thesis,
and, as it turns out, another application of the stochastic wave-function techniques developed in Chapter 3. In this chapter, the fundamental picture of the mechanism forcavity cooling is presented, as well as the particular details, such as Raman transitions,that must be considered when considering the cooling of molecules instead of atoms.
Following this is the development of a fully quantum-mechanical theory of a particle,
driven by a laser, and coupled to a high-finesse optical cavity, with both the atom andthe cavity coupled to the vacuum photon field. The methods of system-reservoir theory
are clearly useful again here, but with some generalization to accommodate the more
complex situation at hand. Certain approximations are introduced to make the problem more amenable to numerical solution, and their validity is discussed. Next, it isshown how the resulting master equation can be simulated using the same stochastic
wave-function framework developed in Chapter 3 to describe the dynamics of multilevelatoms. Finally, some results are shown.
In the end, the technique of solving a fully quantized master equation proves
mostly unfeasible for describing collective effects. Despite the greatly enhanced calculation possibilities provided by the stochastic wave-function methods, the large numberof particles required to model the desired collective effects is simply too large, withthe existing semiclassical calculations [86, 87, 88] having a distinct advantage in thisrespect. The fully quantum methods presented here, however, allow for a more detailed
description of the cooling process, as well as easily describing the time-evolution of thecooling. This turns out to be important in aiding the design of experiments.

89

5.2

Physical picture of the cooling mechanism
This section discusses the basic physical mechanism by which cavity laser cooling

occurs. Atoms and molecules will be considered, both treated as two-level systems, butwith the difference that the two-level molecule can decay from its excited state either
back to the ground state or to one or more other internal states. The latter process
is called a Raman transition, and it is assumed here that a molecule undergoing sucha transition becomes no longer trapped by the experiment, and the molecule is lost.
These other internal atomic states are not included in the model explicitly, but thepossibility of a Raman transition occurring and the molecule being lost will be built in.
The branching ratio for either going back to the ground state or undergoing a Raman
transition is determined by the various internal oscillator strengths for the molecule,and it can be adjusted in the model to suit the molecule being considered.
When an atom interacts with a laser, photons are absorbed and scattered. In
general, two different processes can occur, incoherent and coherent scattering. Incoherent scattering dominates at higher laser intensities which saturate the atom. The
absorbing photons arrive often enough that, by the energy-time uncertainty principle,
the scattered photons do not necessarily have the same energy as the incoming photons.
If, say, two photons are absorbed and scattered in a short enough time, only the sumof the frequencies of the outgoing photons must equal the sum of the two absorbed
photons. Coherent scattering on the other hand, dominates at low laser intensities forwhich the atom is only weakly perturbed. Photons arrive spaced apart far enough intime that the scattered photons must have the same frequency of the incoming photons.
By detuning far from the atomic transition, and using an appropriately small
Rabi frequency, that is, the intensity of the probe laser, the incoherent scattering peaks
are suppressed, and coherent scattering becomes the dominant process. However, coherent Raman transitions to untrapped states remain. In free space, the ratio of the90rate of Rayleigh (good) scattering events to the rate of Raman (bad) scattering eventsis independent of the detuning, and depends only on the various molecular oscillator
strengths. Introducing a cavity can, however, enhance one process over the other, forthe right parameter values. Momentum and energy conservation will be imposed forthe scattering process, and the possibility of Raman transitions will be included.
A molecule is in a cavity oriented along the z-axis and interacts with a transverse
probe laser aligned along the x-axis. The laser frequency is ωL = ckL and the cavity frequency is ωc = ckc . The molecule scatters a photon with wave vector ks . Conservationof momentum requires
p0 = p + ~kL − ~ks ,(5.1)

where p = mv, and v is the velocity of the molecule, the prime referring to the state
after scattering. Energy conservation requiresωs = ωL + ∆ + δr ,

(5.2)

where
∆ = −(kL − ks ) · v −~(kL − ks )2 ,
2m

(5.3)

and δr = ωr − ωa , is the internal energy difference for the Raman transition. The firstterm is the so-called two-photon Doppler effect, and the second term describes recoil
heating due to laser and cavity photons. Rayleigh-scattering processes are described by
the above equations with δr = 0.
Equation (5.2) indicates that the scattered photon has a frequency determinedby the probe-laser frequency, the angle of the scattered photon relative to the initial
velocity of the molecule, and the Raman energy difference. Since the motion along
the x-axis is not of importance here, it will be convenient to set vx = 0. Defining
~k2

2

s
ωrec,L = 2mL as the recoil frequency associated a laser photon and ωrec,s = ~k2m as therecoil frequency associated with a cavity photon, Eq. (5.3) then becomes∆ = ks · v − ωrec,L − ωrec,s .

(5.4)

91
If the scattered photon has a frequency greater than the frequency of the laser photons,
ωs > ωL , then Eq. (5.4) implies that ∆ must be positive for a Rayleigh scattering event(δr = 0), which implies that ks · v > ωrec,L + ωrec,s . This in turn implies that ks · v is
positive, so that the vectors ks and v must point in the same direction along the z-axis.
With this information, it is clear from Eq. (5.1) that the final molecular momentum
along the z-axis p0z must be less than the initial molecular momentum along the z-axis
pz . Thus, for a scattering event in which ωs > ωL and in which the scattered photon isemitted in the z-direction, the molecule’s momentum is decreased.
Up to this point, none of this is terribly useful for cooling since scattered photons
can of course be emitted in directions other than along the z-axis and with frequencies
that are not greater than the laser frequency. However, if an optical cavity with frequency ωc is introduced along the z-axis, it will greatly enhance photon emission into
the cavity mode from an atom or molecule in the cavity. If the cavity frequency is set
such that ωc > ωL , then from the discussion in the preceding paragraph, scattering processes which decrease the particle’s momentum along the z-axis will occur preferentially,
more often than processes that increase its momentum. The effect of this on a gas of
atoms or molecules over repeated scattering events is to cool the gas.
In a cooling experiment, Raman transitions are undesirable since they lead to aloss of molecules. The process of cavity cooling acts to suppress Raman processes. The
enhancement by the cavity of emissions along the z-axis with frequency ωc implies in
Eq. (5.2) that if Rayleigh processes (δr = 0) are resonant with the cavity, then Ramanprocesses (δr 6= 0) will not be, and will not be enhanced.

5.3

Theoretical models
For the purposes of developing the basic physical model, only atoms and notmolecules will be considered for the time being. As before, the atom will be approximated as a simple two-level system. Later, once the underlying theory has been de-92veloped, the behavior of molecules will be integrated into the equations. The modelis intended to describe a situation in which one or more two-level atoms in a cavity,interacting with a single cavity mode along the z-axis, and driven by a transverse pumplaser along the x-axis. The equations of motion will first be developed for a single atom
in a cavity. After that, they will be extended to treat the case of multiple atoms, eachinteracting with the cavity but neglecting interactions between atoms. Master-equation
techniques will be used to describe the fully quantum system. The cavity mode will thenbe approximated as a classical field, a valid approximation in the limit of large numbersof cavity photons. This will keep the number of quantum states needed for a numericalcalculation to a minimum when the cavity is well populated. If the cavity populationis very small and the semiclassical approximation is not valid, the fully quantum equations will remain numerically tractable since only a few modes of the cavity will need
to be included in a calculation. Stochastic wave function simulations similar to those
employed in Chapter 3 will be applied to perform calculations.

5.3.1

Single atom

The master equation describing an atom, coupled to the photon vacuum, in acavity, also coupled to the photon vacuum, with the atom driven transversely by alaser, isρ̇ =i[ρ, Hs ] + Lρ + Kρ,
~

(5.5)

where ρ is the reduced density matrix of the atom-cavity-mode system, L is a superoperator describing spontaneous emission (characterized by the parameter Γ), K is asuperoperator describing damping of the cavity mode (characterized by the parameterκ), and the system Hamiltonian is given byhi
Hs = −~δc a† a − ~δa σ † σ + ~ σ † (g(x)a + Ω) + H.c. ,

(5.6)

93
where a (a† ) is the annihilation (creation) operator for a cavity photon, σ (σ † ) is theannihilation (creation) operator for an internal atomic excitation, δc = ωl − ωc is thelaser-cavity detuning, and δa = ωl − ω0 is the laser-atom detuning.
The equation of motion for the cavity mode can be determined by first writing
down the Heisenberg equation of motion for the cavity-mode operator. This will yielda hybrid set of equations, with the atom-laser system being described by a densityoperator equation of motion (master equation) and the cavity mode described by anoperator equation of motion with fluctuation terms (Langevin equation), in which form
it can be more easily be approximated as a classical field. The process of transitioning
from a master equation to a Langevin equation is discussed in the literature [58, 15, 89].
Skipping these details, the resulting equation of motion for the cavity mode isκi[Hs , a] − a + F1~
2³
κ´
= iδc −
a − ig(x)σ + F1 .
2

ȧ =

(5.7)

To treat the cavity mode classically, it is necessary to find the equation of motion for
the average value of the mode operator a, which is denoted here by α = hai. Utilizing
the characteristic of noise operators that they average to zero, hF1 i = 0, a set of two
coupled equations is obtained, the first of which is³κ´
α − ig(x)σ,α̇ = iδc −
2

(5.8)

and the second of which isρ̇ =i[ρ, Hs ] + Lρ,
~

(5.9)

with the density matrix now only consisting of internal and external atomic states (thecavity states having been traced over), and with the system Hamiltonian given now byhi
Hs = −~δc α2 − ~δa σ † σ + ~ σ † (g(x)α + Ω) + H.c. .

(5.10)

Note that in obtaining Eq. (5.9), the fact that Trcav Kρ = 0 has been used, which isa general property of superoperators of Lindblad form when the states to which they

94
correspond are traced over. The number of photons n in the cavity mode is given byn = α2 . The solution of this system of equations should agree with the solution of
the fully quantum density matrix equation in the limit that the cavity mode is well
populated.
Given that only situations in which the pump laser will drive the atom weakly
are being considered here, it is possible to make the further approximation that theexcited state of the atom can adiabatically eliminated. This will decrease the state
space required for the calculation by a factor of two, thus improving the efficiency. This
will be especially necessary when considering more than one atom at a time, which is
expected to be computationally intensive.
The adiabatic-elimination procedure is easiest to implement using projectionoperator techniques. The goal is to break Eq. (5.9) into four separate equations, one eachfor the sub-density-matrices ρij = Pi ρPj , for i,j = g,e, and where Pi is the projection
operator into the i subspace. These sub-equations are³´ρ̇ee = −Γρee + i [g(x)α + Ω] ρeg σ − σ † ρge
¶
µ
³
´
Γ
− iδa ρeg + i [g(x)α + Ω] ρee σ † − σ † ρgg
ρ̇eg = −2¶µΓ+ iδa ρge − i [g(x)α + Ω] (σρee − ρgg σ)ρ̇ge = −2³´
ρ̇gg = Γρee + i [g(x)α + Ω] ρge σ † − σρeg .

(5.11)(5.12)(5.13)(5.14)The next step is to set ρ̇ee = ρ̇eg = ρ̇ge = 0 (detailed explanations of this procedure canbe found, for example, in Refs [8] and [60]). The term ρee in the two optical coherence
equations can also be neglected (again, see [60]) to obtain
ρeg =g(x)α + Ω †σ ρgg
δa + i Γ2(5.15)

ρge =g(x)α + Ωρgg σ.
δa − i Γ2(5.16)

The condition ρ̇ee = 0 can now be imposed on the ρee equation, and, substituting in the95values of the optical coherences just obtained,
#
"
i
1
1
2
ρee = [g(x)α + Ω]−σ † ρgg σΓδa + i Γ2δa − i Γ2=[g(x)α + Ω]2 †¡ ¢2 σ ρgg σ
δa2 + Γ2

(5.17)= sσ † ρgg σ,
where the saturation parameter of the atom due to interaction with both the laser andthe cavity field has been defined,s=[g(x)α + Ω]2¡ ¢2 .
δa2 + Γ2

(5.18)

Inserting the expression for ρee along with the previously obtained optical-coherence
expressions into the ground-state equation of motion,#"†
σσ † ρgg
2 ρgg σσ†−ρ̇gg = Γsσ ρgg σ + i [g(x)α + Ω]δa − i Γ2
δa + i Γ2#"
δa + i Γ2δa − i Γ22†††= Γsσ ρgg σ + i [g(x)α + Ω]
¡ ¢2 ρgg σσ −¡ ¢2 σσ ρgg
δa2 + Γ2
δa2 + Γ2hi Γ no= iδa s ρgg , σσ † − s ρgg , σσ † + Γsσ † ρgg σ,
2(5.19)

where the square brackets denote a commutator while the curly brackets denote an
anticommutator. Note that the last two terms comprise a decay operator of Lindbladform. This is encouraging, since the adiabatic elimination applied to Lindblad-typesystem should yield another Lindblad-type system. Note also that all terms are of firstorder in the perturbation parameter s. This makes sense because the above adiabaticelimination procedure is equivalent to expanding the equations of motion in powersof the intensity of the laser driving transitions within the atom, then truncating thisexpansion and keeping only terms up to first order.
The ground-state equation of motion that has been obtained can now be solved
using the same wave-function Monte Carlo technique described in Chapter 3, whereat each time step the other coupled equation of the cavity mode amplitude α found96
above must be consistently solved. Note that, as opposed to the atom-laser system
studied in Chapter 3, the equations here consist of two separate entities coupled to thephoton reservoir. These must each be treated separately, with appropriate quantum
jump operators. In any given time step, it is possible for either or both the atom andthe cavity to undergo a quantum jump. The jump operator for the atom is identicalto the operator used in Chapter 3, except that it is much simpler since the atom here
is only a two-level system. The jump operator for the cavity causes the cavity mode,
which is a quantum harmonic oscillator, to decay down the ladder by one step. The
properties of this operator are easily obtained from the atomic case, with the lowering
operator a playing the role of the internal atomic lowering operator (i.e., Pauli spin
matrix) σ.

5.3.2

Multiple atoms

For N atoms, with each atom interacting individually with the cavity, and neglecting interactions between different atoms, the system Hamiltonian becomes
Hs = −~δc a† a +

N
X

Hi ,

(5.20)

i=1

where Hi is the Hamiltonian describing a single atom and its interaction with the pumplaser,hi
Hi = −~δa σi† σi + ~ σi† (ga + Ω) + H.c. .
5.3.3(5.21)

Accounting for Raman losses and other experimental effects

In order to extend the model just developed to describe Raman transitions that
occur in molecules, the following procedure can be employed during the Monte Carlowave-function procedure. At the end of each time step, the saturation of the atom given
in Eq. (5.18), due to the atom’s interaction with both the laser and the cavity mode,can be calculated. Then, using a random number generated in the simulation, along97
with the knowledge of the branching ratio for Rayleigh to Raman scattering processes
obtained from details of the internal structure of the molecule, it can be determined
whether or not a Raman transition has occurred. The probability of this occurrence
is proportional to the value of the saturation parameter. If a Raman transition does
occur, the molecule can be removed from the simulation.
The detailed time-dependent calculations that the method presented here yields
also allow for the ability to include realistic experimental effects that simpler calculations
must ignore. This can be useful when the success of an experiment depends on factors
such as the loading of a pulse of molecules into the cavity region, or the escape of fastmoving molecules out of the cavity region. In principle, it is possible to include all ofthese features in the Monte Carlo procedure. For example, using the many-atom model,
the simulation could begin with a certain number of atoms of a given temperature. These
atoms then begin cooling, while new atoms can be introduced as time evolves.

5.4

Comparison of fully quantum to semiclassical model
The validity of the classical-cavity-mode approximation made earlier can be tested

by performing Monte Carlo simulations using both the fully quantum equations of motion and the semiclassical equations of motion. A simulation is set up in which an
ensemble of atoms is cooled using cavity laser cooling. The possibility of Raman transitions is neglected for the present case. The parameters for the simulation are Ω = 5ωr ,g = ωr , η = 20ωr , δa = −10ωr , δc = −20ωr , Γ = 2ωr , κ = 5ωr , and an initial temperature of around T = 50~ωr /kB . The calculation utilizes the single-atom model developed
above, and propagates 25 stochastic wave functions. The cavity-laser detuning δc is held
fixed.
The ensemble-average population distribution for a series of times during the
propagation is shown in Fig. 5.1. From these plots, it can be seen how the coolingaffects only the population with velocities (both positive and negative) that satisfy the

98resonance condition in Eq. (5.2), and also including velocities in a range of ±κ/2kL
around that point. Population around p = ±8~kL is shifted into lower momentum
states as the scattering proceeds, leading to a buildup of population in the middle ofthe distribution. The cavity-laser detuning could then be tuned slowly as the population
shifted down, squeezing the ensemble towards the center of the graph. Thus, it would
be expected that the limiting temperature should be on the order of κ, which is inaccordance with other cooling theories. This limit has been verified by the simulations.
The average kinetic energy over time for the same calculation is shown in Fig. 5.2.
The black line represents the calculation performed using the fully quantum model, while
the red line represents the calculation performed using the model with the classicalcavity-mode approximation. The results are in good agreement within the error bars of
the calculation, which used only 25 stochastic wave functions. This indicates that thenumerically faster semiclassical model provides a good alternative to the slower fully
quantum model.
Figure 5.3 illustrates the evolution of the population of the cavity mode over time
for the same calculation. The black line is again obtained using the fully quantum model,
and the quantity calculated is the quantum-mechanical average of the cavity photon
number operator n = a† a. And again, the red line is obtained using the semiclassicalmodel, with the quantity calculated being the average of the square of the cavity-mode
amplitude α2 . The fully quantum calculation was performed including only the first 5
cavity states (n = 0 - 4) in the numerical basis. It can be observed in the figure that theresults do not agree well in the short term, with the semiclassical calculation oscillatingwith much greater amplitude than the fully quantum calculation, before the two settle
down in agreement at a value of 1. This is not surprising, though, since a cavity
population of unity is hardly in the high-population limit for which the classical-mode
approximation is valid. Despite this discrepancy, the evolution of the atomic momentum
distribution in Fig. 5.2 agrees quite well. This is encouraging, as the important quantity

99

0.09
0.08

t=1

t=0

0.07

Probability

0.06
0.05
0.04
0.03
0.02
0.01
0
−20

−10

0

10

−20
20

−10

p/h̄k

0

10

20

p/h̄k

0.09
0.08

t=3

t=2

0.07

Probability

0.06
0.05
0.04
0.03
0.02
0.01
0
−20

−10

0

10

−20
20

−10

p/h̄k

0

10

20

p/h̄k

0.09
0.08

t=5

t=4

0.07

Probability

0.06
0.05
0.04
0.03
0.02
0.01
0
−20

−10

0

p/h̄k

10

−20
20

−10

0

10

20

p/h̄k

Figure 5.1: Time evolution of the momentum probability distribution for an ensembleof particles undergoing cavity laser cooling. Time is indicated, measured in units of theinverse recoil frequency ωr−1 .

100

50
49
48

Ē/h̄ωr

47
46
45
44
43
42
41
40
0

5

10

15

20

25

tωr

Figure 5.2: Time evolution of the average kinetic energy of an ensemble of particles
undergoing cavity laser cooling. Energy is given in units of the recoil energy ~ωr andtime is measured in units of the inverse recoil frequency ωr−1 . The black line was obtainedusing the fully quantum model while the red line was obtained using the semiclassicalmodel.

101
for studying cooling processes is the momentum distribution, and not the cavity-mode
properties.

5.5

Conclusion
This section has explored an interesting application of the master equation tech-niques built up elsewhere in this thesis. The generalization of the equations of motion
for an atom and a laser, to include a cavity, adds a great deal of richness to the phenomena that can be explored. Extending the stochastic wave-function methods first applied
to an atom to a cavity interacting with the photon vacuum provided another exampleof the utility of such techniques. In the end, however, the number of atoms required to
explore the collective effects as originally intended proved to be out of reach of thesemethods, given the computing speeds of the present. More efficient models that treat
the center of mass of the atom classically have had more success in observing collectivecooling phenomena. The main practical use of the calculation methods presented here
has been in providing realistic modeling of experimental details.

102

3.5
3

­ † ®
a a

2.5
2
1.5
1
0.5
0
0

5

10

15

20

25

tωr

Figure 5.3: Time evolution of the cavity-mode population. The black line was obtainedusing the fully quantum model, calculating the quantum-mechanical average of the cavity photon number operator n = a† a. The red line was obtained using the semiclassicalmodel, calculating the average of the square of the cavity-mode amplitude α2 . Time ismeasured in units of the inverse recoil frequency ωr−1 .

Chapter 6
Nonlinear Optical Spectroscopy

6.1

Introduction
The goal of this chapter is to bring together some of the ideas developed in thisthesis and apply them to the subject of nonlinear optical spectroscopy. The problem
is formulated in a rather general manner, while the eventual aim is to provide a theoretical analysis of the physics of collisional effects on atomic lineshapes. An important
application of this theory is the investigation of the experimental technique of transient
four-wave mixing, in which a series of short, temporally spaced laser pulses are applied
to an atomic or molecular gas or to a condensed-matter system. This setup can provide
a great deal of information about the interaction properties of the system. Importantly,
details of collisions within the system can be probed on a time scale for which theMarkov approximation, which has been utilized in Chapter 3 and Appendices B and Cof this thesis, is no longer valid. The physics of these short time scales, for which theinteraction properties are affected by the system’s memory of previous events, is referred
to as non-Markovian dynamics.
The current state of ultrafast laser technology has made possible the productionof laser pulses short enough to effectively probe the non-Markovian regime. Recent
experiments at JILA have made use of these advances in order to observe non-Markovian
dynamics in the collisions between atoms in a dense gas of potassium [90]. A common
signature of physics beyond the Markov approximation is revealed by lineshapes that

104
do not conform to a Lorentzian profile.
The field of nonlinear optical spectroscopy is well established and contains a refined set of concepts, formalism, and techniques [91]. The development in this chapter
aims to marry these established methods to a detailed microscopic theory of atomic collisions. Beginning in the 1950’s, physicists have approached the problem of a theoreticaltreatment of the phenomenon called pressure broadening, in which collisions within agas of atoms have the effect of broadening the spectral lines obtained by probing the gas
with light. An atom, by itself in free space, will exhibit a finite-width lineshape whenprobed. This lineshape is a result of the atom’s constant interaction with the surrounding photon vacuum, and it is this coupling that produces the phenomenon known as
spontaneous emission, a subject explored in detail elsewhere in this thesis. Introducing
the interaction with surrounding particles, which can be thought of as another reservoir
in addition to the photon vacuum, causes the lineshapes to be further broadened.
Some of the earliest work on the theory of pressure broadening was done by
Baranger [92, 93, 94], who developed an extensive theory utilizing relaxation theory,in which a single atom was modeled as interacting with a reservoir of other atoms, aswell as the photon reservoir. This provided the basis for further work by others such as
Fano [95, 96, 97] and Ben-Reuven [98, 99, 100, 101]. Fano’s biggest contribution to thisfield was to rework Baranger’s theory into a much more elegant formalism of relaxation
theory, as presented initially in Ref. [96] and in a more pedagogical format in Ref. [97].
Ben-Reuven, whose early work was with Fano, continued along these lines to extend
and apply the theory to pressure-broadening experiments.
The theory up to this point was mostly within what is called the impact approximation, in which it is assumed that particles interact with each other instantaneously
before moving on, and which is equivalent to the Markov approximation. Later work by
Cooper and Burnett [102, 103, 104], Silbey and Harris [105], Szudy and Balis [106, 107],
Allard and Keilkopf [108], and Mukamel [109] built on previous work, applying the105theory to newer experiments as well as considering physics outside the impact approximation.
In the theory presented in the chapter, the modern methods of nonlinear optical
spectroscopy will be used, with the physics of the interparticle interactions describedby extending the relaxation-theory techniques cited above. The nonlinear spectroscopy
techniques offer the advantage of being well suited to describing time-dependent phenomena and experimental techniques such as transient multiwave mixing. The established relaxation-theory techniques offer the advantage of a detailed description of particle scattering using realistic potentials and one or more collision channels. The main
goal of this chapter is to present the basic ideas of this theory, with the hope that they
can be used in the future to better understand the physics of the experiments currently
being conducted.

6.2

Basics of nonlinear optical spectroscopy
Many sources describe the formalism of nonlinear optics, and the basics will besummarized in this section before attempting to incorporate scattering relaxation theory.
In particular, Ref. [91] by Mukamel provides an elegant and detailed development of
the methods of nonlinear optics, and this section will employ the same notation andterminology. Beginning with the the Hamiltonian of the total system,(S)

H = H0(S)

where H0

(B)

+ H0

+ HSB + HSL (t),

(6.1)
(B)

is the Hamiltonian of the atom or molecule (system), H0

is the Hamilto-

nian of the bath, which could be the photon vacuum and/or other perturbing particles,
HSB is the interaction of the system with the bath, and HSL is the interaction of thesystem with all of the probe lasers. The density matrix obeys the Liouville equation,iidρ= − Lρ − Lint (t)ρ.
dt
~
~

(6.2)

106
Here, the following Liouville-space operators have been defined (see Appendix A for abrief review of Liouville space),
h
i
(S)
(B)
LA = H0 + H0 + VSB , A ,(6.3)

Lint A = [HSL (t), A] ,(6.4)

VA = [V, A] ,(6.5)

where V is the dipole operator such that HSL (t) = −E(r, t) · V . Expanding the density
matrix in powers of the electric field yields
ρ(t) = ρ(0) (t) + ρ(1) (t) + ρ(2) (t) + . . .

(6.6)

P (r, t) = Tr [V ρ(t)] = hhV |ρ(t)ii,(6.7)

The polarization,is expanded similarly, and the nth order contribution, P (n) (r, t) = hhV |ρ(n) (t)ii, can beshown to be
P

(n)

Z ∞
(r, t) =
0

Z ∞
dtn

0

Z ∞
dtn−1 . . .

0

dtn S (n) (tn , tn−1 , . . . , t1 )(6.8)

× E(r, t − tn )E(r, t − tn − tn−1 ) . . . E(r, t − tn − tn−1 . . . − t1 ),
with
S

(n)

µ ¶n
i
(tn , tn−1 , . . . , t1 ) =hhV |G(tn )VG(tn−1 )V . . . G(t1 )V|ρ(−∞)ii ,
~

(6.9)

the nth-order nonlinear response function. The retarded Liouville-space Green function
for evolution of the system and bath without interaction with the probe laser is given
by

µ
¶
i
G(t) = θ(t) exp − Lt .
~(6.10)

The first order response function gives the usual linear response of a system. The
present work will be primarily concerned with the first important nonlinear contribution.

107
Since the second-order contribution vanishes for homogeneous systems [91], this will be
the third-order nonlinear response function,
S (3) (t3 , t2 , t1 ) = hhV |G(t3 )VG(t2 )VG(t1 )V|ρ(−∞)ii .

(6.11)

This equation separates the physics due to interactions with the probe-laser field V,
from the physics of the atom interacting with the other atoms and with the vacuum
photon reservoir, described by the Green function G(t). This is convenient, and thetheory developed later derives an expression for this Green function utilizing a detailed
scattering theory framed in a relaxation-theory format. With a method of calculating
G(t) in hand, Eq. (6.11) can be used to determine a variety of features relevant totransient four-wave mixing experiments.

6.3

Phenomenological damping model
Some important concepts are now introduced, in the context of a simple model

damping. The effects of collisions on an atomic lineshape can be divided into two categories, homogeneous and inhomogeneous. Homogeneous broadening is a dynamical
effect that is caused by the repeated collisions an atom undergoes, leading to a fluctuation in the atom’s transition energy. This fluctuation, over time, results in a net
broadening. Inhomogeneous broadening, on the other hand, is a static effect due to thefact that atoms see different interaction-induced frequency shifts depending on their location, and becomes more prominent at low temperatures and in solids. In experiments,
the basic CW photon-echo technique allows for the cancelation of the effects of inhomogeneous broadening. Each atom’s transition frequency, having a static shift, causes thephase to evolve at a different rate when the system is probed. When the second pulse
is applied and the evolution reversed, the disparate atoms all catch up with each otherat the same time, and an echo pulse is emitted. Homogeneous broadening, however, isunaffected by this technique, and remains in the observed signal.

108
The JILA experiment mentioned above [90] uses transient four-wave mixing techniques to probe collision physics on very short time scales. Over such short probe
times, the collision physics that usually produces homogeneous broadening is now approximately static (the so-called quasi-static approximation), with each atom’s relative
distance changing very little over the probe time. Thus, although the scattering is stilla fundamentally dynamical process, the physics being probed is of the inhomogeneous
type. This means that a photon echo technique can be applied to cancel out this new
type of inhomogeneous broadening, and details about the interatomic potentials canby probed. As a result, this type of technique is generally not considered to be a true
photon echo, since there is no fundamentally static interaction shifts in the problem, as
there is in the traditional sense of photon echoes. Nevertheless, the basic idea is thesame.
A basic understanding of collisional damping can be obtained by following thephenomenological approach presented in many introductory quantum-optics textbooks,
a good example of which is Meystre and Sargent [57]. To begin with, take a simple
density-matrix model of a two level atom, ground state |gi and excited state |ei, including spontaneous emission with a rate Γ, i.e., the optical Bloch equations. Since
homogeneous collisional broadening is caused by a fluctuation of the atomic transition
frequency over time as the atom interacts with other atoms, it affects the off-diagonal
density-matrix element. Representing the fluctuation as δω(t) results in an equation of
motion for the coherence density-matrix term ofρ̇ge = − (iω0 + iδω(t) + Γ) ρge ,(6.12)

which has the formal solution
µ
¶
Z t
0
0
ρge (t) = ρge (0) exp −(iω0 + Γ)t −dt iδω(t ) .

(6.13)

0If, as is often the case, the time scale considered is long enough that the fluctuatingterm oscillates many times over that period, the integral will evaluate to zero. (Actually,

109
the integral might yield a nonzero constant, but this could just be absorbed into thebare transition frequency ω0 . This is the case for an atom interacting with the photon
vacuum: the resulting Lamb shift is considered to be a part of the bare atomic energy
splittings.) This will yield an equation describing exponential decay of the coherence
matrix element, and is the case for most descriptions of spontaneous emission.

6.4

Transient multiwave mixing using relaxation theory

6.4.1

Relaxation theory of static pressure broadening

It is now necessary to summarize the essential results of the relaxation theory ofpressure broadening. The work by Fano [96] formalizes and generalizes earlier work [92,93, 94], and the notation and techniques used by Fano will be employed here.
Consider a system, consisting of one, or possibly more atoms or molecules, and abath, consisting of a large number of perturbing atoms. The two can interact with eachother. Note that, for simplicity, spontaneous emission (another form of relaxation) has
been left out for now. However, this does not affect the equations in any nontrivial way,
and it can be easily inserted later. The total Hamiltonian is(S)

H = H0(S)

where H0

(B)

+ H0

+ VSB ,
(S)

is the Hamiltonian of the system, H0(6.14)

is the Hamiltonian of the particle

bath, and VSB is the interaction between the two.
An important part of the development in Ref. [96] is to show that the resolvent
operator for the total evolution can be disentangled into a contribution of the baresystem and bath, and a contribution of the interaction of the two by using a known
operator identity,
·
¸
1
1
1
=
1 + M (ω).
ω − L0 − L1
ω − L0ω − L0(6.15)The operator M (ω) contains all of the physics of the interaction, with everything else

110
simply free evolution, and it can be expressed in various forms,
1
L1
1 − L1 (ω − L0 )−11= L1
1 − (ω − L0 )−1 L1
1= L1 + L1
L1
ω − L0 − L1
¸n
∞ ·
X
1
= L1
L1 .
ω − L0M (ω) =

(6.16)

n=0

The interaction operator M is then expanded in terms of density,
M=

n
X
i

(1)
Mi +

n X
i−1
X
i=2

(2)
Mij +

j−1
n X
i−1 X
X

j

i

j

(3)

Mijk . . . ,(6.17)kwhere n is the number of atoms or molecules. Reference [96] restricts its analysis toinclude only the first term in the above expansion, describing only binary collisions, and(1)

defining m = Mi . The present work will make the same assumption.
Also assumed in Ref. [96], and common to most of the early work on pressure
broadening, is that the total density matrix can be factorized into a product of a systemdensity matrix and a bath density matrix,ρ = ρ(S) ρ(B) .

(6.18)The validity of this assumption is explored in later papers [100, 102] by other authors.
For present purposes, this factorization is assumed to be valid.
Having made the binary-collision approximation, the operator m, which it isuseful to think of as a Liouville-space version of the t-matrix, obeys the equationm(ω) = L + L

1t(ω),ω − L0(6.19)

which is a Liouville-space version of the Lippman-Schwinger equation. Ref. [96] showsthat this equation is solved by½
Z
11
ω − L0 ∞+iη∗dψt(ψ) +m(ω) =∗ t (ψ − ω)
2πi
ψ
−
H
ψ
−
ω
−
H
0
−∞+iη
0
(6.20)¾
1
1
ω − L0∗+t(ψ)t (ψ − ω),
ψ − H0 ψ − ω − H0∗(ψ − H0 ) (ψ − ω − H0∗ )

111
with
² = Imω > η > 0.

(6.21)

This ensures that Im(ψ − ω) < 0.

6.4.2

Extension to transient nonlinear-optics problems

Reference [96] manipulates the expression for the Liouville-space interaction operator given in Eq. (6.20), expanding it in a series of terms involving the two-particle
t-matrix. Although this series of terms is rather complicated, it provides a starting
point for making approximations.
Now, a method will be developed that utilizes m, obtained in some manner, andapplies it to the nonlinear optics formalism earlier in this chapter. This means obtaining
a Green function that can be used in determining the third-order nonlinear response
function in Eq. (6.11). Using Eq. (6.15) and taking the Fourier transform to go to thetime domain, the time-domain Green function for evolution of the system interactingwith the perturber bath is given by
1
G(t) = G0 (t) −
2πi

Z ∞
−∞

dωG0 (ω)m(ω)G0 (ω)e−iωt .

(6.22)

Expanding the frequency-space Green function in a basis of bare-system Liouville-space
eigenstates,1ω − L0 + i²
X |nmiihhnm|,=
ω − ωnm + i²n,mG0 (ω) =(6.23)

112
which yields
G(t) = G0 (t) −X

1|nmiihhn m |
2πi
0
0

0

n,m,n0 ,mX= G0 (t) +

Z ∞
dω
−∞

hhnm|m(ω)|n0 m0 iie−iωt(ω − ωnm + i²)(ω − ωn0 m0 + i²)|nmiihhn0 m0 | θ(t)n,m,n0 ,m0 ,ωnm 6=ωn0 m0µ

×mnm,n0 m0 (ωnm ) −iωnm t mnm,n0 m0 (ωn0 m0 ) −iωn0 m0 te+e
ωnm − ωn0 m0
ωn0 m0 − ωnm¶(6.24)
In obtaining the above expression, integration has been performed along the real line,
skirting the poles and obtaining residues.
If the methods in Ref. [96] are used to calculate m, then it can be applied toEq. (6.24) to determine a time-dependent Green function for use in Eq. (6.11). In
turn, this quantity will provide a powerful means of understanding the physics of nonMarkovian dynamics in experiments.

6.5

Conclusion
The theory in this chapter represents a preliminary idea for merging the elegant

techniques of relaxation theory — applied successfully to static pressure broadening in
the past — to the time-dependent nonlinear optics experiments being conducted today.
In it’s current form, it is an interesting exercise in the methods of Liouville space, systemreservoir interactions, and relaxation theory that have been discussed extensively in thisthesis. Additionally, the ideas here show promise for obtaining a better understandingof non-Markovian dynamics.

Bibliography[1] J. W. Dunn, D. Blume, B. Borca, B. E. Granger, and C. H. Greene. Feshbachresonance cooling of trapped atom pairs. Phys. Rev. A, 71(3):033402, 2005.
[2] J. W. Dunn and C. H. Greene. Predictions of laser-cooling temperatures for
multilevel atoms in three-dimensional polarization-gradient fields. Phys. Rev. A,
73(3):033421, 2006.
[3] J. W. Dunn, J. W. Thomsen, C. H. Greene, and F. C. Cruz. Coherent Quantum
Engineering of Free-Space Laser Cooling. submitted to Phys. Rev. Lett., Arxiv
preprint physics/0610272, 2006.
[4] A. Einstein. Uber die von der molekularkinetischen Theorie der Warme geforderte
Bewegung von in ruhenden Flussigkeiten suspendierten Teilchen. Ann. Phys.
(Leipzig), 17:549.
[5] M. Smoluchowski. Zur kinetischen Theorie der Brownschen Molekularbewegung
und der Suspensionen. Ann. Phys. (Leipzig), 21(756), 1906.
[6] C. W. Gardiner. Handbook of Stochastic Methods. Springer, New York, 1985.
[7] P. Langevin. Sur la theorie du mouvement brownien. C. R. Acad. Sci., 46:530,
1908.
[8] C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg.
Interactions. Wiley, New York, 1998.

Atom-Photon[9] K. Berg-Sorenson, Y. Castin, E. Bonderup, and K. Molmer. Momentum diffusionof atoms moving in laser fields. J. Phys. B, 25(20):4195, 1992.
[10] K. Mølmer, Y. Castin, and J. Dalibard. Monte Carlo wave-function method inquantum optics. J. Opt. Soc. Am. B, 10(3):524, 1993.
[11] Y. Castin and K. Mølmer. Monte Carlo Wave-Function Analysis of 3D Optical
Molasses. Phys. Rev. Lett., 74(19):3772, 1995.
[12] R. Dum, P. Zoller, and H. Ritsch. Monte Carlo simulation of the atomic masterequation for spontaneous emission. Phys. Rev. A, 45(7):4879, 1992.
[13] N. Gisin and I. C. Percival. Wave-function approach to dissipative processes: arethere quantum jumps? Phys. Lett. A, 167(4):315, 1992.

114[14] N. Gisin and I. C. Percival. The quantum-state diffusion model applied to open
systems. J. Phys. A, 25(21):5677, 1992.
[15] H. J. Carmichael. An Open-Systems Approach to Quantum Optics. SpringerVerlag, Berlin, 1993.
[16] J. Dalibard and C. Cohen-Tannoudji. Laser cooling below the Doppler limit bypolarization gradients: simple theoretical models. J. Opt. Soc. Am. B, 6(11):2023,
1989.
[17] P. J. Ungar, D. S. Weiss, E. Riis, and S. Chu. Optical molasses and multilevelatoms: theory. J. Opt. Soc. Am. B, 6(11):2058, 1989.
[18] P. D. Lett, W. D. Phillips, S. L. Rolston, C. E. Tanner, R. N. Watts, and C. I.
Westbrook. Optical molasses. J. Opt. Soc. Am. B, 6(11):2084, 1989.
[19] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell.
Observation of Bose-Einstein Condensation in a Dilute Atomic Vapor. Science,
269(5221):198, 1995.
[20] K. B. Davis, M. O. Mewes, and W. Ketterle. An analytical model for evaporative
cooling of atoms. Appl. Phys. B, 60(2):155, 1995.
[21] J. Herbig, T. Kraemer, M. Mark, T. Weber, C. Chin, H. C. Nägerl, and R. Grimm.
Preparation of a Pure Molecular Quantum Gas. Science, 301(5639):1510, 2003.
[22] S. Dürr, T. Volz, A. Marte, and G. Rempe. Observation of Molecules Produced
from a Bose-Einstein Condensate. Phys. Rev. Lett., 92(2):020406, 2004.
[23] K. E. Strecker, G. B. Partridge, and R. G. Hulet. Conversion of an Atomic Fermi
Gas to a Long-Lived Molecular Bose Gas. Phys. Rev. Lett., 91(8):080406, 2003.
[24] K. Xu, T. Mukaiyama, J. R. Abo-Shaeer, J. K. Chin, D. E. Miller, and W. Ketterle. Formation of Quantum-Degenerate Sodium Molecules. Phys. Rev. Lett.,
91(21):210402, 2003.
[25] H. Feshbach. Theoretical Nuclear Physics. Wiley, New York, NY, 1992.
[26] W. C. Stwalley. Stability of Spin-Aligned Hydrogen at Low Temperatures andHigh Magnetic Fields: New Field-Dependent Scattering Resonances and Predissociations. Phys. Rev. Lett., 37(24):1628, 1976.
[27] E. Tiesinga, B. J. Verhaar, and H. T. C. Stoof. Threshold and resonance phenomena in ultracold ground-state collisions. Phys. Rev. A, 47(5):4114, 1993.
[28] B. Borca, D. Blume, and C. H. Greene. A two-atom picture of coherent atommolecule quantum beats. New J. Phys., 5(1):111, 2003.
[29] T. Busch, B. G. Englert, K. Rzażewski, and M. Wilkens. Two Cold Atoms in a
Harmonic Trap. Found. Phys., 28(4):549, 1998.
[30] D. Blume and C. H. Greene. Fermi pseudopotential approximation: Two particles
under external confinement. Phys. Rev. A, 65(4):043613, 2002.

115[31] B. Borca. Two-Channel Zero-Range-Potential Model for the Feshbach Resonance.
Private communication, 2003.
[32] C. H. Greene. Two-body Ultracold Feshbach Resonance Behavior versus Energy
and Magnetic Field. Private communication, 2004.
[33] M. Greiner, C. A Regal, and D. S. Jin. Probing the Excitation Spectrum of a
Fermi Gas in the BCS-BEC Crossover Regime. Phys. Rev. Lett., 94(7):070403,
2005.
[34] U. Fano. Effects of Configuration Interaction on Intensities and Phase Shifts.
Phys. Rev., 124(6):1866, 1961.
[35] E. L. Bolda, E. Tiesinga, and P. S. Julienne. Effective-scattering-length modelof ultracold atomic collisions and Feshbach resonances in tight harmonic traps.
Phys. Rev. A, 66(1):013403, 2002.
[36] F. H. Mies, E. Tiesinga, and P. S. Julienne. Manipulation of Feshbach resonances
in ultracold atomic collisions using time-dependent magnetic fields. Phys. Rev.
A, 61(2):022721, 2000.
[37] W. Ketterle and D. E. Pritchard. Atom cooling by time-dependent potentials.
Phys. Rev. A, 46(7):4051, 1992.
[38] D. Frese, B. Ueberholz, S. Kuhr, W. Alt, D. Schrader, V. Gomer, and D. Meschede.
Single Atoms in an Optical Dipole Trap: Towards a Deterministic Source of Cold
Atoms. Phys. Rev. Lett., 85(18):3777, 2000.
[39] N. Schlosser, G. Reymond, I. Protsenko, and P. Grangier. Sub-poissonian loadingof single atoms in a microscopic dipole trap. Nature (London), 411(6841):1024,
2001.
[40] M. D. Barrett, J. A. Sauer, and M. S. Chapman. All-Optical Formation of an
Atomic Bose-Einstein Condensate. Phys. Rev. Lett., 87(1):010404, 2001.
[41] P. Rabl, A. J. Daley, P. O. Fedichev, J. I. Cirac, and P. Zoller. Defect-Suppressed
Atomic Crystals in an Optical Lattice. Phys. Rev. Lett., 91(11):110403, 2003.
[42] P. D. Lett, R. N. Watts, C. I. Westbrook, W. D. Phillips, P. L. Gould, and H. J.
Metcalf. Observation of Atoms Laser Cooled below the Doppler Limit. Phys. Rev.
Lett., 61(2):169, 1988.
[43] A. Bambini and A. Agresti. Radiative cooling force in atoms with multiplet
structure. Phys. Rev. A, 56(4):3040, 1997.
[44] K. Mølmer. Friction and diffusion coefficients for cooling of atoms in laser fields
with multidimensional periodicity. Phys. Rev. A, 44(9):5820, 1991.
[45] J. Javanainen. Numerical experiments in semiclassical laser-cooling theory ofmultistate atoms. Phys. Rev. A, 46(9):5819, 1992.
[46] J. Javanainen. Polarization gradient cooling in three dimensions: comparison oftheory and experiment. J. Phys. B, 27(3):L41, 1994.

116[47] A. P. Kazantsev and I. V. Krasnov. Rectification of the gradient force of resonantradiation pressure. JETP Lett., 46:332, 1987.
[48] R. Grimm, Y. B. Ovchinnikov, A. I. Sidorov, and V. S. Letokhov. Observation ofa strong rectified dipole force in a bichromatic standing light wave. Phys. Rev.
Lett., 65(12):1415, 1990.
[49] P. R. Hemmer, M. S. Shahriar, M. Prentiss, D. P. Katz, K. Berggren, J. Mervis,
and N. P. Bigelow. First Observation of Forces on Three-Level Atoms in Raman
Resonant Standing-Wave Optical Fields. Phys. Rev. Lett., 68(21):3148, 1992.
[50] R. Gupta, C. Xie, S. Padua, H. Batelaan, and H. Metcalf. Bichromatic Laser
Cooling in a Three-Level System. Phys. Rev. Lett., 71(19):3087, 1993.
[51] M. Drewsen. Effective stimulated laser cooling of atoms with a three-level λconfiguration by two negatively detuned standing light waves. Phys. Rev. A,
51(2):1407, 1995.
[52] X. Xu, T. H. Loftus, J. W. Dunn, C. H. Greene, J. L. Hall, A. Gallagher, andJ. Ye. Single-Stage Sub-Doppler Cooling of Alkaline Earth Atoms. Phys. Rev.
Lett., 90(19):193002, 2003.
[53] C. Fort, A. Bambini, L. Cacciapuoti, FS Cataliotti, M. Prevedelli, GM Tino, and
M. Inguscio. Cooling mechanisms in potassium magneto-optical traps. Euro.
Phys. J. D, 3(2):113, 1998.
[54] U. Schünemann, H. Engler, M. Zielonkowski, M. Weidemüller, and R. Grimm.
OpticsMagneto-optic trapping of lithium using semiconductor lasers.
Communications, 158(1):263, 1998.
[55] R. Maruyama, R. H. Wynar, M. V. Romalis, A. Andalkar, M. D. Swallows, C. E.
Pearson, and E. N. Fortson. Investigation of sub-Doppler cooling in an ytterbium
magneto-optical trap. Phys. Rev. A, 68(1):011403, 2003.
[56] D. Lucas, P. Horak, and G. Grynberg. Sisyphus cooling of rubidium atoms on theD2 (F = 1 − F 0 = 1) line: The role of the neighbouring transitions. Euro. Phys.
J. D, 7(3), 1999.
[57] P. Meystre and M. Sargent III. Elements of Quantum Optics. Springer-Verlag,
Berlin, 1998.
[58] C. W. Gardiner and P. Zoller. Quantum Noise. Springer-Verlag, Berlin, 2004.
[59] H. Spohn. Kinetic equations from Hamiltonian dynamics: Markovian limits. Rev.
Mod. Phys., 52(3):569, 1980.
[60] C. Cohen-Tannoudji. In J. Dalibard, J.-M. Raimond, and J. Zinn-Justin, editors,
Fundamental Systems in Quantum Optics, pages 1–164, Amsterdam, 1992. NorthHolland.
[61] Y. Castin, H. Wallis, and J. Dalibard. Limit of Doppler cooling. J. Opt. Soc. Am.
B, 6(2046):128, 1989.

117[62] H. J. Kluge and H. Sauter. Levelcrossing experiments in the first excited 1 P1states of the alkaline earths. Z. Phys., 270(4):295, 1974.
[63] K. Blum. Density Matrix Theory and Applications. Plenum, New York, NY,
1981.
[64] H. Katori, T. Ido, Y. Isoya, and M. Kuwata-Gonokami. Magneto-Optical Trappingand Cooling of Strontium Atoms down to the Photon Recoil Temperature. Phys.
Rev. Lett., 82(6):1116, 1999.
[65] T. Kuwamoto, K. Honda, Y. Takahashi, and T. Yabuzaki. Magneto-optical trapping of Yb atoms using an intercombination transition. Phys. Rev. A, 60(2):R745,
1999.
[66] P. Gill, editor. Proceedings of the 6th Symposium on Frequency Standards andMetrology, St. Andrews, Scotland. World Scientific, 2001.
[67] M. Takamoto, F. L. Hong, R. Higashi, and H. Katori. An optical lattice clock.
Nature (London), 435(7040):321, 2005.
[68] H. J. Metcalf and P. van der Straten. Laser Cooling and Trapping. SpringerVerlag, New York, 1999.
[69] M. Kasevich and S. Chu. Laser Cooling Below a Photon Recoil with Three-Level
Atoms. Phys. Rev. Lett., 69(12):1741, 1992.
[70] A. Aspect, E. Arimondo, R. Kaiser, N. Vansteenkiste, and C. Cohen-Tannoudji.
Laser cooling below the one-photon recoil energy by velocity-selective coherent
population trapping: theoretical analysis. J. Opt. Soc. Am. B, 6(11):2112, 1989.
[71] H. W. Chan, A. T. Black, and V. Vuletić. Observation of Collective-EmissionInduced Cooling of Atoms in an Optical Cavity. Phys. Rev. Lett., 90(6):063003,
2003.
[72] C. A. Sackett, J. Chen, J. J. Tollett, and R. G. Hulet. Dipole-Force Laser Coolingin Multilevel Atoms. Laser Physics, 4:861, 1994.
[73] L. Tan, C. K. Zhong, C. L. Zhang, and Q. Zhang. Two-Photon Doppler Cooling of
Rubidium Atoms to the Recoil Limit Using Cascade Transitions. Chin. J. Phys.,
44(3), 2006.
[74] G. Morigi and E. Arimondo. Two-photon and EIT-assisted Doppler cooling in athree-level cascade system. Arxiv preprint quant-ph/0610223, 2006.
[75] G. Morigi, J. Eschner, and C. H. Keitel. Ground State Laser Cooling Using
Electromagnetically Induced Transparency. Phys. Rev. Lett., 85(21):4458, 2000.
[76] C. F. Roos, D. Leibfried, A. Mundt, F. Schmidt-Kaler, J. Eschner, and R. Blatt.
Experimental Demonstration of Ground State Laser Cooling with Electromagnetically Induced Transparency. Phys. Rev. Lett., 85(26):5547, 2000.

118[77] I. Marzoli, J. I. Cirac, R. Blatt, and P. Zoller. Laser cooling of trapped three-levelions: Designing two-level systems for sideband cooling. Phys. Rev. A, 49(4):2771,
1994.
[78] M. Fleischhauer, A. Imamoglu, and J. P. Marangos. Electromagnetically induced
transparency: Optics in coherent media. Rev. Mod. Phys., 77(2):633, 2005.
[79] W. C. Magno, R. L. Cavasso Filho, and F. C. Cruz. Two-photon Doppler coolingof alkaline-earth-metal and ytterbium atoms. Phys. Rev. A, 67(4):043407, 2003.
[80] N. Malossi, S. Damkjær, P. L. Hansen, L. B. Jacobsen, L. Kindt, S. Sauge, J. W.
Thomsen, F. C. Cruz, M. Allegrini, and E. Arimondo. Two-photon cooling ofmagnesium atoms. Phys. Rev. A, 72(5):051403, 2005.
[81] R. CiuryÃlo, E. Tiesinga, and P. S. Julienne. Optical tuning of the scattering lengthof cold alkaline-earth-metal atoms. Phys. Rev. A, 71(3):030701, 2005.
[82] Y. Takasu, K. Maki, K. Komori, T. Takano, K. Honda, M. Kumakura,
T. Yabuzaki, and Y. Takahashi. Spin-Singlet Bose-Einstein Condensation of TwoElectron Atoms. Phys. Rev. Lett., 91(4):040404, 2003.
[83] V. Vuletić and S. Chu. Laser Cooling of Atoms, Ions, or Molecules by Coherent
Scattering. Phys. Rev. Lett., 84(17):3787, 2000.
[84] V. Vuletić, H. W. Chan, and A. T. Black. Three-dimensional cavity Doppler
cooling and cavity sideband cooling by coherent scattering. Phys. Rev. A,
64(3):033405, 2001.
[85] A.T. Black, H.W. Chan, and V. Vuletić. Observation of Collective Friction Forces
due to Spatial Self-Organization of Atoms: From Rayleigh to Bragg Scattering.
Phys. Rev. Lett., 91(20):203001, 2003.
[86] P. Domokos and H. Ritsch. Collective Cooling and Self-Organization of Atoms ina Cavity. Phys. Rev. Lett., 89(25):253003, 2002.
[87] S. Zippilli, G. Morigi, and H. Ritsch. Suppression of Bragg Scattering by Collective
Interference of Spatially Ordered Atoms with a High-Q Cavity Mode. Phys. Rev.
Lett., 93(12):123002, 2004.
[88] P. Domokos and H. Ritsch. Mechanical effects of light in optical resonators. J.
Opt. Soc. of Am. B, 20(5):1098, 2003.
[89] H. J. Carmichael. Statistical Methods in Quantum Optics: Master Equations andFokker-Planck Equations. Springer, 1999.
[90] V. O. Lorenz and S. T. Cundiff. Non-markovian dynamics in a dense potassium
vapor. Phys. Rev. Lett., 95(16):163601, 2005.
[91] S. Mukamel. Principles of Nonlinear Optical Spectroscopy. Oxford University
Press, USA.
[92] M. Baranger. Simplified Quantum-Mechanical Theory of Pressure Broadening.
Phys. Rev., 111(2):481, 1958.

119[93] M. Baranger. Problem of Overlapping Lines in the Theory of Pressure Broadening.
Phys. Rev., 111(2):494, 1958.
[94] M. Baranger. General Impact Theory of Pressure Broadening.
112(3):855, 1958.

Phys. Rev.,[95] U. Fano. Description of States in Quantum Mechanics by Density Matrix andOperator Techniques. Rev. Mod. Phys., 29(1):74, 1957.
[96] U. Fano. Pressure Broadening as a Prototype of Relaxation.
131(1):259, 1963.

Phys. Rev.,[97] U. Fano. Lectures on the many-body problem, vol. 2, 1964.
[98] A. Ben-Reuven. Impact Broadening of Microwave Spectra. Phys. Rev., 145(1):7,
1966.
[99] A. Ben-Reuven. Symmetry Considerations in Pressure-Broadening Theory. Phys.
Rev., 141(1):34, 1966.
[100] A. Ben-Reuven. Resonance broadening of spectral lines. Phys. Rev. A, 4(6):2115,
1971.
[101] Y. Prior and A. Ben-Reuven. Nonimpact theory of four-wave mixing and intracollisional dynamics. Phys. Rev. A, 33(4):2362, 1986.
[102] K. Burnett, J. Cooper, R. J. Ballagh, and E. W. Smith. Collisional redistributionof radiation. I. The density matrix. Phys. Rev. A, 22(5):2005, 1980.
[103] K. Burnett and J. Cooper. Collisional redistribution of radiation. II. The effects
of degeneracy on the equations of motion for the density matrix. Phys. Rev. A,
22(5):2027, 1980.
[104] K. Burnett and J. Cooper. Collisional redistribution of radiation. III. The equationof motion for the correlation function and the scattered spectrum. Phys. Rev. A,
22(5):2044, 1980.
[105] R. Sibley and R. A. Harris. Tunneling of Molecules in Low-Temperature Media:
An Elementary Description. J. Phys. Chem., 93:7062, 1989.
[106] J. Szudy and W. E. Baylis. Unified Franck-Condon Treatment of Pressure Broadening of Spectral Lines. J. Quant. Spectrosc. Radiat. Transfer, 15:641, 1975.
[107] J. Szudy and W. E. Baylis. Profiles of Line Wings and Rainbow Satellites Associated with Optical and Radiative Collision. Phys. Rep., 266(3-4):127, 1996.
[108] N. Allard and J. Kielkopf. The effect of neutral nonresonant collisions on atomic
spectral lines. Rev. Mod. Phys., 54(4):1103, 1982.
[109] S. Mukamel. Nonimpact unified theory of four-wave mixing and two-photon processes. Phys. Rev. A, 28(6):3480, 1983.

120[110] C. R. Vidal, J. Cooper, and E. W. Smith. Hydrogen Stark Broadening Calculations
with the Unified Classical Path Theory, volume 116. National Bureau of Standards, Boulder, CO, 1970.
[111] C. Cohen-Tannoudji. In R. Balian, S. Haroche, and S. Liberman, editors, Frontiersin Laser Spectroscopy, volume 1, pages 3–104, Amsterdam, 1977. North-Holland.
[112] I. I. Sobelman. Atomic Spectra and Radiative Transitions. Springer-Verlag, Berlin,
1992.

Appendix A
Liouville Space

The primary purpose of the Liouville-space formalism is to provide a convenient
notation for expressing operations of the type often encountered in relaxation problems
and using density matrices. According to Ref. [91], the term “Liouville space” was firstused by Fano in his 1963 publication on pressure broadening [96], and his work throughthe years made extensive use of these techniques. Fano’s paper, in fact, is a very good
illustration of the utility of the Liouville-space formalism. In it, equations are developed
to describe the interaction of a single particle with a surrounding gas of particles, whichis represented as a reservoir, and thus requires the use of a density matrix for theatom. By using the Liouville-space formalism, Fano is able to derive an equation forthe operator m, given in Eq. (6.19) from Chapter 6 of this thesis, which contains allof the physics of the system-reservoir interaction. This equation has the form of the
Lippman-Schwinger equation obtained in two-particle scattering theory, and, indeed,just as m is the Liouville-space version of the t-matrix, Eq. (6.19) is the Liouville-spaceversion of the Lippman-Schwinger equation. Analogous equations such as this illustrate
the power that using such a notation can offer.
The Liouville space formalism is described in a variety of places in the literature [91, 63, 96, 97, 102, 110], and the intention of this appendix is to provide a basic
motivation for the formalism, and summarize its main tenants, for reference throughout the thesis. For a more detailed account, the reader is referred to any one of the122
aforementioned references.
Ordinary quantum mechanics, using wave functions, is formulated in Hilbert
space. The density matrix, a generalization of the wave function to account for non-purestates, can also be used in Hilbert space. This is, in fact, how the Liouville equation as
it is usually written,iρ̇(t) = − [H, ρ(t)] ,
~

(A.1)

for the density operator ρ and the Hamiltonian H. This equation cannot be expressed asa linear system, since the commutator produces terms in which ρ appears left multiplying
the Hamiltonian. In Liouville-space notation, an operator called a Liouvillian or asuperoperator is introduced. This operator L, when acting upon the density operator,
provides the required commutator, and thus the Liouville equation written in Liouvillespace notation isiρ̇(t) = − Lρ.
~

(A.2)

The Liouville equation is now written in a form that is equivalent to the Schrödinger
equation. This equivalence turns out to be very useful because many operations performed in Liouville space will appear to be formally identical to familiar operations
performed in Hilbert space. Equation (A.2), when expanded in a basis, denotes the
N × N density matrix now as a vector of length N 2 , and the superoperator L appears
as a N 2 × N 2 matrix. This is now a linear system that can, for example, by plugged in
to a linear solver on a computer. However, this advantage has come at a cost, since thesize of the matrices involved are much larger than the Hamiltonian matrix represented
in Hilbert space.
Superoperators have four indices when referring back to the original Hilbert-space
operations, two that denote left multiplication of the density operator, and two that
denote right multiplication of the density operator. Explicitly, the operator L given

123
above is expressed in terms of the Hamiltonian as∗
Lij,kl = Hik δjl − Hjl
δik .

(A.3)

This four-index notation is often referred to as tetradic notation [91, 110].
Abstract quantum-state vectors in Liouville-space notation are often, especially inthe more modern literature, written with two angle brackets as |aii. The completeness
relation in Liouville space isX|ijiihhij| = 1,

(A.4)

i,j

and vectors can be expanded in basis in a manner analogous to Hilbert space,|aii =

X|ijiihhij| |aii.

(A.5)

i,j

The formalism of Liouville space is used extensively in Chapter 6 of this thesis.
There, the formalism allows a great deal of simplification of the resulting equations,and makes clear the analogy to similar but less general equations in Hilbert space. The
methods of Green functions is utilized in Liouville-space form. Liouville space concepts
are also very useful when performing numerical calculations, and have been used inevery calculation involving density matrices in this thesis. The author has constructeda suite of Fortran subroutines that allow conversion back and forth between Hilbert
space and Liouville space.

Appendix B
Spontaneous-Emission Relaxation Operator for Atoms with Complex
Internal StructureThis appendix is concerned with the derivation of the spontaneous-emission relaxation operator for an atom with hyperfine degeneracy, taking into account multiple
hyperfine excited-state manifolds. This type of operator is important for atoms having
multiple excited-state manifolds which are spaced in energy on the order of or smaller
than the natural linewidth, the topic addressed in Chapter 3. When these atoms are
driven with a laser near the transition frequency — as in, for example, laser cooling —coupling occurs for each of the excited-state manifolds and their effects must be considered. To understand the dynamics of such an atom in a polarization-gradient laser
field, it is thus important to accurately describe the coherent effects of each substate ineach hyperfine manifold, interacting with the vacuum photon field.
In addition to providing a detailed derivation of a quantity relevant to our numerical calculations, this appendix is meant to illustrate in a general manner the development of a relaxation operator. Through this rather elaborate example, it is shownhow the description of a small system, interacting with a large reservoir approximated(or exactly described) as made up of a large number of harmonic oscillators, can be
reduced to a simple, intuitive, parameterized description, as long as certain criteria(such as the Markov approximation) are met. In this manner, spontaneous emission, aphenomenon experimentally observed to behave according to a single parameter Γ, is125
built up from a microscopic description of the atom interacting with the vacuum photon field. From this description, the phenomenologically determined Γ can be given afundamental expression in terms of a series of microscopic interactions.
To begin, this problem is mapped onto the general class of system-reservoir interactions. In the present case, the atom is the system and the vacuum photon field isthe reservoir. The full Hamiltonian for the system plus the reservoir is
H = HA + HR + VA−R ,(B.1)

where HA is the bare atomic Hamiltonian, which includes the internal states of the atom
as well as its center of mass, HR is the Hamiltonian of the radiation field, and VA−R isthe interaction between the two. The interaction is described within the electric-dipole
approximation as
VA−R = −D · E(R)X=−(−1)q Dq E−q (R)q=0,±1= −iX(−1)q Dq

q=0,±1X X
i²(i) ⊥ki

(B.2)
r

i
~ωi h (i) iki ·R(i) † −iki ·R
²
a
e
−
²
a
e
,
i
q
q
i
2²0 L3where D is the electric dipole operator for the atom, and E is the electric-field operator
for the photon field. On the second line, the interaction has been expanded into its
spherical components and on the third line the radiation field has been written in an
explicit quantum-mechanical form, quantized in a box of length L, in terms of raisingand lowering operators, a† and a, respectively. The sum over i indicates a sum over fieldmodes, and for each ki there are two mutually perpendicular polarization vectors ²̂(i) .
To simplify the formalism, the atomic center of mass will initially be ignored, setting
the position coordinate to be the origin, R = 0. This will not affect the derivation
fundamentally, and in the end the equations will be easily generalized to include thecenter-of-mass degrees of freedom.
The full system-plus-reservoir density matrix ρ evolves in time according to the126
Liouville equation,ρ̇(t) =i[ρ(t), H] .
~

(B.3)

This equation describes the evolution of all degrees of freedom of the problem, including
the vacuum photon field. This is unnecessarily detailed and will be impractical to solve
in most situations. The goal is to find an equation of motion for the reduced systemdensity matrix, defined by taking the trace over all reservoir degrees of freedom of thefull density matrix,
σ = TrR ρ.

(B.4)

To this end, it is useful to utilize some of the the well-established methods of quantum
system-reservoir theory. It is impractical to present a detailed discussion of this methodology here, but these techniques can be easily found in quantum-optics literature, anda particularly thorough presentation is given by Cohen-Tannoudji in Refs. [8] and [111].
Summarizing, it is assumed that the reservoir is unaffected by the interactions with thesystem and that the reservoir has a short memory (the Markov approximation). With
these assumptions, an equation of motion is obtained for the system reduced density
operator σ with temporal coarse graining on a time scale much larger than the timescale of reservoir fluctuations but shorter than the time scale of evolution of the system,i1σ̇(t) = [σ(t), HA ] − 2
~
~

Z ∞
dτ
0

X
(−1)q
q

n
h
i
o
× gq (τ ) Dq e−iHA τ /~ D−q eiHA τ /~ σ(t) − e−iHA τ /~ Dq eiHA τ /~ σ(t) D−q + H.c. .
(B.5)
In the previous equation, gq (τ ) is the two-time correlation function for vacuum fluctuations of the photon reservoir and is given byZ~cω 3 −iωτgq (τ ) =edω
3²0 (2π)3c4hi= TrR σR Ẽq (τ )Ẽq (0) .

(B.6)

127
where c is the speed of light and ²0 is the permittivity of free space, and where thevariables with tildes are operators in the interaction representation,Ẽq (t) = eiHR t/~ Eq e−iHR t/~ .

(B.7)

It is assumed that the reservoir is initially the vacuum state, so that σR = |0i h0|. From
this it can be seen thatgq (τ ) =X|hν | Eq | 0i|2 e−iων t ,

(B.8)

ν

where the kets and bras refer to reservoir states. Note that g(τ )∗ = g(−τ ). The
correlation time of the reservoir τC is defined such that g(τ ) → 0 for τ À τC .
In addition to the above approximations, the secular approximation will also be
made, which requires that the equation of motion for each density-matrix element σ̇ij
have only terms involving density-matrix elements σkl on the right-hand side such that
|ωij − ωkl | ¿ γ, where γ is the order of magnitude of the system-reservoir coupling. The
following discussion focuses on a system with a ground state coupled to multiple excitedstates that are separated in energy of the order or smaller than γ. Thus, the groundexcited energy splitting |ωge | À γ will be a non-secular frequency, while ωei ej ∼ γ will
be a secular frequency.
The particular atomic system under consideration consists of an ground statewith electronic angular momentum J = Jg = 0 and an excited state with J = Je = 1.
The electronic angular momentum is coupled to the nuclear spin quantum number I,
resulting in a ground state with total angular momentum Fg = I, and three excitedstates with {Fei } = {I − 1, I, I + 1}. Note, however, that this derivation can be easily
extended to arbitrary angular momentum schemes. It is useful to decompose the systemdensity operator as illustrated in Fig B.1,
σ̇(t) = σ̇gg (t) +Xi,jσ̇ei ej (t) +X[σ̇ei g (t) + σ̇gei (t)] ,

(B.9)

i

where σij (t) = Pi σ(t)Pj ; Pi is a projection operator onto the i-th hyperfine manifold,P
Pi = Mi |JIF Mi i hJIF Mi |; and Mi is the substate label for the i-th manifold.

128

gg

ge1 ge2ge3

e1g e1e1 e1e2

e1e3

e2g e2e1 e2e2

e2e3

e3g e3e1 e3e2

e3e3

Figure B.1: The partitioning of the density operator for an atom with multiple coupled
excited-state manifolds, each potentially having multiple substates.

129
It will be useful to first work out some relations involving sums over electric-dipole
matrix elements related to reduced matrix elements and raising and lowering operators:
X

hMg |Dq |Mei i |Mg i hMei |Mg ,MeiX=

hJg IFg Mg |Dq |Je IFei Mei i |Jg IFg Mg i hJe IFei Mei |Mg ,MeiX=

Fg +Mg +Jg +I+Fei +1(−1)

Mg ,Mei
1 Fei   Jgq Mei Fei



 Fg
×
−Mg

q
(2Fg + 1)(2Fei + 1)

Fg



I

Je


1

(B.10)

hJg ||D||Je i |Jg IFg Mg i hJe IFei Mei |hJg ||D||Je i
=A(i)
,
q √
2Je + 1
andX(−1)q hMei |D−q |Mg i |Mei i hMg |Mg ,Mei=X(−1)q hJe IFei Mei |D−q |Jg IFg Mg i |Je IFei Mei i hJg IFg Mg |Mg ,Mei=

X(−1)Fei +Mei +Je +I+Fg +1+q

Mg ,Mei



 Fei
×
−Mei

q
(2Fg + 1)(2Fei + 1)
1 Fg   Je Fei−q Mg Fg Jg



I

1

hJe ||D||Jg i |Je IFei Mei i hJg IFg Mg |† hJe ||D||Jg i√= − A(i),q
2Je + 1

(B.11)(i) †where Aq(i)and Aq are the atomic raising and lowering operators defined in Eq. (3.6)

and the symmetry properties of the three-J and six-J symbols have been exploited (seeRef. [112]).
The equation of motion for each type of sub-density-matrix in Eq. (B.9) will now
be examined. Consider first the equation for the ground-state sub-density-matrix σgg (t),130
which evolves according to
Z ∞ X
1σ̇gg (t) = − 2
dτ
(−1)q
~ 0
q




X
X
(i)
(i)
× g(τ ) Pg Dq Pei D−q Pg σ(t)Pg e−iω0 τ −Pg Dq Pei σ(t)Pej D−q Pg eiω0 τ 

i
i,j



X
X
(j)
(i)
−iω0 τ iω0 τPg Dq Pei σ(t)Pej D−q Pg e,
Pg σ(t)Pg Dq Pei D−q Pg e−+g(−τ )
i,ji

(B.12)(i)where ω0 = ωei − ωg . It will first be shown that the first and third terms combine togive zero. This first term can be rearranged as
−

1
~2

Z ∞
dτ
0

X
X
(i)
(−1)q g(τ )Pg Dq Pei D−q Pg σ(t)Pg e−iω0 τ
q

i

Z ∞
(i)
1 XX
=− 2
dτ g(τ )e−iω0 τ , (B.13)(−1)q Pg Dq Pei D−q Pg σ(t)Pg~ q0i

and the third term can be rearranged as
1
− 2
~

Z ∞
dτ
0

X
X
(i)
(−1)q g(−τ )Pg σ(t)Pg Dq Pei D−q Pg eiω0 τ
q

i

Z ∞
(i)
1 XX
q
=− 2
(−1) Pg σ(t)Pg Dq Pei D−q Pgdτ g(−τ )eiω0 τ . (B.14)~ q0i

The portion involving atomic projection operators will now be considered. It is shownthat the projection-operator factors in front of the first and third terms are identical,
and can thus be factored out. Expanding out the projection-operator portion of the

131first term,(−1)q Pg Dq Pei D−q Pg σ(t)Pg =X|Jg IFg Mg i hJg IFg Mg | Dq |Je IFei Mei i

Mg ,Mg0 ,Mei× (−1)q hJe IFei Mei | D−q |Jg IFg Mg0 i hJg IFg Mg0 | σ(t)Pg|hJg ||D||Je i|2 (i) (i) †Aq Aq Pg σ(t)Pg
2Je + 1|hJg ||D||Je i|2(i) †=Pg σ(t)Pg A(i)q Aq
2Je + 1X
Pg σ(t) |Jg IFg Mg i hJg IFg Mg | Dq |Je IFei Mei i==(B.15)

Mg ,Mg0 ,Mei× (−1)q hJe IFei Mei | D−q |Jg IFg Mg0 i hJg IFg Mg0 |=(−1)q Pg σ(t)Pg Dq Pei D−q Pg ,where the final line is the projection-operator portion of the third equation. In obtaining
(i)

(i) †

the previous result, the fact that the product of raising and lowering operators Aq Aq

is proportional to the ground-state internal atomic projection operator Pg , and thus
commutes with the ground-state sub-density-matrix σgg has been utilized. Combining
the first and third terms of Eq. (B.12), the projection-operator terms can be factoredout,
1 X X |hJg ||D||Je i|2 (i) (i) †Aq Aq Pg σ(t)Pg− 2~ q
2Je + 1
i

Z ∞

´
³
(i)
(i)
dτ g(τ )e−iω0 τ + g(−τ )eiω0 τ

0

1 X X |hJg ||D||Je i|2 (i) (i) †=− 2
Aq Aq Pg σ(t)Pg~ q
2Je + 1iZ ∞−∞X X |hJg ||D||Je i|
1
(i) †
A(i)
q Aq Pg σ(t)Pg
3
3
3²0 (2π) c ~ q
2Je + 1
2=−i

X X |hJg ||D||Je i|
1
(i) †
A(i)
q Aq Pg σ(t)Pg
3²0 (2π)3 c3 ~ q
2Je + 12=−(i)dτ g(τ )e−iω0 τZdω ω 3Z ∞(i)dτ e−i(ω+ω0 )τ

−∞

Z
(i)

dω ω 3 δ(ω + ω0 ).

i(B.16)
It is clear that the contributions from the first and third terms of Eq. (B.12) vanish,
since the argument of the delta function can never be zero. Physically, this reflects thefact that the atom cannot make a transition from the ground state to the excited stateby emitting a photon.

132This leaves only the second and fourth terms,
1 XX
σ̇gg (t) = 2(−1)q Pg Dq Pei σ(t)Pej D−q Pg
~ q
i,j
Z ∞ ³
´
(i)
(j)
dτ g(τ )eiω0 τ + g(−τ )e−iω0 τ×0XX(−1)q Pg Dq Pei σ(t)Pej D−q Pg
=
q

i,j

µZ
¶ (B.17)
Z
(i)
(j)
~
3 −iωτ iω0 τ
3 iωτ −iω0 τdω ω ee+ dω ω e edτ
3²0 (2π)3 c30XX1=(−1)q Pg Dq Pei σ(t)Pej D−q Pg
3
3
3²0 (2π) c ~ q
i,j
µZ ∞
¶
Z
Z ∞
(i)
(j)
× dω ω 3dτ e−iωτ eiω0 τ +dτ eiωτ e−iω0 τ .
1
× 2
~

Z ∞

0

Now, using the fact that

0

R∞
0

dτ e−iωτ = πδ(ω) + iP ω1 , where the P operator indicatesthat the principal part must be taken when an integration is performed. It is seenfrom the last equation in Eq. (B.17) that the integrand of the ω integral will consist(i)of two delta functions, on each for ω − ω0(j)and ω − ω0 , along with two principal-

part terms, which combine to be of the form P(

1
1
(i) −
(j) ).
ω−ω0
ω−ω0

The principal-partexpression gives the additional energy splitting between the ith and the jth excitedstate manifolds due to interaction with the photon reservoir accompanied by the usual
pathological divergences which will be ignored in the present case. Neglecting thisenergy splitting and considering only the delta-function terms, as well as making the(i)approximation that the frequencies {ω0 } are approximately equal,σ̇gg (t) ∼=XX1(−1)q Pg Dq Pei σ(t)Pej D−q Pg
3
3
3²0 (2π) c ~ q

Z

i,j

XX
ω0 3
(−1)q Pg Dq Pei σ(t)Pej D−q Pg ,=3²0 (2π)3 c3 ~ qi,j(i)dω ω 3 δ(ω − ω0 )(B.18)

where ω0 is the average of the atomic Bohr frequencies ω0 . Expanding the projection-

133
operator portion,(−1)q Pg Dq Pei σ(t)Pej D−q Pg
X=(−1)q |Mg i hMg | Dq |Mei i hMei | σ(t) |Mej i hMej | D−q |Mg0 i hMg0 |

Mg ,Mg0 ,Mei ,Mej

X

=

£
¤
[hMg |Dq |Mei i |Mg i hMei |] σ(t) (−1)q hMej |D−q |Mg0 i |Mej i hMg0 |

Mg ,Mg0 ,Mei ,Mej

− hJg ||D||Je i hJe ||D||Jg i (i)†Aq σ(t)A(j)q
2Je + 1|hJg ||D||Je i|2 (i)†=Aq σ(t)A(j)q ,2Je + 1=

(B.19)
where the fact that hJg ||D||Je i† = (−1)Je −Jg hJe ||D||Jg i has been used. Equation (B.12)
finally becomesσ̇gg (t) =|hJg ||D||Je i|2 X X (i)ω0 3†
Aq σ(t)A(j)q3
3
3²0 (2π) c ~ 2Je + 1
q
i,j

∼
=γ

XX
q

(B.20)(j) †A(i)q σ(t)Aq ,i,jwhere γ = γJe →Jg is the transition rate from the Je excited state to the Jg ground state
and which are assumed here to be equal to the transition rate for all of the excited-state
hyperfine manifolds. It is clear that the approximation that the various internal atomic(i)frequencies ω0 are all equal to some average frequency ω0 amounts to assuming thateach excited-state manifold relaxes with the same spontaneous-emission rate.
Turning now to the excited-state sub-density-matrices,Z ∞ X
1σ̇ei ej (t) = −iδωei ej Pei σ(t)Pej − 2
dτ
(−1)q
~ 0
q
(
"
#
X
(k)
(i)
× g(τ )Pei Dq Pg D−q Pek σ(t)Pej eiω0 τ − Pei Dq Pg σ(t)Pg D−q Pej e−iω0 τk+g(−τ )"X#)
Pei σ(t)Pek Dq Pg D−q Pej e(k)−iω0 τ

− Pei Dq Pg σ(t)Pg D−q Pej e

(j)
iω0 τ

,

k

(B.21)

134
where δωab = ωa − ωb . The second and fourth terms will combine (having identical(i)projection-operator prefactors) to give two delta functions with arguments of ω + ω0(j)and ω +ω0 , respectively, which cannot be satisfied and thus give zero upon integration,
and also a principal part contribution of the form P(

1
1
(i) −
(j) ), which are neglected
ω+ω0
ω+ω0

as before. What remains is the first term and the third term,σ̇ei ej (t) = − iδωei ej Pei σ(t)Pej−

1 XX(−1)q Pei Dq Pg D−q Pek σ(t)Pej
~2 q
k

−

1 XX
(−1)q Pei σ(t)Pek Dq Pg D−q Pej
~2 q
k

Z ∞

(k)

dτ g(τ )eiω0 τ(B.22)

0

Z ∞

(k)

dτ g(−τ )e−iω0 τ .

0Simplifying the projection-operator prefactors,(−1)q Pei Dq Pg D−q Pek σ(t)Pej=

X(−1)q |Mei i hMei | Dq |Mg i hMg | D−q |Mek i hMek | σ(t)Pej

Mg ,Mei ,Mek=X

[hMei |Dq |Mg i |Mei i hMg |] [(−1)q hMg |D−q |Mek i |Mg i hMek |] σ(t)Pej

Mg ,Mei ,Mek=|hJg ||D||Je i|22Je + 1X

αFg ,Fei ,Mg ,Mei ,Jg ,Je ,I |Mei i hMg |Mg ,Mei ,Mek× αFg ,Fek ,Mg ,Mek ,Jg ,Je ,I |Mg i hMek | σ(t)Pej|hJg ||D||Je i|2 (i) † (k)=Aq Aq σ(t)Pej ,2Je + 1
(B.23)

135
and(−1)q Pei σ(t)Pek Dq Pg D−q Pej
X=(−1)q Pei σ(t) |Mek i hMek | Dq |Mg i hMg | D−q |Mej i hMej |Mg ,Mei ,MekX=

Pei σ(t) [hMek |Dq |Mg i |Mek i hMg |]Mg ,Mek ,Mej£¤
× (−1)q hMg |D−q |Mej i |Mg i hMej |=|hJg ||D||Je i|22Je + 1

X

(B.24)Pei σ(t)αFg ,Fek ,Mg ,Mek ,Jg ,Je ,I |Mek i hMg |Mg ,Mei ,Mek× αFg ,Fej ,Mg ,Mej ,Jg ,Je ,I |Mg i hMej |=|hJg ||D||Je i|2† (j)Pei σ(t)A(k)Aq ,q
2Je + 1the excited-state density-matrix equations of motion can be written asσ̇ei ej (t) = −iδωei ej Pei σ(t)Pej −´
γ X X ³ (i) † (k)† (j)Aq Aq σ(t)Pej + Pei σ(t)A(k)
A
q
q
2 q
k

= −iδωei ej Pei σ(t)Pej −

γ XX
2

q

³
´
† (l)
(k) † (l)Pei A(k)
A
σ(t)
+
σ(t)A
A
Pej ,
q
q
q
q

k,l

(B.25)
where again the principal-part contributions have been neglected and where in the lastline a trivial summation index has been added that will be useful later when combining
the various sub-density-matrix decay terms.
The final equations to be considered are the optical-coherence sub-density-matrices,
Z ∞ X
1
dτ
(−1)q
σ̇ei g (t) = −iδωei g Pei σ(t)Pg − 2
~ 0
q




X
X
(j)
(i)
Pei Dq Pg D−q Pej σ(t)Pg eiω0 τ −Pei Dq Pg σ(t)Pej D−q Pg e−iω0 τ × g(τ ) 

j
j



X
X
(j)
(j)
+g(−τ ) Pei σ(t)Pg Dq Pej D−q Pg eiω0 τ −Pei Dq Pg σ(t)Pej D−q Pg e−iω0 τ  .
jj(B.26)

136
It is important to mention here that the indices on the frequencies in the exponentials(i)are correct in the above equation. The frequency ω0 (not the summation index) should(j)only appear in the exponential in the second term, while ω0 ’s (the summation index)
should appear in all of the other terms. The second term and the fourth term are bothnon-secular, since they involve density-matrix elements σgei , and thus these terms canbe neglected. The third term is anti-resonant, i.e., it will contribute delta functions(i)with arguments ω + ω0 which cannot be satisfied. Again, principal-part contributions
from these terms have been neglected, yielding
σ̇ei g (t) = − iδωei g Pei σ(t)Pg −1 XX(−1)q Pei Dq Pg D−q Pej σ(t)Pg~2 qj= − iδωei g Pei σ(t)Pg −

Z ∞

(j)

dτ g(τ )eiω0 τ

0

γ XX
† (l)
Pei A(k)Aq σ(t)Pg .
q
2 q
k,l

(B.27)
Similarly,σ̇gei (t) = iδωgei Pg σ(t)Pei −

γ XX
† (l)
Pg σ(t)A(k)Aq Pei .
q
2 q(B.28)

k,lUsing Eq. (B.9), the equation of motion due to spontaneous emission for the fulldensity operator can be constructed,
σ̇(t) =´XX
γ X X ³ (i) † (j)
i
† (j)
(j) †
[σ(t), HA ] + γA(i)−Aq Aq σ(t) + σ(t)A(i).
q σ(t)Aq
q Aq
~
2 q
q
i,j

i,j

(B.29)
Defining the spontaneous emission relaxation operator,
Lsp [σ] = γ

XX
q

†

(j)
A(i)
−
q σ(t)Aq´
γ X X ³ (i) † (j)† (j)Aq Aq σ(t) + σ(t)A(i)
A
,
q
q
2 q

(B.30)i,ji,jthe equation of motion can be written asσ̇(t) =

i
[σ(t), HA ] + Lsp [σ].
~(B.31)

The atomic center-of-mass dependence that was ignored for simplicity in theprevious derivation can now be reintroduced. This amounts to adding an integral over137
momentum states in 3D that should have been included when the atomic projection
operators were inserted. Including this, the full relaxation operator takes the formshown in Eq. (3.4).

Appendix C
Semiclassical Master Equation for Atoms with Complex Internal
Structure

Treating the atom’s center-of-mass degree of freedom as a classical variable rather
than a quantum operator — often referred to simply as the “semiclassical approximation” despite being somewhat ambiguous — yields a set of dynamical equations called
the optical Bloch equations. This appendix provides a detailed derivation of the optical
Bloch equations corresponding to the situation discussed in Chapter 3: an atom withmultiple hyperfine manifolds of arbitrary spacing, interacting with an arbitrary threedimensional polarization-gradient laser field. These equations are then solved by direct
numerical means in Chapter 3 to obtain information such as the friction and diffusion
coefficients for laser cooling. This information, although often technically inferior to thesolutions of fully quantum equations, is usually much easier to calculate and can provide helpful qualitative insight into the dynamics from the use of the intuitive classicalconcepts of force, friction, and diffusion.

C.1

Equations without spontaneous decay

C.1.1

Atom-field interaction

The field can be written, in general, as
E(z, t) = E+ (z)e−iωL t + c.c.

(C.1)

139
For the σ+ σ− configuration, the positive frequency component is
E+ = E0 ²̂eikL z + E00 ²̂0 e−ikL z ,(C.2)

where E0 and E00 are the amplitudes of the two opposing, circularly polarized beams,
respectively, and ²̂ and ²̂0 are the respective polarization vectors, given by´
1 ³²̂ = ²̂+ ≡ − √ î + iĵ ,2´
1 ³²̂0 = ²̂− ≡ √ î − iĵ .
2

(C.3)
(C.4)

The interaction between the atom and the coherent field in the electric dipole
approximations is described by
V (t) = −D · E(z, t)(C.5)= E0 D(+1) e−i(ωL t−kL z) + E00 D(−1) e−i(ωL t+kL z)− E0 D(−1) ei(ωL t−kL z) − E00 D(+1) ei(ωL t+kL z) .
Using the Wigner-Eckart theorem, and noting that S commutes with D,
D(q) =

X|j, m, l, sihj, m, l, s| D(q) |j 0 , m0 , l0 , s0 ihj 0 , m0 , l0 , s0 |j,j 0 ,m,m0 ,l,l0 ,s,s0




0
j
1
j


0 0
=
|j, m, l, sihj 0 , m0 , l0 , s| (−1)j−m  hj, l, skDkj , l , si0j,j 0 ,m,m0 ,l,l0 ,s−m q mX=

X

0|j, m, l, sihj 0 , m0 , l0 , s| (−1)j−m+l+s+j +1

j,j 0 ,m,m0 ,l,l0 ,s







1
p
 j
 l j s0× (2j + 1)(2j + 1) hlkDkl0 i .
−m q m0 j 0 l0 1


j0

(C.6)
The properties of the 3-j symbol provide the information that the triangle relation mustbe satisfied between j,j 0 and 1, and that m0 − m = q. The 6-j symbol indicates thatthe following sets must each satisfy triangle relations: {j, l, s}, {j 0 , l0 , s}, {j, j 0 , 1}, and140{l, l0 , 1}. The first two are trivial, the third is redundant, and the fourth tells us that l0
can equal l + 1, l, and l − 1, however, parity considerations rule out the l0 = l possibility.
Now specialize to the case, useful to laser cooling, of a ground state jg coupled toan excited state multiplet je = jg − 1, jg , jg + 1, resulting from a system with l = lg inthe ground state and l = le = lg + 1 in the excited state, coupled with a spin of s = 1.
The dipole matrix element can then be written asX·(q)|je , mg + q, lg + 1ihjg , mg , lg |D = Degje ,mgq
× (−1)je +jg −mg +lg −q+1 (2jg + 1)(2je + 1) 






je
1 jg  lg + 1 je 1



lg 1
−mg − q q mg  jg

− |jg , mg , lg ihje , mg − q, lg + 1|




¸
q

1
je  lgjg
1
 jg
je +jg −mg +lg
× (−1)(2jg + 1)(2je + 1) 
,





−mg q mg − q
je lg + 1 1(C.7)
which can be further simplified to yield

D(q) =Deg

X(−1)je +jg −mg +lg

q(2jg + 1)(2je + 1)



 lg

jg



1


je lg + 1 1
·
1 jg 
 je
× |je , mg + q, lg + 1ihjg , mg , lg | (−1)q−1 −mg − q q mg1
je  ¸ jg− |jg , mg , lg ihje , mg − q, lg + 1| 
 ,−mg q mg − q
je ,mg

(C.8)

or
(q)

(q)

D(q) = Deg (L+ − L− ),(C.9)

where the Deg is the reduced dipole matrix element, given byDeg ≡ hle kDklg i = − hlg kDkle i ,(C.10)

141
and the raising and lowering operators are defined as
(q)

L+ ≡

X

|je , mg + q, lg + 1ihjg , mg , lg |je ,mg






q

1 jg  lg + 1 je 1
 je
je +jg −mg +lg −q+1
× (−1)(2jg + 1)(2je + 1) 
,





−mg − q q mg
jg
lg 1
(q)

L− ≡

X

|jg , mg , lg ihje , mg − q, lg + 1|

je ,mg






q
jg
1
1
je   lg
 jg
je +jg −mg +lg
× (−1).
(2jg + 1)(2je + 1) 
−mg q mg − q je lg + 1 1(C.11)

142
(q)

(−q)

Note that (L± )† = L∓(q)(L+ )† =X, as can be seen for the raising operator:|jg , mg , lg ihje , mg + q, lg + 1|

je ,mgq
1 jg  lg + 1 je 1
 je
je +jg −mg +lg −q+1
× (−1)(2jg + 1)(2je + 1) ,
−mg − q q mg  jg
lg 1
=

X



|jg , mg , lg ihje , mg + q, lg + 1|

je ,mgq
× (−1)je +jg −mg +lg −q+1 (2jg + 1)(2je + 1) =X|jg , mg , lg ihje , mg + q, lg + 1|

je ,mgq
× (−1)je +jg −mg +lg −q+1 (2jg + 1)(2je + 1) 

=

X






je
1 jg  je lg + 1 1
,





−mg − q q mg
lg
jg
1






je
1 jg  lg
jg
1
,





−mg − q q mg
je lg + 1 1|jg , mg , lg ihje , mg + q, lg + 1| (−1)je +1+jg

je ,mg






q

je
jg
1
 jg 1 lg
je +jg −mg +lg −q+1
(2jg + 1)(2je + 1) ,× (−1)

mg q −mg − q je lg + 1 1=X

|jg , mg , lg ihje , mg + q, lg + 1| (−1)2(je +1+jg ) (−1)−q+1

je ,mg× (−1)je +jg −mg +lg

(−q)

= L+


q
 jg
(2jg + 1)(2je + 1) 
−mg





1
je   lg
jg
1
,



−q mg + q je lg + 1 1,(C.12)

where the fact that (−1)2(je +1+jg ) = 1 and that (−1)−q+1 = 1 for q = ±1 has beenused. The procedure for the lowering operator is similar.

143
Equation (C.5), written in terms of raising and lowering operators, becomes
· ³
´
³
´
(+1)
(−1)
(+1)
(−1)
V (t) = Deg E0 L+ − L−e−i(ωL t−kL z) + E00 L+ − L−e−i(ωL t+kL z)
¸
³
´
³
´
(−1)
(−1)
(+1)
(+1)
i(ωL t−kL z)
0i(ωL t+kL z)− E0 L+ − L−e
− E0 L+ − L−e.
(C.13)In this form, the rotating wave approximation can easily be made by removing terms
(q)

(q)

of the form L± e±ωL t , and keeping only terms of the form L± e∓iωL t . The resulting
expression for the interaction is(+1) −i(ωL t−kL z)

V = Deg [E0 L+e(−1) −i(ωL t+kL z)+ E00 L+e

(C.14)(−1)(+1)+ E0 L− ei(ωL t−kL z) + E00 L− ei(ωL t+kL z) ],or, if E0 = E00 ,V =i~Ω h (+1) −i(ωL t−kL z)
(−1)
(−1)
(+1)
L+ e
+ L+ e−i(ωL t+kL z) + L− ei(ωL t−kL z) + L− ei(ωL t+kL z) ,2
(C.15)

where Ω ≡ 2Deg E0 /~ is the reduced Rabi frequency.
C.1.2
C.1.2.1

Transformation to moving, rotating frame
General properties of time-dependent transformations

For a Hamiltonian H, the equation of motion is given byi~ρ̇ = [H, ρ] .

(C.16)The goal is to find the equation of motion in some other coordinate system, with a timedependent transformation given by T (t). The density operator in the new coordinate
system is given byρ̃ = T (t)ρT † (t),(C.17)

144
which has a time derivative of´d ³
T ρT †dt
dT †dT= i~T ρ̇T † + i~ ρT † + i~T ρ
dt
dt
h
i
= i~T ρ̇T † − [H, T ] ρT † − T ρ H, T †hi= i~T ρ̇T † − [H, ρ̃] + T HT † , ρ̃ .

i~ρ̃˙ = i~

(C.18)

The above result can then be used in Eq. (C.16) with a transformation applied,i~T ρ̇T † = T [H, ρ] T † ,(C.19)hii~ρ̃˙ = T HT † , ρ̃ − [H, ρ̃] + T [H, ρ] T † .

(C.20)

giving

The second two terms can be manipulated,
T [H, ρ] T † − [H, ρ̃] = T HρT † − T ρHT † − HT ρT † + T ρT † Hhi= −T ρ H, T † − [H, T ] ρT †¸·
dT † dT †+ρT= i~ T ρ
dt
dt
·
¸
†
dT
dT †
= i~ ρ̃T+
T ρ̃ ,
dt
dt(C.21)·¸hidT † dT ††˙i~ρ̃ = T HT , ρ̃ + i~ ρ̃T+T ρ̃ .
dt
dt(C.22)

resulting in

Also, noting that
µ
¶
i
1 h
−
H, T †i~³´
1=−
T HT † − Hi~´
1 ³
H − T HT †=i~
1=
[H, T ] T †i~dT= − T †,
dt

dT †
=T
T
dt

(C.23)

145
yields
· †
¸
h
i
T
dT †
†
˙
i~ρ̃ = T HT , ρ̃ + i~
T ρ̃ − ρ̃ T
dt
dt
¸
h
i · dT= T HT † , ρ̃ + i~ T † , ρ̃ ,dt

(C.24)orhii~ρ̃˙ = H̃, ρ̃ ,H̃ ≡ T HT † + i~

C.1.2.2

dT †
T .
dt

(C.25)(C.26)

Application to an atom in the σ+ σ− laser configuration

Since the coupling of the atom to the vacuum field is being ignored for the time
being, the total Hamiltonian is
H = H0 + V (t),(C.27)

where H0 is the internal atomic Hamiltonian and V (t) represents the interaction between
the atom and the coherent external laser field. If the transformation to the center-ofmass frame of an atom moving with velocity v is made, z must be replaced with vt. In
this moving frame, the interaction term is of the form
V (t) ∼ d · ²̂y = dx sin kvt + dy cos kvt.

(C.28)

It is useful to transform to a rotating coordinate system in order to simplify the interaction term. Consider a transformation operator given by
T (t) = exp(−ikvtJz /~).

(C.29)

Using the commutation relations [Jz , dx ] = i~dy and [Jz , dy ] = −i~dx , the transformedinteraction operator takes the simplified formṼ (t) = T (t) [dx sin kvt + dy cos kvt] T † (t) = dy .

(C.30)

146
The time-independent form of the interaction operator in the rotating frame allows
steady-state solutions to be obtained when working in this reference frame. With theresults of the previous section, the equation of motion ishii~ρ̃˙ = H̃, ρ̃ .

(C.31)The inertial term in Eq. (C.26), with the transformation operator given in Eq. (C.29),
becomesdT †T = kvJz ,dt(C.32)

H̃ = H0 + kvJz + Ṽ (t),(C.33)i~and sosince H0 commutes with Jz . In the rotating frame, the interaction operator from
Eq. (C.15) becomes
·
~Ω
(+1)
(−1)
T (t) L+ e−i(ωL t−kL z) + L+ e−i(ωL t+kL z)
Ṽ (t) =
2
(−1)
(+1)
+ L− ei(ωL t−kL z) + L− ei(ωL t+kL z)¸T † (t)·~Ω (+1) −i(ωL t+kL vt−kL z)
(−1)
=
L+ e
+ L+ e−i(ωL t−kL vt+kL z)
2

¸

(−1)
(+1)
+ L− ei(ωL t+kL vt−kL z) + L− ei(ωL t−kL vt+kL z)(C.34)
Here it is useful to compare to the form of the interaction term obtained by
Cohen-Tannoudji and Dalibard [16], in which cooling is performed on a J = 1 groundstate (with sublevels |gi i) and a J = 2 excited state (with sublevels |ei i). For thissituation, Eq. (C.34) becomes
·µ
¶
~Ω
1
1
Ṽ (t) = √|e+2 ihg+1 | + √ |e+1 ihg0 | + √ |e0 ihg−1 | e−i(ωL t+kL vt−kL z)
2 3
2
6
¶
µ
1
1
+ |e−2 ihg−1 | + √ |e−1 ihg0 | + √ |e0 ihg+1 | e−i(ωL t−kL vt+kL z)
2
6(C.35)µ¶1
1i(ωL t+kL vt−kL z)+ |g+1 ihe+2 | + √ |g0 ihe+1 | + √ |g−1 ihe0 | e
2
6
µ
¶
¸
1
1
i(ωL t−kL vt+kL z)+ |g−1 ihe−2 | + √ |g0 ihe−1 | + √ |g+1 ihe0 | e
,
2
6

.

147
or, in the notation of Cohen-Tannoudji and Dalibard,
·
¸
~Ω
1
1
√
√
√
Ṽ (t) =|g+1 ihe+2 | +|g0 ihe+1 | +|g−1 ihe0 | ei(ω+ −kL z) + h.c
2 3
2
6
·
¸
~Ω
1
1
+ √ |g−1 ihe−2 | + √ |g0 ihe−1 | + √ |g+1 ihe0 | ei(ω− +kL z) + h.c,
2 3
26

(C.36)

where
ω± = ωL ± kL v.

(C.37)√
From Eq. (C.36), it appears that Eq. (5.1) of Ref. [16] is missing a factor of 1/ 3 andthe ω± should be reversed. Furthermore, Eq. (5.1) of Ref. [16] is misleading because thetransformation to the moving frame has not been made. In applying this transformation (z → vt) to my Eq. (C.34), the velocity terms in the exponentials cancel exactly,
resulting inṼ (t) =C.2i~Ω h (+1) −iωL t
(−1)
(−1)
(+1)
L+ e
+ L+ e−iωL t + L− eiωL t + L− eiωL t .
2(C.38)

Optical Bloch equations
The equations of motion developed up to this point can now be extended toinclude spontaneous emission. The form of the spontaneous-emission relaxation operator is derived in detail in Appending B and the results will be used here, with thecenter-of-mass degree of freedom neglected due to the semiclassical nature of the current
equations. There will be three classes of equations:
³
´ h
i
(m)
i~ρ̃˙ ei , gj = H̃, ρ̃  (m)eihii~ρ̃˙ (gi , gj ) = H̃, ρ̃,gj(gi ,gj )

,

,

³
´ h
i
(m) (n)i~ρ̃˙ ei , ej
= H̃, ρ̃  (m) (n)  .
ei,ej(C.39)(C.40)(C.41)

The first class describes optical coherences, the second describes ground-state populations and ground-state coherences, and the third describes excited-state populations
and excited-state coherences.

148
C.2.1

Class I: optical coherences

Redefining the density matrix elements corresponding to optical coherences as
³
´
³
´
(m)
(m)
ρ̃ ei , gj → ρ̃ ei , gj eiωL t ,
³
´
³
´
(m)
(m)
ρ̃ gi , ej
→ ρ̃ gi , ej
e−iωL t .

(C.42)(C.43)All other elements are unchanged. This results in a change to the optical-coherence
equations of motion of
³
´
³
´
´
³
(m)
(m)
(m)
ρ̃˙ ei , gj = iωL ρ̃ ei , gj + ρ̃˙ (old) ei , gj eiωL t ,
´
³
´
³
´
³
(m)
(m)
(m)
e−iωL t .
ρ̃˙ gi , ej
= −iωL ρ̃ gi , ej
+ ρ̃˙ (old) gi , ej

(C.44)(C.45)

³
´
³
´ i
h
i
(m)
(m)
ρ̃˙ ei , gj = iωL ρ̃ ei , gj − eiωL t H0 + kL vJz + Ṽ , ρ̃  (m) ~ei ,gj
³
´ i
h
i
(m)
= i [ωL − ωe(m) + ωg − kL v(i − j)] ρ̃ ei , gj − eiωL t Ṽ , ρ̃  (m) ~ei ,gj
h
i
h
i ³
´ i
(m)
= i δω (m) − ∆ − kL v(i − j) ρ̃ ej , gi − eiωL t Ṽ , ρ̃  (m)  .
~ei ,gj(C.46)
³
´
³
´ i
h
i
(m)
(m)
ρ̃˙ gi , ej
= −iωL ρ̃ gi , ej
− e−iωL t H0 + kL vJz + Ṽ , ρ̃  (m) ~
gi ,ej
³
´ i
h
i
(m)
− e−iωL t Ṽ , ρ̃  (m) = i [−ωL − ωg + ωe(m) − kL v(i − j)] ρ̃ gi , ej~
gi ,ej
h
i ³
´ i
h
i
(m)
= i ∆ − δω (m) − kL v(i − j) ρ̃ gi , ej
− e−iωL t Ṽ , ρ̃  (m)  .
~gi ,ej(C.47)

hieiωL t Ṽ , ρ̃ (m)
(ei

149
¯
¯ E D
¯
¯ Ei
hD
¯
¯
(m) ¯
(m) ¯
ei ¯ Ṽ ρ̃ ¯ gj − ei ¯ ρ̃Ṽ ¯ gj
,gj )·~Ω iωL t D (m) ¯¯ (+1) −iωL t ¯¯ E D (m) ¯¯ (−1) −iωL t ¯¯ E=eei ¯ L+ eρ̃ ¯ gj + ei ¯ L+ eρ̃ ¯ gj
2
¯
¯
¯ E D
¯ E¸
D
(m) ¯
(+1) −iωL t ¯(m) ¯(−1) −iωL t ¯− ei ¯ ρ̃L+ e¯ gj − ei ¯ ρ̃L+ e¯ gj·~Ω D (m) ¯¯ (+1) ¯¯ E D (m) ¯¯ (−1) ¯¯ E
ei ¯ L+ ρ̃ ¯ gj + ei ¯ L+ ρ̃ ¯ gj
=
2
¯
¯ E D
¯
¯ E¸
D
(m) ¯
(+1) ¯(m) ¯(−1) ¯− ei ¯ ρ̃L+ ¯ gj − ei ¯ ρ̃L+ ¯ gj .
= eiωL t(C.48)The raising operators in the four matrix elements above will produce terms of theform hgi |ρ̃|gj i and hei |ρ̃|ej i, both of which are not redefined according to Eqs. (C.42)
˙ i , e(m) ), an exponential of e−iωL tand (C.43). Similarly, in the equation of motion for ρ̃(gjwill multiply the interaction commutator, and this time will cancel the eiωL t exponentials
in the lowering operator portions of the interaction operator.
At this point, further confirmation can be made of the error in Ref. [16]. E.g.,
if Eq. (C.46) is used to determine the coefficient of ρ̃ (g0 , g0 ) in the equation of motionfor ρ̃˙ (e1 , g0 ), one finds that Cohen-Tannoudji and Dalibard’s result differs by lacking a√factor of 1/ 3, as we noted before. Using the paper by Bambini and Agresti (Ref. [43]),
Eq. (2.8b), agreement with my result is found. (Note: Bambini and Agresti define their
reduced dipole matrix element as hlg kdkle i, which differs from mine by a sign.)
C.2.2

Class II: ground-state populations and ground-state coherencesiihṼ , ρ̃.
ρ̃˙ (gi , gj ) = −ikL v(i − j)ρ̃ (gi , gj ) −~(gi ,gj )(C.49)

150
C.2.3

Class III: excited-state populations and excited-state coherences

i
³
´
³
´ ih
(m) (n)(m) (n)Ṽ , ρ̃  (m) (n) ρ̃˙ ei , ej
= i [ωe(n) − ωe(m) − kL v(i − j)] ρ̃ ei , ej−~
ei ,ejihi ³´ ih(m) (n)Ṽ , ρ̃  (m) (n)  .
= i δω (m) − δω (n) − kL v(i − j) ρ̃ ei , ej−~
ei ,ej(C.50)