Implosions as a probe for fluctuations in Bose-EinsteincondensatesbyMichael Van de Graaff
B.A., Physics, Lawrence University, Appleton, WI. 2014

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
2022
Committee Members:
Eric A. Cornell, ChairJohn Bohn
Jun YeJose D’Incao
Carl Lineberger

ii
Van de Graaff, Michael (Ph.D., Physics)
Implosions as a probe for fluctuations in Bose-Einstein condensates
Thesis directed by Prof. Eric A. Cornell

Quantum and thermal effects play important roles in the fluctuations of ultracold Bose gases,
yet the fluctuations remain difficult to probe experimentally in 3D bulk systems. We seek toamplify these fluctuations via a quench from positive to negative scattering length. Following sucha quench, initial fluctuations in the gas, whether seeded by a lattice potential, thermal fluctuations,or quantum fluctuations, are amplified on a time-scale which is fast compared to the bulk collapseof the cloud. We perform experiments on clouds seeded with an initial, well-understood density
perturbation due to an optical lattice and on clouds with no externally imposed fluctuations. We
compare our observations to a variety of theoretical models, which predict the static structure
factor as a function of time and momentum in the collapsing gas.

Dedication

To my parents.

iv

Acknowledgements

I owe thanks a great number of people in my life for all the help and support without
whom I would never have completed this degree. That list is very long, and I cannot fit them
all here. My parents always cultivated my sense of curiousity and have supported me financiallyand emotionally throughout my life. I have many fond memories of attending the public physics
lectures at the University of Chicago with my dad, who was himself a physics major in college.
That he passed away before he could see his son finally acheive the PhD I sought from the young
age of ten or eleven is one of the greatest sadnesses of my life. I miss him very much.
I have lived an exceptionally fortunate life in almost every respect. It was not only that the
mentors, teachers, and professors in my life have genuinely cared about my happiness and future
success, but also that they were and are extremely competent in their roles. I was taught by goodteachers, not simply well meaning ones. I was raised by wise parents, not simply loving ones.
I began working at JILA on July 14, 2014. For my first two years I worked with fellow
graduate students Ming-Guang Hu and then Dhruv Kedar and was advised by Debbie Jin and Eric
Cornell. Tradgically Debbie passed away in 2016 and the world was robbed of a great scientist.
I remember the fear and uncertainty I experienced in the aftermath, not knowing what the restof my graduate career would look like. Fortunately, it went about as well as could be expected
and I continued working with Eric and my new co-asdvisor Jun Ye on new experiments with thepotassium-39 machine which had previously been Debbie’s Fermi gas lab. From 2016 to early 2020
I had the pleasure of working with a number of talented students, particularly fellow graduate
student Roman Chapurin, Xin Xie, and Noah Schlossberger, as well as undergraduates Carlos

v
Lopez-Abadia, Jared Popowski, and Bjorn Sumner. However, it was not long after the pandemic
began in earnest that I found myself the only student in the lab, working on a new topic quite
different from the spectroscopy based measurements I was used to.
Although I have been the only experimental graduate student working on the topics of thisthesis for the last two years, I have been aided by fellow student Eli Halperin, who has excellent
theoretical and numerical work on implosions critical to this thesis. He has rigourously attacked
the subject from numerous directions in an effort to understand our observations, and without his
work this thesis would be sorely lacking. Eli is also a funny and awesome guy whom I shall miss as
I leave JILA. Eli’s advisor John Bohn has also been a source of advice for me; his involvment has
also given me a sense of connection and continuity with Debbie Jin, my first advisor and John’slate wife, who’s signature cannmot be on this thesis.
Thank you, everyone. I couldn’t have done it without you.

vi

Contents

Chapter
1 Introduction

1

1.1

The world of ultracold interacting Bose gases . . . . . . . . . . . . . . . . . . . . . .

1

1.2

Timeline of scientific investigation of interacting Bose gases . . . . . . . . . . . . . .

3

1.2.1

Before 1995 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.2.2

Interactions and early BEC . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.2.3

Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.3

A brief overview of my career . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.4

Vision for a fluctuation microscope . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10

1.5

1.4.1

Wavefunction implosion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

1.4.2

Mode stability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.4.3

Calibrating fluctuations measurement . . . . . . . . . . . . . . . . . . . . . . 17

Contents of this Thesis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2 Apparatus

21

2.1

Overview of cooling procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21

2.2

Lasers for cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

2.3

Magnetic field control and trapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.3.1

High current coils. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.3.2

Low current science coils . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

vii
2.4

State Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.5

Hyperfine transitions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.6

Optical Dipole trap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.6.1

2.7

Trap characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35

Optical lattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3 Image Apparatus and Analysis
3.1

3.2

Absorption Imaging: the Basics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
3.1.1

Ideal Beer’s Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42

3.1.2

Absorption Imaging in the Real World: Bad Light . . . . . . . . . . . . . . . 43

3.1.3

Absorption Imaging in the Real World: Saturation . . . . . . . . . . . . . . . 44

3.1.4

Absorption imaging In the real world, final form . . . . . . . . . . . . . . . . 47

Side Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2.1

3.3

3.4

Design of vertical imaging system . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.3.2

Top camera calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

Fringe removal with PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
Nature of fringes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

Fitting the Images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
Column densities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Potential problems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

4 Fluctuation Theory of BEC
4.1

4.2

. . . . . . . . . . . . . . . . . . . . . . . . . 48

3.3.1

3.5.1
3.6

Calibrating Side probe intensity

Top imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

3.4.1
3.5

40

65

Criterion for BEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
4.1.1

Uniform systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.1.2

Trapped systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Presence and absence of BEC in 1, 2 and 3 dimensional ideal gases . . . . . . . . . . 68

viii

4.3

4.2.1

General facts of Bose systems . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.2.2

Foreshadowing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.2.3

Uniform ideal system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

4.2.4

Harmonically trapped ideal systems . . . . . . . . . . . . . . . . . . . . . . . 75

Weakly interacting Bose gases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.3.1

Dilute and cold gases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 79

4.3.2

Higher order approximations: Bogoliobov approximation and excitation spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

4.4

4.3.3

LHY corrections . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

4.3.4

The GP equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

BEC, Superfluidity, and Fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
4.4.1

Landaus criterion for superfluidity . . . . . . . . . . . . . . . . . . . . . . . . 89

4.4.2

The role of the condensate phase . . . . . . . . . . . . . . . . . . . . . . . . . 91

4.4.3

Quantum Hydrodynamics and the phase correlation function . . . . . . . . . 92

4.5

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

4.6

The Static Structure Factor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96

4.7

4.6.1

Static Structure factor for a > 0 . . . . . . . . . . . . . . . . . . . . . . . . . 96

4.6.2

Evolution of structure factor at a < 0 . . . . . . . . . . . . . . . . . . . . . . 98

Related work in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

5 Controlled Demolition

102

5.1

Optical Lattice overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5.2

Lattice geometries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5.3

5.2.1

X-lattice geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 102

5.2.2

Y-lattice geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

Optical lattice characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.3.1

Absolute lattice depth measurements of deep lattices . . . . . . . . . . . . . . 104

ix

5.4

5.5

5.3.2

Deep Lattice calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 106

5.3.3

Shallow lattice calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Lattice implosion experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111
5.4.1

Lattice sanity checks . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.4.2

Persistent lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5.4.3

Vanishing lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

Contrast correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

6 Intrinsic Fluctuations
6.1

128

Discussion of sources of error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128
6.1.1

Shall I compare BEC to summer’s day? . . . . . . . . . . . . . . . . . . . . . 137

6.2

Sample preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

6.3

Computation of the structure factor . . . . . . . . . . . . . . . . . . . . . . . . . . . 140

6.4

Contrast correction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

6.5

Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144

6.6

6.5.1

Calibrated wavenumbers . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

6.5.2

Uncorrected results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148

Other potential sources of error . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155
6.6.1

Miscalibrated conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 155

6.6.2

Spin purity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

6.6.3

Revisiting S0exp (k)/S0th (k) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 156

7 Summary and Outlook

159

x
Bibliography

162

Appendix
A Density of States

171

A.1 Density of states for a free particle . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171
A.2 Density of States for harmonic confinement . . . . . . . . . . . . . . . . . . . . . . . 173
B Lattice potentials

175

B.1 Interference of two laser fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 175
B.2 Partial derivation of Raman-Nath approximation . . . . . . . . . . . . . . . . . . . . 176
C Thermal distribution of Bosons

177

C.0.1 Time of flight analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 178

xi

List of Tables

Table
5.1

Conditions for seeded implosions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 116

6.1

Conditions for unseeded implosions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139

xii

List of Figures

Figure1.1

Growth and dissipation of a density perturbation in a uniform gas . . . . . . . . . . 11

1.2

Illustration of amplification . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

1.3

Force balancing in GPE . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

1.4

Dispersion relation for positive and negative a . . . . . . . . . . . . . . . . . . . . . . 16

1.5

Real space and Fourier space density distributions. . . . . . . . . . . . . . . . . . . . 18

2.1

Energy level diagram for 39 K. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2

39 K S
1/2 Manifold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.3

39 K P
3/2 manifold.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.4

39 K P
1/2 manifold. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.5

Magnetic field control enclosure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.6

Real and imaginary parts of the intra-species scattering length . . . . . . . . . . . . 28

2.7

|1, −1i → |2, −2i transition measurements . . . . . . . . . . . . . . . . . . . . . . . . 30

2.8

A photo of the science cell area. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

2.9

H1 and H2 beams breadboard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

2.10 Top Breadboard . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.11 Top probe and ODT V beam path through science cell . . . . . . . . . . . . . . . . . 34
2.12 Calibration of Z-axis trap frequency . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.13 Calibration of X and Y axis trap frequencies . . . . . . . . . . . . . . . . . . . . . . . 37

xiii
2.14 Y and Z final lattice optics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
3.1

Saturation of OD due to bad light . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.2

Side imaging system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

3.3

Optical depth vs. probe intensity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.4

image of beam used for calibratiing probe intensity. right, camera signal vs probe
power measured via powermeter. linearity is good

. . . . . . . . . . . . . . . . . . . 50

3.5

Example absorption images . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.6

Top imaging system used for in-situ imaging. . . . . . . . . . . . . . . . . . . . . . . 52

3.7

The objective for the high magnification imaging system . . . . . . . . . . . . . . . . 53

3.8

Left, out of focus, but fringes on top of cloud. middle, cloud size vs objective Z.
right focused bulk cloud, but fringes far away.

. . . . . . . . . . . . . . . . . . . . . 54

3.9

OD reconstruction with PCA . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5.1

Geometry of X-lattice and Y-lattice

5.2Time-of-flight images of excited lattice modes . . . . . . . . . . . . . . . . . . . . . . 107

5.3

Lattice calibration using Bragg diffraction . . . . . . . . . . . . . . . . . . . . . . . . 108

5.4

Lattice peaks in Fourier space . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

5.5

Calibration of the in-situ X-lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 111

5.6

Perturbation response to turning off the lattice. . . . . . . . . . . . . . . . . . . . . . 113

5.7

Implosion with persistent X-lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

5.8

Dispersion relation and timescales for persistent lattice implosions . . . . . . . . . . 115

5.9

X-lattice growth vs initial depth . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115. . . . . . . . . . . . . . . . . . . . . . . . . . . 103

5.10 Vanishing Y-lattice Perturbation depth vs implosion time . . . . . . . . . . . . . . . 117
5.11 Vanishing X-lattice Perturbation depth vs implosion time . . . . . . . . . . . . . . . 118
5.12 Vanishing Y-lattice Perturbation depth vs implosion time . . . . . . . . . . . . . . . 119
5.14 Phantom lattice leaking into Y-lattice . . . . . . . . . . . . . . . . . . . . . . . . . . 121
5.17 Ideal Fourier signal from numerics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125

xiv
5.18 X-lattice signal fraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126
5.19 Ratio of ideal sensitivity to observed sensitivity for the X-lattice Fourier component 126
6.1

Structure Factor vs scattering length . . . . . . . . . . . . . . . . . . . . . . . . . . . 130

6.2

FFT without atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 134

6.3

FFT with atoms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

6.4

Fourier space signal with and without atoms . . . . . . . . . . . . . . . . . . . . . . . 136

6.7

Observed structure factor vs implosion time for the k = 1.7 ring in the high temperature, high density data set . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 145

6.8

s . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146

6.9

Initial structure factor vs. wavenumber. . . . . . . . . . . . . . . . . . . . . . . . . . 147

6.10 Initial structure factor vs. wavenumber. . . . . . . . . . . . . . . . . . . . . . . . . . 149
6.12 Low temperature High density implosions . . . . . . . . . . . . . . . . . . . . . . . . 151
6.13 High temperature High density implosions . . . . . . . . . . . . . . . . . . . . . . . . 152
6.14 Low temperature Low density implosions . . . . . . . . . . . . . . . . . . . . . . . . 153
6.15 High temperature Low density implosions . . . . . . . . . . . . . . . . . . . . . . . . 154

Chapter 1
Introduction

1.1

The world of ultracold interacting Bose gases
The field of ultracold atoms has flourished tremendously in the last several decades, par-ticularly after the development of laser cooling, magnetic trapping, and evaporative cooling. The
first observation of Bose-Einstein condensation (BEC) was achieved right here in JILA in 1995 ina gas of 87 Rb and was quickly followed by a veritable explosion of work done with BECs which has
continued to this day [6, 41]. Among the many things learned during this time was the importanceof interatomic interactions. At very low energies, the interactions between particles may be characterized by a parameter called the s-wave scattering length: as . I will generally omit the subscripts and use a to denote the s-wave scattering length. To oversimplify, when a > 0, the interaction
between atoms is repulsive, and when a < 0 they are attractive. The magnitude of a determines
the interaction strength between particles [98, 101].
The interaction strength of various atomic species varies widely [32]. It is impossible tocalculate a from first principles for all but the simplest atoms. For instance, in a single atom of
87 Rb there are 37 electrons, 37 protons, and 50 neutrons, and that’s only mentioning particles

known about in the 1930s. Yet interactions are of extreme importance to the experimentalist
hoping to understand and exploit them. Certain species, such as 87 Rb, have interactions which
are unvarying over a large range of experimental parameters and which are also favorable forthe elastic interactions necessary for evaporative cooling, which is frequently required to achieve
quantum degeneracy. Other atoms, such as 85 Rb, have much stronger interactions which make2them prone to inelastic collisions, an impediment to evaporative cooling[32].
Yet the interactions are what make the few-body and many body physics questions interesting. A truly ideal gas is a pretty boring thing. With ultracold atoms, and with the right choice ofatom, we are fortunate enough to possess a tool allowing us to tune the strength |a| and sign of theinteractions via a Fano-Feshbach resonance (which I will abbreviate to Feshbach resonance for thesake of brevity)[32]. These Feshbach resonances arise due to the proximity of a bound molecularstate in a closed spin channel to the energy of two colliding particles in an open spin channel.
[52, 53, 50, 51, 49]
With the ability to tune the scattering length a in an ultracold dilute gas comes the opportunity to explore quantum mechanical effects which can only be accessed at such low temperatures.
Previous work in my lab has investigated the Efimov effect, which is a three-body effect resulting in
an infinite series of three-body bound states in a system with unitary interactions (a → ±∞)[45].
Observing this effect requires a system of three particles which interact with pairwise divergent
interactions in a coherent way. Phenomena like the Efimov effect (and the Feshbach resonanceitself) require these very cold, dilute conditions for the subtle quantum-mechanical effects to beobserved.
Today’s dilute ultracold atomic gasses are certainly not the only systems in which quantum
mechanical effects are important. Superconductivity was first observed in 1911 when mercury was
sufficiently cooled [93, 43]. Superfluid Helium, despite its sweltering temperature, is even moresimilar; liquefied 4 He exhibits viscosity-less flow up to a toasty 2.17 K. Clearly, one need not becold for quantum mechanics to be relevant. If 2.17 Kelvin does not seem hot to you, we can increase
the temperature very slightly to 300 K, where everyday metals contain a degenerate Fermi gas ofelectrons. Increasing our temperature a bit more to 106 K and we have the degenerate quantum
matter of a neutron star. It may seem unfair to call both superfluid helium at 2 K and neutronstars at 106 K hot, but since the condensates in our lab are on the order of 10−8 K, I think its fine.

3

1.2

Timeline of scientific investigation of interacting Bose gases

1.2.1

Before 1995

The prehistory of Bose gas physics began in 1924 when Satyendra Nath Bose derived Plank’sradiation law using only counting statistics and the quantized nature of photons [18]. Prior to Bose’s
work, all derivations of Planck’s law required appealing to classical physics; Bose’s derivation was
entirely quantum in nature. Bose sent his results to Albert Einstein following difficulty in publishing
his findings, and Einstein, recognizing the importance of this new insight, translated Bose’s work
into German and had it published on Bose’s behalf [46]. It is worth appreciating the fact that
Bose’s derivation fit quite comfortably on a single page, and in notation so similar today’s notationthat I found it only slightly difficult to follow the derivation even in German; it should be notedthat I do not speak a word of German.
Einstein extended Bose’s argument to an ideal gas of massive particles and clearly elucidated
the phenomena of condensation for an ideal gas at low temperatures, emphasizing how this follows
from Bose’s method of defining the state of the gas [47]. Einstein summarized Bose’s method“Following Bose: A states of the gas is microscopically defined by specifying how many molecules
are sitting in each [unit cell of configuration space]”. In contrast “According to the hypothesisof the statistical independence of the molecules: A state is microscopically defined by specifying
for each molecule the [unit cell of configuration space] in which is sits.” This is how Einstein
elucidated the difference between gases composed of indistinguishable (following Bose, wherein themolecules are not treated as ‘mutually statistically independent’) versus distinguishable particles(the molecules are treated as ‘mutually statistically independent’). Einstein’s detailed treatmentof both approaches in [47] makes clear how the subtle difference in defining the state of a gas has
dramatic impact on its properties.
However, while the work of Bose and Einstein feels very familiar to an ultracold atomic
physicist in the early 2020s, the phenomena of Bose-Einstein condensation would first be recognized within a system very unlike the dilute gases envisioned by Einstein and generated in today’s

4
ultracold atom labs: superfluid Helium. Helium had first been liquefied in 1908 [92], but efforts
to solidify it indicated that helium remained a liquid under its own pressure all the way down
to zero temperature. Since the entropy of a zero temperature system must be zero, this impliedthat somehow the liquid helium was undergoing an unknown ordering process. In 1937-1938, thephenomenon of superfluidity was experimentally verified by two groups independently [3, 75]. Thatsame year Fritz London and L. Tisza proposed that the explanation for this ’liquid degeneracy’
may be closely related to the condensation predicted by Bose and Einstein the previous decade[56, 123].
However, while interest in superfluid Helium and superfluidity in general remained highin the decades that followed, Bose-Einstein condensation was less fashionable. When Landaudeveloped his own two-fluid model in 1941, he did not consider BEC to be a promising explanationof Helium’s superfluidity [79]. Indeed the disagreement between Landau and Tisza was rather strong
as demonstrated by back-to-back letters to the editor of Physical Review in 1949 [80, 121]. Yet
the connection between superfluidity and Bose-Einstein condensation would continue to develop
via work of Bogoliobov which was published in 1947 and to this day remains an indispensable
framework for understanding elementary excitations of superfluids and Bose-Einstein condensates[17]. Bogoliobov’s transformation recasts the real, interacting, particles into a new basis of noninteracting quasiparticles. Diagonalization of the Hamiltonian in the quasiparticle basis provides
the dispersion relationship for weakly interacting superfluid systems. Underappreciated at the time
was the applicability of Bogoliobov’s treatment to Bose-Einstein condensation.
In the 1950’s work by Feynmann, Matsubara, Onsager, and Penrose revisited the subjectof superfluid helium and contended that London’s hypothesis was essentially correct; Bose condensation was fundamental to superfluidity. Onsager and Penrose calculated that even at zerotemperature, only 10% of the atoms were in the ground state due to depletion caused by the interatomic interactions [54, 96, 84]. Around this same time work in superconductivity was converging
on similar ideas, exemplified in the famous paper on superconductivity by Bardeen, Cooper, andSchrieffer in 1957 [11].

5
1.2.1.1

Spin-polarized hydrogen

In 1959, Hecht [64] considered systems of spin-polarized bosonic gases of hydrogen. and theprospect of spin-polarized and seemed puzzled that” “...it seems that no one has pointed out theexciting possibility that such atomic gases and liquids would show superfluid properties since H
and T atoms are bosons.” It seems that this observation was slightly ahead of its time; similar
observations were made by Stwalley and Nosanow in 1976 [115]. Experimental progress towards
stabling a sample of spin polarized hydrogen reached an important milestone in 1980 with the heroic
experimental efforts of Silvera and Walraven [111] to achieve such a stable sample using a dilution
refrigerator coated with a thin film a superfluid 4 He, which was required as otherwise the hydrogen
would instantly adsorb to any surface. However, progress towards BEC on that front slowed dueto the dueling requirements for high density and low temperatures and the enhanced losses suchdensities and temperatures caused through recombination processes. Attempts to reach BEC in
Hydrogen would continue, motivated in part by the determination that spin-polarized Hydrogen
would remain a gas at all temperatures, ensuring that a BEC would be robust against collapse[116]. Despite this fact, Hydrogen would have to wait a few years after other systems to finallybe condensed [48, 39]. In the efforts to cool Hydrogen, the technique of evaporative cooling wasdeveloped [83]. Evaporative cooling remains the method of choice for cooling bulk atomic gases toultracold temperatures and is the final stage of cooling used in for the work in this thesis.
1.2.1.2

Alkali atoms

While spin polarized hydrogen was being pursued, Hydrogen’s neighbors below it on the
periodic table became increasingly popular subjects of experiments on laser cooling and trapping,
despite the fact that the ground state of low temperature alkali systems is a solid. Beginningin the mid 1970s, it was proposed that lasers interacting with atoms could be used to slow andcool atoms [129, 63]. Tremendous progress was made in the next two decades in the experimental
development of techniques used to cool atoms such as Doppler and sub-Doppler cooling of neutral

6
atoms and ions [128, 100, 88, 34, 35, 117, 91, 8]. Ultimately these efforts led to the first observationof Bose Einstein condensation in a dilute atomic gas of 87 Rb here at JILA in 1995. This was thefirst really unambiguous demonstration of Bose condensation to demonstrate the phenomena as
Einstein had conceived it: in a gas, albeit non-ideal. These systems are not technically in theabsolute ground state, which is indeed a solid, but the extremely low densities slow the rate ofrecombination processes which nucleate solidification sufficiently to allow for long-lived metastable
BEC.
1.2.2

Interactions and early BEC

With the success of condensing 87 Rb here at JILA and 23 Na at MIT [7, 41], attention turned
to understanding these fascinating new systems, and the role of interactions was paramount. The
scattering lengths of the alkalis were not as well understood as it is now, although rapid progressin understanding interactions was swiftly made. For ultracold alkali atoms, the interactions can be
encapsulated almost entirely by a single parameter, the s-wave scattering length a [32]. It had been
predicted earlier that negative scattering lengths precluded stable condensation except in small
samples [106]. This was the reason for the difficulty in condensing 7 Li [21, 20, 31, 38]. Negative
scattering lengths also pose a problem for efficient evaporation generally. The rate of inelastic
three-body collisions is generally higher for negative a than for positive a of equal magnitude [19].
For these reasons positive s-wave scattering lengths are required for production of large BECs
by way of evaporative cooling. Early experiments demonstrated in spectacular fashion how a BEC
is incapable of persisting when the scattering length is tuned negative earning the name bosenova
for the similarity between condensate collapse and core collapse supernovae of very massive stars[105, 44]. These experiments produced a small, persistent remnant since for a sufficiently small
atom number the zero point energy stabilizes the sample against collapse.

7
1.2.2.1

General overview of the role of interactions from 2000 to present

Following the realization of BEC, interest was renewed in understanding these systems, particularly in the context of trapped gases where BEC was accessible. The behavior of BEC at
negative scattering lengths was explored in the context of one dimensional samples [97, 24, 23, 25].
These studies found that a 1D condensate remained stable at negative a and formed stable solitons.
Early experiments in one dimensional 7 Li samples produced such soliton trains following a quench
to attractive interactions [77, 114]. These experiments confirmed the prediction that during thecollapse, phase differences between the solitons determine their interactions; neighboring solitons
with a phase difference of π between them repel each other and prevent the local density from exceeding the critical density for total collapse. The origin of these phase fluctuations was attributed
to a combination of factors. Fluctuations of the initial condensate were proposed to seed these
relative phases in [2]. Alternatively, in [22] interference of the condensate order parameter during
the first stages of collapse were found to create soliton trains wherein the solitons behaved both
repulsively (with π relative phase) and attractively (with no relative phase) over the course of an
implosion process with multiple phases.
In three dimensions, such as the 85 Rb bosenova experiments [44] the picture was less clear.
The observations reported in [44] indicated that long after the quench to negative a, the numberof condensed atoms was well above the predicted critical number. This surprising result was explained by the fact that the implosion process was sufficiently violent enough to allow multiple
bright solitons to nucleate with relative phases such that the solitons interact repulsively despite
the attractive interatomic interactions [37]. This resembles the mechanism discussed above in 1D.
Subsequent numerical studies supported this notion [94, 95]. However, later numerical calculations
did not find a repulsive phase between the implosion remnants as originally postulated when fluctuations were included [42]. The question of how to interpret the results of experiments like the
JILA bosenova [44] has remained salient and continues to be pursued [14, 15, 120].
Solitons have also been generated in two dimensional system. In contrast to 1D and 3D, the8two dimensional solitons were stabilized by rapidly oscillating the scattering length [108, 33, 61, 31].
However, our interest in this thesis is not in the nature of solitons in various dimensions, but inunderstanding fluctuations in bulk 3D BEC.
1.2.3

Fluctuations

The nature of fluctuations in Bosonic systems of varying dimension is the subject of much of
Chapter 4, where we shall derive quantitative results for uniform and harmonically trapped systems
for d = 1, 2, 3. In lower dimensional systems, the fluctuations are larger, to the point that theypreclude the existence of a true, phase-coherent BEC. This fact has been understood for many
decades. Mermin and Wagner [86] provided a a proof that long-range order is impossible at any
finite temperature in a one- or two- dimensional isotropic spin-S Heisenberg model. A consequenceof this fact is that Bose-Einstein condensation is impossible in such systems [85, 67]. This follows
from the fact that the systems considered by Mermin and Wagner belong to the same universality
class as Bose gases [68].
In 3D, the phase fluctuations (in a uniform gas) for large distances are bounded from above
by a constant, as was derived for finite temperature by Kane and Kadanoff [74]. The one-bodydensity matrix therefore asymptotically approaches a finite value. In one and two dimensions, themagnitude of the phase fluctuations is unbounded for large distances, and therefore the off-diagonal
portion of the one-body density matrix decays to zero. Interestingly, this decay is exponential for
one-dimension but algebraic for two-dimensions. This turns out to be related to the BerezinskiiKosterlitz-Thouless (BKT) transition [78, 13]. We will not be discussing the BKT transition in thisthesis.
The fact that the phase fluctuations remain large over large distances in 1D and 2D mathematically follows from the nature of the ∝ p−1 momentum dependence of the phase operator
for the condensate order parameter and the convergence, or lack thereof, of a momentum integral (Eq. 4.130). The analogous calculation for density fluctuations always converges (Eq. 4.135),
because the momentum dependence ∝ p of the density operator is linear. Hence, the density fluc-9
tuation effects will tend to be smaller than those of the phase fluctuations. Fluctuations in lower
dimensions are larger than in 3D as a consequence of the volume element of the integral.

1.3

A brief overview of my career
Physics has always been a passion of mine, and ultracold atoms in particular have always

been a topic of fascination. I was fortunate enough to work with laser cooling and trapping ofatoms as an undergraduate at Lawrence University, and continued in this direction when I arrived
at JILA, originally in the bi-group of Debbie Jin and Eric Cornell. My first few years at JILA werespent working with Bose-Fermi mixtures of 87 Rb and 40 K and studying their interactions, which
nature has allowed us to tune simply by changing the magnetic field the atoms experience. We also
worked to create very smooth optical potentials amenable to rapid rotation of these ultracold gases.
I owe a great deal to my fellow graduate student Ming-Guang Hu, who was an excellent teacher andenjoyable to work with. Ming-Guang and I were joined by Dhruv Kedar in 2015. While the grandvision of the lab in B231 was aimed toward rotating Bose-Fermi mixtures with the ultimate goal ofinvestigating quantum Hall physics in ultracold atoms, we took a detour to perform experiments
on Bose polarons in the strongly interacting regime [69]. This work was tragically interrupted by
Debbie’s illness and death in 2016. The world was robbed of a brilliant, kind scientist, and I wasrobbed of a wonderful advisor. However, I’m not sure how the aftermath of this tragedy could havegone better for my fellow graduate students and I; our voices were heard, listened, and acted on byprofessors and advisors who clearly cared deeply about helping Debbie’s students succeed. Dhruv
joined Jun Ye’s experiments on Strontium optical lattice clocks, and I was sad to cease working
with him, he was a wonderful labmate.
Ultimately, I found myself switching to slightly different, but very related work in what
was previously the Fermi Gas lab in Debbie’s group. I learned to operate my second BEC creation
machine of graduate school, and with a new species of atom, 39 K, which would allow us to investigate
and expand the limits of our understanding of few-body physics. I worked with my fellow graduate
students Roman Chapurin and Xin Xie, as well as with undergraduates Jerod Popowski, Carlos

10
Lopez-Abadia, and Bjorn Sumner. Things proceeded well in the new lab and we were able toexperimentally and theoretically placed Efimov physics on a more robust universal foundation than
had yet been realized. Noah Schlossberger joined our lab, but in what felt like no time at all, heswitched to the eEDM Gen-III experiment. From late 2016 to early 2020, the JILA Resonantly
Interacting Boson lab (JRIB) focused on few-body physics. We produced detailed measurementsof the Feshbach resonance in the |1, −1i state of 39 K and its associated Efimov spectrum, providing
a through picture of this effect and testing the limits of universality [28, 132].
In late 2019, an collection of chronic medical conditions I have lived with since childhood
coalesced into an obstacle I could no longer ignore and I required some time off and a restructuringof my schedule to accommodate my treatment. This was compounded by the global Covid-19
pandemic and the most pathological series of objective lenses that myself and my advisors haveever seen. By the middle of 2020, Xin Xie and Roman had both graduated and I was the solemaster of my laboratory, which certainly has its pro and cons. A good deal of time and effort was
expended in the analysis and trouble shooting of the high resolution imaging, a problem which was
not satisfactorily solved until a replacement lens arrived in fall 2021 which reduced aberrations to
a tolerable level. From there, it was a short two months until data collection was completed. The
road to results has been long and difficult, and the path taken has diverted somewhat from theoriginal vision for an fluctuation microscope. To understand our successes and struggles in thisproject, a prologue of sorts for the original scientific vision is useful.

1.4

Vision for a fluctuation microscope
Fluctuations in BEC are important to understanding systems of ultracold bosons, particularly

when the interactions between the atoms are strong. Unfortunately, they are also very difficult toobserve in 3D. The goal of this work is to overcome the difficulties of accurately characterizing subtle
density perturbations that are the signatures of these fluctuations by amplifying those perturbations
via a quench to attractive interactions. The general concept is straightforward: an initially stabletrapped condensate undergoes a rapid interaction quench to negative as which suddenly changes

11

(a)

(b)

Figure 1.1: Illustration of the evolution of a local density perturbation in a uniform gas. A initial
density bump dissipates for repulsive interatomic interactions as shown in (a). A initial density
bump grows in magnitude and contracts in size following a quench to strong attractive interactions
as in (b). The curves are for qualitative and illustrative purposes only.

the mean field potential experience by the atoms from repulsive to attractive. The inward force on
the cloud increases as the cloud contracts, and local density maxima continue to drive this inward
movement. The densest regions of the cloud pull harder and harder on all the atoms in the vicinity.
Consider an infinite uniform BEC sitting happily at a stable positive scattering length andsuppose there exists a single isolated density bump, perhaps the results of a externally imposed
local potential well which is switch off at t = 0. The aligned forces of kinetic energy pressure
and repulsion due to as > 0 each push particles away from higher density regions. The density
profile will thus tend to dissipate to match the background density as show in Fig. 1.1a. If theinteractions between the particles was instead quenched at t = 0 to a sufficiently attractive value,
the evolution of the density bump takes a very different turn. An atom at the edge of the bump will
still be pushed outwards due to the kinetic pressure of the wavefunction, but for a sufficiently strong
interatomic attraction, the net force on the particle will be inward and the bump will grow until
other processes such as molecular recombination become relevant. The growth of the perturbation
following such a quench is show in Fig. 1.1b.
Of course, even if the wavefunction is perfectly uniform at t = 0, any measurement of the12

Figure 1.2: Illustration of how density fluctuations are amplified in a uniform gas, the data are
randomly generated and for illustrative purposes only. A measurement of the density profile att = 0 (blue line) would yields shot noise on top of the uniform density. If the density profile were
measured at a later time (orange line) following a quench to attractive interactions, the density
perturbations from shot noise will have concentrated in the dense regions, thereby being amplified
compared to there initial amplitude. This is illustrated by the orange line magnifying the preexisting
structure indicated by the blue line.

system will inevitably contain shot noise. If the atomic density is measured at t = 0, the result
would more closely resemble the t = 0 curve in Fig. 1.2. However, density shot noise at t = 0 may
be small and difficult to measure. Consequently, if the condensate is allowed to evolve for sometime after the quench, those same fluctuations will be amplified enough to be measurable.
Above T = 0 K the cloud will also have fluctuations due to the temperature, and interactions
cause additional quantum fluctuations apart from the ubiquitous shot noise. The density perturbations are simply the coordinate space signature of the elementary excitations of the Bose fluid.
In momentum space, these excitations are simply plane waves traveling in all different directions.
In the case of a trapped BEC, thermal fluctuations, quantum fluctuations, and shot noise are all
still present, and a quench to negative scattering length will still serve to amplify local densities
perturbations if the density is sufficiently high or the scattering length sufficiently large. Whilein the uniform gas there is no external potential such that the evolution of the wave function is13
governed solely by kinetic and interaction dynamics, in a trapped gas there exists an additional
inward force. The density also drops to zero at the edge of the cloud, so the force due to theinteraction potential is exclusively inward at the clouds edge. In essence, the cloud will implodeglobally following a quench to negative a in addition to the local dynamics.
1.4.1

Wavefunction implosion

While the logic that a sudden change of interaction from repulsive to attractive would leadto the BEC imploding is intuitively reasonable, we can be more quantitative here in elucidating
the physics behind implosions without delving yet as deeply into the theory as we shall in Chapter
4. Yet a few equations are not unwarranted here. The behavior of BECs is commonly consideredusing the time-dependent Gross-Piteavskii equation (GPE):


~2 2
∂
2
−
∇ + Vext (r, t) + g|Ψ(r, t)| Ψ(r, t) = i~ Ψ(r, t).
2m
∂t

(1.1)

The GPE, is a non-linear Schrodinger equation incorporating interactions between the particles[102]. It provides a good model to expand upon the discussion of 1.4 by analyzing each of its terms.
~
The first term − 2m
∇2 Ψ represents what I have been calling the kinetic pressure of the gas. Where
2

the wavefunction has large spatial curvature this term is large, corresponding to a large wave vector
k and hence large momentum ~k. Such a wavefunction will not generally remain localized in theabsence of other confining factors, as things with large kinetic energy tend to move. Hence, a very
tightly contained Ψ will spread out, as though responding to an internal pressure, and this tendency
will be stronger the more tightly contained and thus rapidly varying Ψ is.
The second term on the LHS of Eq. 1.1 is the external potential. In the context of my work
trapping atoms, this is most often a trapping potential, provided in the lab by magnetic or opticalfields, but is not always. Gravity is a non-trapping potential contributing to Vext (non-trapping in
the lab case, certainly trapping on the scale of the planet). In our experiments, Vext is nearly always
modeled as a harmonic potential with its minimum defining the origin of our coordinate space. The
optical lattice potentials which we impose to seed density perturbations are an exception to this

14
rule.
The third term is the star of the GPE show; a flamboyant new character to Schrodinger’s
cast, ready to cause drama and incapable of minding its own business. This is because the third
term represents the interactions between atoms. We can rewrite the interaction term as g|Ψ|2 ∝ anwhere the density is of course n(r) = |Ψ|2 and the interaction parameter g ∝ a is linear in the
scattering length. Whether this term aids the kinetic pressure in pushing the atoms apart or theexternal potential in forcing them together depends on the scattering length a, which is among
the most important quantities throughout this work. When a is negative, the interaction term gn
produces a confining potential, whereas a positive a provides a repulsive potential. This term isalso called the mean field (MF).
With this in mind, we consider the balance of these forces for a trapped Bose cloud, and inparticular consider how the kinetic energy, confinement potential, and mean-field interactions willchange when the cloud departs from its equilibrium size. Before the quench, if the cloud is toolarge, the kinetic energy and mean-field repulsion each lessen while the trapping potential becomessteeper; the end result being the cloud contracts. If the cloud is too small, the opposite happens
and the cloud will expand back towards equilibrium. This is the balance of forces at work in Fig.
1.3 on the left.
When the scattering length is suddenly quenched to a negative value, the balance of forces may
no longer be stable. The wavefunction Ψ is unchanged during the quench if it is fast enough, andthe trapping potential remains the same, so they act as they did before the quench. The interaction
potential is a very different matter and the cloud will begin to contract as the outward pressure of
the mean-field is suddenly replaced with an inward pull. But in contrast to the repulsive case, thisproduces positive feedback since the mean-field potential scales with the density, and as the cloud
implodes it becomes denser and denser. For weak enough interactions, the negative feedback of theincreasing influence of the kinetic term can eventually halt this collapse and a smaller equilibrium
cloud may become stable, but for dense clouds or strong interactions, the positive feedback of the
imploding cloud’s mean-field attraction outpaces the negative feedback response of the kinetic and15

Figure 1.3: The balance of forces changes with a quench to attractive interactions. Before thequench all three components provide negative feedback to perturbations away from equilibrium.
After the quench, the mean-field potential instead provides positive feedback

potential terms in the GPE. To be clear the potential Vext is always inwards, but weakens as thecloud shrinks, providing the negative feedback.
1.4.2

Mode stability

While the GPE is certainly a quantum mechanical expression, the argument and physicalintuition provided by the previous Section feel quite classical in nature to me personally. A complimentary picture for understanding implosions is provided by the theory of elementary excitations,first developed by Bogoliobov [17]. The details of Bogoliobov’s theory are explained in more detailin Chapter 4, but the results of important here is the dispersion relation for the excitations upona Bose-condensed gas:"gn 2
(k) =
k +
m



k2
2m

2 #1/2(1.2)

This dispersion relation is derived in Chapter 4 starting from assumptions of weak interactions
in a uniform condensate at positive scattering length, but remains reasonably accurate even forsystems without these attributes [62, 72] For positive a, 2 (k) > 0 for all densities, and the k-modesk
may be grouped into two regimes: the high-k particle-like excitations where the kinetic energy ~2m

2 2

dominates over the interaction energy, and the low-k phonon-like excitations. For a > 0 the modes

16

Figure 1.4: Dispersion relation for positive and negative a. On the positive side the low-k modes
exhibit linear momentum dependence like those of phonons while at large k the modes exhibit thefree-particle-like quadratic behavior.

are all stable. For a < 0 the situation is very different. As mentioned, the same dispersion relation
applies even for a < 0, and this fact has important consequences for how these modes evolve intime. A quantum state ψ will evolve according toEψ(t) = ψ0 e−i ~ t .

(1.3)

k
Expressing the energy in terms of frequency E = ~ω and considering Eq. 1.2 when |g|n > ~2m

2 2

we see how the frequency associated with (k) is purely imaginary for a range of wavenumbers.
The effect on the evolution of ψ according to Eq. 1.3 is to induce the exponential growth of thesemodes. Now the connection to density fluctuations and how implosions should amplify them is
becoming apparent. Density fluctuations are equivalently expressible in terms of superpositions of
momentum modes, and those modes, or rather a subset of them, will grow following a quench toa < 0; and in a way which depends on the post-quench scattering length, the cloud density, andthe wavenumber of the mode. The growth of the modes goes hand in hand with the amplification

17of the associated density perturbations.
This is the motivation for our use of an interaction quench to amplify the subtle fluctuationsof our BEC such that they may be directly observed. The discussion has been necessarily limited
and incomplete, and indeed as our experiments were performed and analyzed, our understanding
improved, and the simple physical picture presented in these Sections became more complicated.
1.4.3

Calibrating fluctuations measurement

The fluctuations we seek to amplify are detected via in-situ absorption images of condensates
following a quench and subsequent evolution of the condensate where the fluctuations grow during
the implosion process outlined above. The in situ density distribution is then analyzed to computethe structure factor S(k) defined as:
S(k) =

1h|δn(k)|2 i.
N

(1.4)

The structure factor is simply the average squared value of the Fourier transformed density for aparticular k. Experimentally, this is determined by taking an ensemble of in-situ density measurements, subtracting the mean density profile from each individual realization, and computing the
FFT of the remainder. It is therefore critical to understand how structure in the clouds density is
captured by the imaging system. Previous work in this lab was essentially always concerned withbulk properties of the cloud, and in-situ imaging was unnecessary for those projects. In contrast,
understanding and calibrating the newly installed in-situ imaging is essential to the business ofmeasuring fluctuations. The design and characterization on bench tests was carried out by Xin
Xie, and further information may be found in her thesis [131], but Chapter 3 of this work is dedicated to our continuation of her thorough work. Generically, the task of in-situ measurement of ahigh optical depth BEC is a difficult one, and we understood that to have confidence in our results
demanded an empirical method of calibrating the transfer function. This is the basis for our latticecalibration.
The calibration lattices are the other important addition to the apparatus in this work (the

18

(a)

(b)

Figure 1.5: Real space and Fourier space density distributions. The red circle in (b) corresponds
to the Y-lattice, with a lattice spacing of 2.5 µm while the green circle marks the X-lattice peak in
Fourier space. In (b) the central peak corresponding to the culk distribution of the cloud is masked.

first being the high-resolution imaging just mentioned). With optical lattices, we have the ability
to write density fluctuations directly onto the cloud and know with high confidence precisely how
deep or shallow the lattice is. The lattice potential can be characterized via establish methods
using short lattice pulse to excited momentum modes in units of 2~kL which are then measured
in time-of-flight. The transfer function of the imaging can then be calibrated, and by decreasing
the depth of the lattice to a very small fraction of the chemical potential, we may calibrate oursensitivity to structure at the wavelength of the lattice. For this reason two lattices with differing
periodicity are used in this work.
This turned out to be very important for a responsible methodology, as the imaging system
did not perform as expected for several reasons, some of which are still not fully understood.
Hints of this are already present in Figure 1.5, where the signal for the Y-lattice seems to be
distorted and very near the edge of the resolved region, despite bench measurements and the designindicating this spatial frequency should be well within our resolution. When the terms X-latticeor Y-lattice are used, they refer to the direction of that lattice’s wave vector. The details of the19lattice characterization are presented in Chapter 5

1.5

Contents of this Thesis
The rest of this thesis will build upon the vision present in Section 1.4. In Chapter 2 I

review the apparatus, focusing on those components which are relevant to the production of BECs
but without going into the level of detail such an apparatus chapter often does because I am thelast person to use this apparatus. Those components of the apparatus which are critical for thescience of this thesis have been given their own Chapters. These include the optical lattices andthe high-resolution imaging system.
Chapter 3 will contain the more detailed description of the imaging systems as it represents amajor component of the experiment that is mostly new to this work. The principles of absorption
imaging will be discussed in detail. Certain features of data analysis are intimately connectedwith the imaging procedure and as such will be reviewed in this chapter as well. These include
the expected density profiles, principal component analysis, and thermometry. The Fourier space
analysis of the imaging will be deferred to later Chapters.
Chapter 4 will review the relevant theory of BEC in general and fluctuations in particular.
The ideas present in Section 1.4 will be discussed in greater detail. The structure factor and its
time evolution following a quench will be expanded upon and the arguments of Section 1.4 will beboth supported and challenged as a result of a more thorough theoretical analysis.
Chapter 5 will discuss the control lattices. The calibration of the lattice will be describedin detail and compared with a priori expectations, with tension between expectations and measurements being reviewed. Predictions of the theory given in Chapter 4 will be compared to ourcontrolled implosions with well characterized lattices. The effects of lattice potentials persisting
during the implosion time are found to be in agreement with theory, and an unforeseen applicationof implosions to characterizing our confining potential will be explored.
Chapter 6 will review our experiments of implosions in the absence of a lattice potential. The
predictions of Chapter 4 will be applied to our efforts using implosions to amplify fluctuations. The

20
successes and failures of our experiments to measure the static structure factor will be evaluated andthe final results summarized. We will return in some detail to the question of how to analyze thedata in Fourier space in this Chapter as well as part of the discussion on experimental challenges.
Chapter 7 will contain the conclusion and outlook of this thesis. Additionally, several appendices containing results which are useful and necessary but which do not warrant placement in themain text are included.

Chapter 2
Apparatus

The work in this dissertation departs from the few-body focused studies of Feshbach andEfimov physics led by my fellow graduate students Roman Chapurin and Xin Xie. However, withthe exception of the vertical imaging system and optical lattice, the apparatus and cooling procedure
follow those outlined in their theses with very little modification. Given my position as the lastperson to use this lab, this apparatus chapter will never be used as a reference in the way that
Roman’s and Xin’s theses were used by me. As such, I will not go into nearly the same level ofdetail as can be found in their theses and will be referring the reader to their theses for a substantial
portion of this chapter [27, 131].

2.1

Overview of cooling procedure
While our atoms eventually end up at ultracold temperatures well below one µK, they start

above room temperature, and over half each experimental cycle consists in cooling the atoms to theultracold regime. Our apparatus has three vacuum chambers: MOT1, MOT2, and the science cell.
The cooling process begins with two three-dimensional magneto-optical traps (MOTs). The energy
level diagram of the three manifolds relevant to laser cooling of 39 K is shown in Fig. 2.1, and the
Zeeman spectra for the S1/2 , P3/2 , and P1/2 are shown in Figs. 2.2, 2.3, and 2.4 respectively. The
collection chamber (MOT1) cools the atoms from well above room temperature down many ordersof magnitude, perhaps down to about 1 ∼ 10 mK; the temperature in the collection MOT is not
something we have measured since I joined the lab in late 2016. Cooling is a secondary purpose

22

Figure 2.1: Energy level diagram for 39 K. The detuning of each hyperfine sublevel with respect tothe fine structure energy is given in parentheses in units of MHz. The dotted line indicates the finestructure energy in absence of hyperfine interactions. The g-factors for each level are computed
analytically from the quantum numbers as in [87]. All calculations for this figure and for Figs. 2.2,
2.4, and 2.3 were performed by myself in Matlab using reference data for 39 K found in [122] orusing NIST CODATA values.

Figure 2.2: Ground state manifold. This manifold is the lower line for both the D1 and D2
transitions. Repump D2 (D1) lasers couple the F = 1 states to the P3/2 (P1/2 ) manifold while trap
D2 (D1) lasers couple the F = 2 states to the P3/2 (P1/2 ) manifold.

23

Figure 2.3: P3/2 manifold. This manifold is the upper line of the D2 transition.

Figure 2.4: P1/2 manifold. This manifold is the upper line of the D1 transition.

24of MOT1 compared to the more vital task of collecting as many atoms as possible. The MOT1chamber possesses a getter which provides potassium vapor to the chamber which is heated toprevent adsorption to the walls thereby allowing reasonable collection rates. During initial loading
the atoms are continuously launched from MOT1 through a narrow tube to MOT2 by a push beam
propagating along the tube and are guided also by a hexapole field consisting of six permanent
magnets on the outside of the tube.
The narrow tube also serves as a narrow enough restriction that the pressure in MOT1 issignificantly higher than in MOT2, allowing maximum atom collection in MOT1 and reducing
background pressure in MOT2 which might hinder cooling. Unlike MOT1, MOT2 is not loaded
from the ambient vapor pressure but from MOT1 by way of a beam of atoms launched from MOT1
to MOT2; therefore the alignment of the push beam from MOT1 to MOT2 is critical. The D2
trap and repump light are derived from the same lasers for both MOT1 and MOT2. The Zeeman
spectrum for the ground state manifold is shown in Fig. 2.2 and that of the excited state of theD2 transition is shown in Fig. 2.3.
For the first 20 seconds of the experimental cycle, MOT1 is on at full power, pushing as
many atoms as possible to MOT2; loading beyond this time yields diminishing returns. Following
the MOT1 to MOT2 loading, a 20 ms endMOT1 stage is performed by detuning the D2 beams
further from resonance, providing transitional cooling to the next stage; the endMOT is similar tothe compressed MOT stage described next, but using only D2 light. Following the endMOT, theMOT1 light is turned off and MOT2 undergoes a 2 ms compressed MOT (CMOT) stage where
the D2 trap light is turned off while the D2 repump light is jumped to a new frequency. The D1
trap and repump beams are turned on at this time and the magnetic field is increased as well (thecompression part). Following the CMOT stage comes the grey molasses stage (GM) where the D2
repump light is turned off leaving only the D1 light on the atoms. This stage lasts about 6 ms andcools the atoms down to 10µK. Further details of the frequencies, gradients, and timings may befound in previous theses [27, 131].
1

The term endMOT is inherited from the graduate students before me.

25
The GM stage ultimately brings the temperature of the atoms in the MOT2 chamber to
10µK. The atoms are then shuttled via a mechanically translated quadrupole trap (QT). The
MOT2 primary coils are mounted on a cart for this purpose and are also called the cart coils. The
loading into the cart QT raises the temperature to 100µK. When the cart arrives at the science
cell, current in the cart coils is ramped down while that in the science coils ramps up to form thescience QT. When the cart QT has been ramped off, the cart is returned to the MOT2 chamber.
The optical dipole trap (ODT) is then loaded for about six seconds before the QT is turned off,the science coils are switched from anti-Helmholtz to Helmholtz configuration, and a bias field
is applied. The bias field is necessary to bring the scattering length for the |1, −1i atoms to asuitable positive value. The all-optical evaporation procedure then takes place, taking several moreseconds to bring the atoms to the desired temperature, typically < 300µK. During evaporation, theatoms are transferred from the initial loading ODT beam to a crossed ODT formed by a vertically
propagating beam and a horizontally propagating elliptical beam. These two beams provide aslightly oblate harmonic trap. Details of the ODT will be covered in a later section.

2.2

Lasers for cooling
The diode laser system is described in great detail in Roman’s thesis [27]. With the exceptionof some rearranging of a few imaging optics to provide a slightly more convenient alignment andpower balancing procedure, the optical layout is identical to Fig. 2.15 of [27], with the only
difference being the addition of a few mirrors and ND filters immediately before the top and sideprobe AOMs at the top of said figure.

2.3

Magnetic field control and trapping
The MOT1 primary and shim coils are driven and servoed by a rack-mounted system which,in contrast to every other coil in the experiment, do not have their control voltage set by thecomputer, but rather by manual potentiometers. The MOT2 shim coils are a different story,
indeed these coils take different values for the MOT1 to MOT2 collection, endMOT, CMOT, GM,26

Figure 2.5: Magnetic field control enclosure. I constructed this enclosure as part of several efforts in
our lab to reduce magnetic field noise. The box contains the H-bridge relays for the science coils in
the rear of the box (upper left, difficult to see). It also contains the following for the cart coil circuitand the science coil circuit: current sensors, IGBTs, cooling plates for the IGBTs, servo electronics,
sense resistors and their TECs, and inputs/outputs for water, analog and digital control, and finally
power.

and QT load stages. The optimum shim current values were optimized every few months, mostoften using the atom number and cloud size in the science QT as the diagnostic observables.
2.3.1

High current coils

The high current coils mounted on the cart and science cell are driven by two large HP
power supplies capable of outputting 400 A. The science coils circuit contains an H-bridge forswitching between Helmholtz and anti-Helmholtz configurations. Fluxgate transducers paired withhigh stability sense resistors mounted on actively stabilized TECs provide in-loop and out-of-loop
current measurements. The servo outputs for the science coils and the cart coils each connect to theinsulated gate bipolar transistor’s (IGBT) gate controlling the coils current. The science coils may

27
be operated in Helmholtz or anti-Helmholtz configurations, and because the required current is 10
times larger when the anti-Helmholtz configuration is used, two separate servo systems completewith current monitors and servo electronics, are used. Which servo controls the science IGBT isdetermined via a digital switch. The current monitors, sense resistors, sense resistor TECS, servos,
H-bridge, IGBTs, and water-cooled IBGT mounting plates are contained in a large aluminum box I
constructed in early 2019, consolidating several free standing components into a shielded enclosure,shown in Fig. 2.5. Combined with other mitigation efforts like moving sources of B-field noise awayfrom the atoms, and improved filtering of the servos, we were able to dramatically improve ourfield stability, ultimately reaching about 2-5 mG RMS shot-to-shot variation; see previous theses
for details [131, 27].
2.3.2

Low current science coils

The science cell also contains several smaller coils used for a variety of purposes. Two pairsof shim coils, one parallel to the horizontal imaging axis which is used to properly bias the field
for low field σ − imaging, and the other perpendicular to that axis which was occasionally used
to provide a slosh-inducing gradient on the atoms. Vertically, a single coil coaxial with the large
bias coils provides a gradient suitable for Stern-Gerlach measurements. A final, specialized coilfor rapid B-field quenches, called the fast-B coil, is also positioned coaxial to the bias coils. In
the early vision for this experiment, the quench to attractive interactions was to be achieved bya very rapid jump of the magnetic field. The fast-B coil requires a much more finely tuned andsophisticated servo including a feed-forward capability to ultimately allow ramps of over 10 G in
5-10µs. Ultimately, the quench to negative scattering length was instead achieved via state transfer
to the |2, −2i weakly interacting state. In this thesis, ket notation |F, mF i refer to the hyperfine
angular momentum quantum numbers. More details on the fast-B apparatus and protocols can befound in Xin Xie’s thesis [131].

28

2.4

State Control
Our atoms are prepared in the |1, −1i sublevel of the 42 S1/2 ground state of 39 K; see Fig.

2.2. This is a low-field seeking state, as is required for confinement and transport in the QTs. Prior
to loading the cart QT in the MOT2 chamber, the repump D1 light is extinguished while the D1
trap light remains on and rapidly pumps the atoms out the of F = 2 manifold. Only the mF = −1
sublevel is loaded into the cart QT and therefore only |1, −1i atoms are loaded into the ODT. Spinpurity is required to sustain efficient evaporation. The |1, −1i state’s Feshbach resonance near 32.6
G provides the tunability of the collisional cross section needed for evaporation, as seen in Fig.
2.6a. The other spin states either posses a negative scattering length or else have a substantial
imaginary component as shown in Fig. 2.6b.

(a) Real part of a(b) Imaginary part of a

Figure 2.6: Real and imaginary parts of the intra-species scattering length. The atoms are preparedin the |1, −1i state which allows tuning of the scattering length via the Feshbach resonance near
32.6 G. The imaginary part of the scattering length is related to two-body inelastic collisions.
Implosions are performed by flipping the atoms into the |2, −2i state, which is also the imaging
state. The |2, −2i state scattering length is flat over the depicted B-field range at 33.4 a0 .

29

2.5

Hyperfine transitions
Previous work in our lab used a wider variety of spin states than was required for this work.

These transitions are driven by a collection of antenna tuned to the transition frequency for a given
pair of states and the magnetic field. Details regarding the tuning, construction, and mountingof these antenna may be found in Roman and Xin’s theses [27, 131]. The transition between spinstates may be driven with a resonant microwave or RF field or via adiabatic rapid passage (ARP),
also known as a Landau-Zener transition [135]. The most important transition for this thesis is
the |1, −1i → |2, −2i microwave transition. The |1, −1i state is the preparation state used in
evaporation and is the initial state of the atoms priori to the quench to negative interactions. The|2, −2i state is both the state to which we quench as it possesses a small negative scattering length(see Fig. 2.6a); it is also the lower state of the cycling |2, −2i → |3, −3i transition used for imaging
the cloud. When speed is not an important factor, such as for imaging expanded and therefore
low density clouds in TOF, we use an ARP via a sweep of the frequency across the resonance.
When the initial and final frequencies are sufficiently detuned from resonance, the atoms will follow
the dressed states which asymptotically approach the bare states away from resonance. ARPs are
more robust and stable against drifting B-fields as the initial and final detuning can be very large
compared to any the change in resonant frequency due to changing B-fields. While ARPs arehighly reliable and only rarely require optimization, they also take a considerable amount of timeas the frequency must be swept slowly with respect to the resonant Rabi frequency. Our ARPs
require about 200 µs to complete, too long for the quench experiments in this thesis. Ideally, ourquench from positive to negative scattering length to amplify implosions would be instantaneous.
An instant quench is obviously impossible, but the slow timescale of the ARP may be dramatically
improved upon by instead using a resonant π pulse at the cost of increased sensitivity to inevitable
magnetic field noise. For the |1, −1i → |2, −2i transition used to initiate implosions, the pulse
lasts 6-8 µs, sufficiently rapid with respect to the implosion dynamics to be negligible. Examplesof resonance position and Rabi-frequency measurements are shown in Fig. 2.7a.

30

(a)

(b) Two measurements of duration of a π pulse.

Figure 2.7: Measurements of the center frequency and pulse duration for the |1, −1i → |2, −2i transition. In (a) are two measurements of the resonance position for the |1, −1i → |2, −2i transition.
In (b) are two measurements of duration of a 2π pulse. The 2π pulse yields the minimum observed
optical depth, this time is divided by two to give the duration of a π pulse. The residual signal inthe π pulse minima in (b) also appears at t = 0 microseconds and is related to the extremely highoptical depth of the in situ cloud.

2.6

Optical Dipole trap
The optical trapping fields are provided by three high powered beam lines focused within

the science cell. The ODT light is derived from a 50 W amplifier made by Nufern which is seededby low power (15-50 mW) 1064 nm laser. This laser is divided into three beams (H1, V, and H2)each coupled to a photonic crystal fiber (PCF). The optical layout of the 1064 nm system up to
and including the three PCF input couplings is enclosed in black metal paneling for safety, andis unchanged from previous theses [27, 131]. The PCF fibers route the H1, V, and H2 to theirrespective breadboards and each beam is ultimately aligned to the center of the science cell. A
photo of the area surrounding the science cell is given in figure 2.8. The two horizontal trapping
beams, H1 and H2, are mounted on the side breadboard shown in the schematic of Fig.2.9. TheODT loading beam, also referred to as H1, is focused to a waist of 28 µm and precisely aligned
to a favorable position within the science cell just off center of the QT. If H1 is too close to thezero of the QT, Majorana spin flips transfer many atoms into states with poor collisional behavior

31

Figure 2.8: A view of the area surrounding the science cell. This image is intended to give thereader an idea of the space constraints of the experiment; optical layouts are given in Figs. 2.9,
2.10, 2.11, and 2.14 which will provide much more clarity and understanding. The copper bias coilsof the cell are visible just to the right of center. The H1 and H2 beams are enter the science cellfrom the breadboard in the bottom-left portion, see Fig. 2.9. The breadboard in the top half of
the image holds the V beam input optics and many of the optics of the X,Y, and Z lattice, see Fig.
2.10. The periscope diverting light from the H1 input optics to the top breadboard may be seen at
the far left, where a 1” mirror diverts the beam upward and a 2” mirror mounted on a 1.5” steel
post steers light onto the top breadboard.

and essentially zero atoms are loaded. The atoms are transferred from H1 to the vertical (V), andsheet (H2) ODT beams during the evaporation procedure. Atom number and temperature at theend of evaporation is sensitive to the relative alignment of these three beams. All three beams are
independently intensity stabilized with a JILA AOM intensity servo and in-loop photodetector. An
additional out-of-loop photodetector allows monitoring of the beams’ power for each beam.

32

Figure 2.9: Layout of the H1 and H2 beam optics. The breadboard for the input optics (up-beamof the science cell) is mounted such that the beam high is level with the atoms and is depicted ingrey. Optics on the input breadboard are positioned roughly to scale. The H1 beam is used only
during loading of the ODT from the science QT and is focused to a 28 micron waist. The H1 beam
propagates at a 45 degree angle to the X-axis. During loading the H1 beam is ramped off as theH2 beam is ramped on. The H2 beam is focused in the vertical direction but not the horizontal
direction by the +500 mm cylindrical mirror mounted on a longitudinal translation stage to alignthe waist to the atoms. The figure shows the behavior of the Z-size of the beam from the +500
cyl lens to the collimating lens after the science cell. The H2 beam propogates along the X-axis inthe −x̂. For both H1 and H2 (and V, see Fig. 2.11) the beams are collimated and sent to servo,
monitoring, and alignment optics, although most light is dumped into water-cooled beam dumps.
The optics down-beam of the atoms are not shown to scale

33

Figure 2.10: Top breadboard. This breadboard holds the input optics for the V beam and thevertical/top probe light. It also receives the lattice light from the H1/H2 breadboard via a periscopefrom the H1 optics. The lattice light is shifted by an AOM to offset the X, Y, and Z lattice beamfrequencies. The Y and Z lattice beam lines are fiber coupled and routed to the Y and Z lattice
optics of Fig. 2.14. The X lattice light is delivered alongside the V beam light and directed down
by a short-wave pass dichroic mirror through which the top probe light is transmitted. The ODT
V beam and the top probe beam are overlapped with each other down through the science cell andatoms until separated by a long-wave pass dichroic below the science cell as shown in Fig. 2.11.

34

Figure 2.11: Top probe light and V ODT beams. The inset box shows the downward propagationof the probe and V beam through the science cell, outside the inset the beams are propagating
horizontally. All beam sizes are exaggerated for clarity, particularly the probe light which is light
blue-grey and is collimated at the atoms. To compensate for the +125 mm lens which focuses the
V beam and steers the X-lattice beams (not shown), the probe light is not collimated when it iscombined with the V beam by the SWP dichroic depicted in green at the top of the inset. The
probe light is resonant with the |2, −2i → |3, −3i cycling transition and is circularly polarized fora σ − transition. Light scattering by the atoms is depicted by the blue dashed lines and is collectedby the +30 mm objective. The unscattered light is divergent after the objective and due the theexaggerated size of the probe beam at the atoms is depicted as diverging dramatically after theobjective. The unscattered light is much less divergent than indicated. The bias field is produced
by Helmholtz coils above and below the science cell which are not show and the bias field oriented
in the ẑ direction. The V beam of the ODT is shown in red. Below the science cell and the objective
the probe light is reflected by the LWP dichroic mirror (red) to the eyepiece and camera. The V
beam is transmitted through the dichroic and sent to the servo, monitoring, and alignment optics.
After the LP-dichroic path lengths are not to scale and several mirrors are omitted.

35
2.6.1

Trap characterization

Accurate knowledge of the trapping potential is critical in determining the density of the
BEC. In fact, the properties of the condensate relevant to this thesis are determined from six
basic quantities: scattering length, total atom number, temperature, and the trapping frequency
in the x̂, ŷ, and ẑ directions. Temperature measurements typically require knowledge of the trap
frequencies as well, although this requirement is relaxed for long time-of-flight measurements. We
experimentally determine ωx , ωy , and ωz by inducing a sloshing motion in the clouds center-ofmass and measuring their position after a variable amount of time. The confining potential stems
from the dipolar interaction between the atoms and the 1064 nm laser forming the optical dipole
trap. The 1064 nm light is red-detuned from the D1 and D2 transitions to the trapping potential
is attractive with respect to beam intensity [57]. The ODT is is deepest at the intersection of
the focused ODT beams, which are Gaussian in profile and approximate a parabola near the trap
minimum.
The vertical slosh is induced by briefly turning off the H2 beam, allowing the atoms to fallaway from their equilibrium position before the potential is restored and the cloud oscillates about
its initial position as shown in Fig. 2.12b. However, for H2 powers near the lower limit required
for sustaining against gravity, this method has its challenges. The power of H2 (and the other two
beams) is determined by a waveform generated by an arbitrary waveform generator; because thesmallest programmable time step is far too long to directly program the chopping of H2, which
generally lasts <500µs, the chop is performed by digitally chopping the AOM RF for the desired
impulse time. As a result, the servo rails and when the AOM is switched back on the power spikes
to its maximum possible power for a time before oscillating as the servo reestablishes the desired
setpoint. This provides a slosh impulse much more powerful than desired at low power settings andportions of the cloud may slosh right on out of the trap or into regions of H2 which are anharmonic.
This leads to difficulty in characterizing the trap at low trap frequencies. Less forceful sloshes may
be attained by simultaneously chopping the servo setpoint and the H2 AOM, since the servo does

36
not rail and only small residual power oscillations which very quickly damp away occur following
the restoration of the trap. This allowed us to reduce the impulse time to < 150 µs and measure
trap frequencies as low as 2π × 45 and with more accuracy and repeatability. The calibration of ωz
is shown in Fig. 2.12a.
The horizontal trap frequencies are measured in situ using the vertical imaging system. To
induce slosh in the x̂, and ŷ directions an additional beam (one of the X-Lattice beams, to be
described below) is intentionally misaligned and ramped up slowly to shift the cloud’s equilibrium
position. The perturbing beam is then extinguished and the cloud returns to and oscillates about
its equilibrium position as shown in Fig. 2.13b. The misalignment of the lattice beam is carefully
done with a micrometer to ensure it may be returned to its proper alignment. We are able tomeasure the horizontal trap frequencies in both principal axes simultaneously. The calibration
curves for ωx and ωy (Fig. 2.13a) are fit with a different functional form than the vertical frequency
ωz (Fig. 2.12a) due to the absence of gravitational effects on the horizontal confinement.

2.7

Optical lattices
In order to understand how implosions will amplify intrinsic density fluctuations we intro-

duced density fluctuations of our own by way of a weak optical lattice. Three lattices, one in eachof the x̂, ŷ, and ẑ directions, each formed by two running wave beams intersecting at a small angle
were constructed to provide the control perturbations. The design and alignment of the lattices
posed a tricky engineering problem due to the limited optical access to the science cell, see Fig.
2.8 All three lattices must be detuned with respect to each other and to the H2 and V beams.
Following the transfer of atoms from the H1 loading beam to the crossed ODT formed by H2
and V, a flipper mount places a half wave plate prior to the PBS in the H1 beam path on thebreadboard containing the H1 and H2 beam paths prior to their entry into the science cell; see
Fig. 2.9. With the lattice diversion waveplate introduced, the PBS reflects the 1064 nm light to a
periscope directing the light to the top breadboard, which is shown in Fig. 2.10. Here the beam’s
polarization is cleaned with waveplates and a PBS, before passing through a beam reducer. The

37

(a)

(b)

Figure 2.12: (a) Calibration of the vertical trap frequency vs the setpoint; the beam power is linear
in the setpoint. Below a certain beam power the ODT cannot hold the cloud up against gravity. The
functional form of the frequency dependence is altered by gravity as the trap frequency depends onthe curvature, i.e. the second derivative of the potential. Measurements from two months in 2021
are shown in (a). The difference between the two calibration curves is due mainly to intentionalchanges in alignment of the H2 beam. (b) Example of a vertical slosh measurement. Each datum
represents one experimental run as the measurements are destructive.

(a)

(b)

Figure 2.13: (a) Calibration of the horizontal trap frequency vs the setpoint; the beam poweris linear in the setpoint. The slosh frequency scales as the square root of the power. The X
frequency fit contains only one free parameter and intersects the origin. The Y frequency contains
a contribution from the H2 beams transverse profile and so the Y frequency fit intersects at the Yaxis at 21 Hz. (b) Example of a horizontal slosh measurement. Each experimental cycle yields an
X and Y position.

38
beam then passes through a 70:30 BS and the 30% transmitted is coupled to a PMSM fiber and islaunched into the Z-lattice optics depicted in Fig. 2.14 The 70% reflected light immediately passes
through an additional 70:30 BS and its 30% transmitted light passes through a broadband AOM
and the first sideband coupled to an additional PMSM fiber which is launched into the Y-lattice
optics also shown in Fig. 2.14. The light which is reflected from the second 70:30 NPBS passes
through an AOM. The first sideband of this light passes through a 50:50 NPBS and a 90:10 NPBS,
such that the two reflected beams have similar powers. The portion of the X-lattice light which istransmitted through the 90:10 NPBS is aligned to the lattice monitor PD used for calibration andmonitoring. The two beams which form the X-lattice are incident upon a three inch mirror with a
45 degree hole bored through it to allow the V beam of the ODT to pass through; the V beam isthen flanked by the parallel X-lattice beams and they are subsequently diverted down toward thescience cell by a SWP dichroic mirror through which the top probe beam passes, see Fig. 2.10. The
X-lattice beams are steered by the same +125 mm lens above the science cell shown in the inset
of Fig. 2.11, but are not shown in that figure. As mentioned above, one of the X-lattice beams isused to induce slosh in the X and Y directions; this is done by misaligning the micrometer on thethree inch mirror which parallelizes the X-lattice beams with respect to the V ODT beam. This
misalignment is done in a controlled and reversible way by the high-tech trick of taking a photoof the micrometer before turning it. The geometry of the lattice beams will be discussed in moredetail in Chapter 5.

39

Figure 2.14: Optics for Y and Z lattices.

Chapter 3
Image Apparatus and Analysis

The primary tool which we have to learn about our atomic system is absorption imaging(AI). Absorption imaging allows us to directly measure the spatial distribution of our atomic
cloud. How the atoms move over time and their various responses to different conditions allows
us to learn other properties, such as temperature and the oscillation frequency of our trap among
many others. The basic premise of AI is simple: we briefly shine laser light upon our atomic
cloud, which absorbs many of the photons; then a bit later we again briefly shine the laser after
the atoms are gone. Simplistically, the difference between the pictures, first when the atoms are
there and thus leave a shadow, second when they are not and have no shadow, tells us where theatoms are. A background picture with the laser off is also taken to subtract out background light.
The probe frequency, atomic spectrum, optical depth, background light, polarization, scattering off
dust, timing, focus, and many other things influence the resulting absorption images. The detailsof AI are well established [81, 103, 131, 27], but given the central role AI in our experiment andthe fact that all images are processed by analysis software written by myself, I believe a thoroughdiscussion on the details is warranted.

3.1

Absorption Imaging: the Basics
As the name suggests, absorption imaging relies on probe light being absorbed by the atoms.

As the light passes through the atomic cloud, it is absorbed and re-emitted by the atoms in somerandom direction. If we have very few atoms, most of the photons from our probe laser will pass41
unaffected through the cloud. If we have many atoms, then many photons are absorbed. The
intensity of the probe thus depends on, among other things, how many atoms through which it
passed. Therefore, as the probe laser beam goes through the cloud, it will become dimmer. Thisattentuation of the probe may be described by the mathematical formula known as Beer’s Law or
the Beer-Lambert Law,
I = I0 e−OD ,(3.1)

where OD is the optical depth (OD) of the atom cloud the laser beam has passed through and I0 is
the initial intensity of the beam. The OD of the atomic cloud is one of the most important pieces ofdata in the lab. While the OD is indeed calculated from a set of images, not a raw image in and of
itself, it is nevertheless treated in lab like a simple picture (when things are going right). From thecalculated ’OD image’ we extract everything we can know about the cloud. Boiling it all down, theonly thing that determines the measured OD is how much light was absorbed/scattered as it passed
through the imaging system. If one knows the properties of your light, atoms, trap, and optics
very well, then the only free parameter is how many atoms the light had to travel through. This isexpressed as a column integrated density n(x, y). Assuming the laser travels in the z direction withinitial intensity I0 before the cloud, then (assuming a collimated probe) the intensity just after the
cloud is:I = I0 e−n(x,y)σ ,(3.2)

where σ is the absorption cross section of our atom [55]. The cross-section contains the information
about how likely the photons are to be absorbed by the cloud. Since n(x, y) has units of inversearea), it is clear that σ must itself have units of area. We will discuss the cross section more later.
Clearly we can see from Eqs. 3.1 and 3.2 thatOD = n(x, y)σ.

(3.3)

Therefore given we understand our atom’s response perfectly as encapsulated in σ, determining the
OD is really determining the density! So now how specifically do we determine the OD? If we had
no background light whatsoever then all of the light would come from the probe beam. Ignoring42
the background light for now, we take two images and consider the probe intensity for the same
pixel in those two images. That light has intensity Ishadow in the shadow frame when the atoms
are present. When the atoms are not present during the light frame, that light has intensity Ilightsince no light was absorbed. In the language of Eq. 3.1 its obvious that I0 = Ilight and therefore
Ishadow = Ilight e−OD . Solving then yields:
OD = ln

IlightIshadow(3.4).

However, in the real world, there will always be some background light. This can be caused by actual
background light from non-laser sources or it can be background currents or thermal fluctuations on
the camera detector. Let us assume that the response of the camera is a contribution linear in theprobe intensity plus a constant background term, then to isolate the laser light in the light frameand shadow frame we simply note that Ilight laser = Ilight − Idark and Ishadow laser = Ishadow − Idark
where Ilight , Ishadow , and Idark are apparent intensities that the camera tells us in some units(frequently just ‘counts’), and Ishadow laser , Ilight laser are the corrected intensities actually due tothe laser (still in the cameras units). But remember, in determining the density we only care about
the OD due to the atoms absorbing the laser, thus:
OD = ln

Ilight laser
Ilight − Idark= ln.
Ishadow laser
Ishadow − Idark(3.5)

On the right hand side Ishadow , Ilight , Idark are the actual camera image files that I work with on
the computer. Presuming that we know σ quite well we can use Eq. 3.3 to get the density andhence the number [103]Z ∞Z ∞
N=
∞

3.1.1

1
n(x, y)dxdy =σ∞

Z ∞Z ∞
ODdxdy.
∞

(3.6)

∞

Ideal Beer’s Law

Consider a laser beam propagating along z incident on a collection of atoms. Let us assume(for now) that each atom has an identical cross-section σ0 . If we consider the atoms as a series of43thin sheets with some atoms in each sheet with thickness dz, then the change in intensity of thelaser field is only on how many atoms lie within that sheet:(3.7)∆I = Iin Natoms = Iin na (x, y, z)σdz.

Where the density of the atoms is na (x, y, z). Rewriting this in terms of an absorption coefficientκ = n(x, y, z)σ0dI
= −κ(na (x, y, z))I(z)
dz

(3.8)

shows how the absorption coefficient is independent of the intensity. This is the meaning behind
our constraint that the cross-section σ is independent of I(z). Omitting the transverse coordinates,
we proceed by rearranging and multiplying by an integrating factor e 0 κ(z )dz
Rz

0

0

Rz
dI(z) R z κ(z 0 )dz 000e0+ κ(z)I(z)e 0 κ(z )dz = 0.
dz

(3.9)

The product rule may be applied in reverse. Integrating and rearranging yields:

Rz
d 
0
0
I(z)e 0 κ(z )dz = 0
dz
Rz

0

(3.10)

0I(z)e 0 κ(z )dz − I(0) = 0
Rz0I(z) = I(0)e− 0 κ(z )dz
Identifying ODBeer =0000 κ(z )dz yields Beer’s Law in its simplest form

Rz

I(z)I0 e−ODBeer .

(3.11)

This only applies when the absorption coefficient κ is independent of I [55]. This relation definesthe ‘ideal’ OD: ODideal = ODBeer = σ0 ncol (x, y). However, life is tough in the aluminum siding
business, and so we must grapple with the realities of absorption imaging in the real world.
3.1.2

Absorption Imaging in the Real World: Bad LightIn the real world, some of the light incident on our CCD is composed of photons which are of
the wrong frequency, polarization, or origin, origin meaning they did not pass through the atoms

44
and imaging system as intended. Inevitably, some of this light shows up on our CCD, and wedenote the intensity of this bad light Ibad . This affects our measured OD and puts a limit on the
maximum OD our imaging can accurately measure [81, 16]:
ODmeas = ln

IL0 + Ibad
IS0 + Ibad

and IL0 is the intensity good photons of the probe beam without atoms present, and IS0 is theintensity of the good photons when atoms are present. The dark frame (ID ) is omitted for simplicity
here. The maximum OD of the imaging system, i.e. the OD when every good photon is absorbed
is thus
ODmax = ln

IL0 + Ibad.
Ibad

We want only the optical depth from the good photons:
ODgood = ln

IL0
.
IS0

(3.12)

We can solve for this in terms of the measured OD and the maximum OD, the latter of which ismeasured by imaging a cloud with a much higher OD than the maximum. This correction yields
ODgood = ln

1 − e−ODmax
e−ODmeas − e−ODmax.

(3.13)

It is best to operate comfortably below the maximum achievable OD to minimize errors in measuring
the OD of clouds near the maximum OD; if ODmeas ' ODmax the denominator in Eq. (3.13)
becomes vary small, and small sources of noise can dramatically influence the inferred ODgood .
3.1.3

Absorption Imaging in the Real World: Saturation

The following derivation follows sections 7.6 and 7.6.2 of [55]. The attenuation of the probe
beam as it passes through the atoms is generallydI= −κ(ω, z)I(z) = −n(z)σ(ω, I(z))I(z).
dz(3.14)

The cross section’s dependence on the frequency ω of the probe light is [55, 36, 66, 127, 125, 126, 124]:
σ(ω) = σ0

Γ2 /4.
(ω − ω0 )2 + Γ2 /4(3.15)

45

Figure 3.1: Saturation of OD due to bad light. This image does not accurately reflect the true
optical depth of the cloud because of bad light which leaks through. The maximum apparent OD
of 2.5 is due to 8% bad light.

The FWHM of the lineshape Γ is the decay rate of atoms from the excited state back to the groundstate, ω0 is the resonant frequency of the transition ~ω = Eexcited − Eground , and the resonantcross-section σ0 in terms of the transition wavelength λ0 is
σ0 =

3λ20
.
2π

(3.16)

Saturation of the transition also gives the cross-section intensity dependence. Reasoning from
population dynamics and Einstein coefficients, we can introduce a terms to the absorption coefficientwhich accounts for the saturation of the transition. The absorption coefficient is thenκ(ω, I) =nσ(ω)
1 + I/Is (ω)(3.17)

46
where
Is (ω) =

~ωΓ
.
2σ(ω)

(3.18)In my experience, the term saturation intensity is, almost without exception, used to refer to thevalue of Is on resonance [55]. Then
Isat ≡ Is (ω0 ) =

πhc
Γ
3λ

(3.19)

is the saturation intensity. For the cycling transition |2, −2i → |3, −3i in 39 K, this is 1.75 mW/cm2 .
Expanding this out into it’s full glory yields
κ(ω, I) = nσ0

Γ2 /4(ω − ω0 )2 + 14 Γ2 (1 − I/Isat )

nσ0
= 2 2
4δ /Γ + 1 + I/Isat(3.20)

where δ = ω − ω0 . On resonance, Eq. 3.14 then becomesdI
nσ0 I(z)=−dz
1 + I(z)/Isat )ds0
nσ0 s0 (z)=−dz
1 + s0 (z)(3.21)

where s0 (z) = I/Isat [16]. Integrating this can be done as follows:Z soutZ z
1 + s0ds0 = −nσ0
dz
s0
sin
0
ln

sout
+ (sout − sin ) = nσ0 z = −ODBeer
sin(3.22)ODBeer = ODmeas + (sin − sout )ODBeer = ODmeas + sin (1 − e−ODmeas )Where ODBeer is the hypothetical OD when the absorption coefficient has no intensity dependence.
We may express this as
ODBeer = ODideal = ln

L−D L−S
+
,
S−D
Isat

(3.23)

where S, L, and D are the rectangular arrays of the shadow, light, and dark frames respectively
in units of cameras counts. The saturation intensity Isat is also expressed in camera counts, andis dependent on the exposure duration of the probe. When effects like atomic acceleration due tophoton recoil may be ignored Isat will be linear in the pulse duration. For the side imaging system
Isat = 125 counts/µs, while for the top imaging system Isat = 9 counts/µs.

47
3.1.4

Absorption imaging In the real world, final form

Putting everything in this section together to understand exactly how to compute the atomic
column density, let’s write the equations with all terms included.
Iout (z) = Iin e−ODmeas
ODideal = ODBeer = −ncol σ0
IL − Idark
IL + Ibad − Idark= ln 0IS − Idark
IS0 + Ibad − Idark
IL − Idark
ODgood = ln 0
IS0 − Idark
IL + Ibad − Idark
ODmax = ln 0
Ibad − Idark
1 − e−ODmax
ODgood = ln −ODmease− e−ODmax
IL
IL − ISODideal = ODgood +(1 − e−ODgood ) = ODgood +
Isat
Isat

ODmeas = ln

(3.24)

Where IL = ILmeas = IL0 + Ibad is the measured light frame instensity and IL0 is the intensity of
only the photons of the correct polarization and intensity (and similarly for IS , IS0 ).

3.2

Side Imaging
The side imaging configuration allows us to characterize cloud during time of flight expansion.

We apply saturation corrections. All images taken with the side camera are absorption images. In
this context, “side” and “side imaging” refers to the imaging performed with the probe beam which√
propagates horizontally through the science cell with k̂side = −(x̂ + ŷ)/ 2 . A diagram of the side
imaging optics is given in Fig. 3.2. A more detailed description may be found in [27], from where
Fig. 3.2 was taken. The side imaging is used to measure the atom number and temperature in timeof-flight expansion and is always used on resonance with the |2, −2i → |3, −3i cycling transition.
Use of a cycling transition is important for absorption imaging; if an atom ever decays into a stateother than the |2, −2i ground state, the probe light will no longer couple (strongly) to an excitedstate and the optical cross section of that atom will drop to nearly zero. When a cycling transition
is used, the excited state has only one decay path available and the atom will always decay back

48

Figure 3.2: Side imaging system, figure courtesy of Roman Chapurin [27] (Fig. 3.10 in [27] reprinted
with permission).

to the initial ground state.
3.2.1

Calibrating Side probe intensity

Our side imaging intensity can be calibrated in a few ways. One is to rely on the camera
specifications and optical layout to extrapolate the probe intensity at the atoms. Another way is
to measure many clouds with stable and repeatable number and vary the intensity; the resulting
behavior of ODmeas can be fit to a model yielding the saturation intensity in terms of camera
signal, as shown in Fig. 3.3. This method is the most convenient, but fails at low probe intensities
as the signal is lost in the noise. A third method, and the most fool-proof (a necessity when I’m
in the lab) is to directly compare camera response to probe with known intensity. This is done by
reducing the beam size such that it is small enough to comfortably fit within the power meter’s49

Figure 3.3: Optical density vs probe intensity. We fit to a model with the saturation intensity incounts is a parameter. Here the saturation intensity is found to be 602(22) counts for the 20 µs
pulse The cyan points are masked because the OD is washed out when signal to noise is low

photo detector and the CCD region of interest. This can be done with an iris placed sufficiently
close to the camera. The power of the beam is then measured with a power meter for various
probe intensities and corresponding images are taken. It’s critical that exactly the same amountof light be measured by our camera and the power meter; if any light misses the camera, then one
cannot trust the comparison to the power meter unless one can be confident in exactly how muchlight is hitting the camera. The number of counts per pixel per microsecond per mW of poweris then computed from these data. The total counts as a function of measured power provides asanity check on the cameras linearity, but this information can also be obtained by taking pixelby-pixel ratios of images with varying intensities. Intensity dependent behavior pixel ratios are a
red flag for camera saturation. Fortunately no saturation was detected in either the top or side
imaging cameras. For the side imaging, the saturation intensity (assuming no bad light) is 1250

50

Figure 3.4: image of beam used for calibratiing probe intensity. right, camera signal vs probe power
measured via powermeter. linearity is goodcounts for our standard 20us pulses, and we typically operate at 15-40% pf this value and apply
the appropriate correction.

(a)

(b)

Figure 3.5: Example absorption images taken with the side imaging system (a) and the top imaging
system (b). The side image (a) is taken after 17 ms of time-of-flight expansion. The top image (b)is taken in-situ. Images are false color.

51

3.3

Top imaging
The top imaging system has much higher magnification and resolution. The term “topimaging” is so named because the probe beam propagates from the top of the chamber downward
through the atoms (k̂top = −ẑ). The setup and characterization of the vertical imaging optics issignificantly more involved than with the side imaging. Much of the design and characterization
was carried out by Xie Xin and is covered in her thesis [131]. The system is outlined in Fig. 3.6,
with the objective shown in Fig. 3.7. The objective had to be modified from Xin Xie’s original
design due to a pathological aberration stemming from the custom machine work required to satisfy
space constraints. We shall first describe the original objective in the main text before describing
the discovery and investigation of aberrations which forced us to replace the primary lens of theobjective.
3.3.1

Design of vertical imaging system

Our design of the higher-resolution imaging system used for this work is guided by the scaleof the implosions, which is initially determined by the geometry of the optical lattice used to seed
the density perturbations. However, we are also constrained significantly by the geometry of thescience cell, its surrounding coils and antennae, and especially the cart which moves beneath theobjective for each shot. In order to maintain a high level of control along as many degrees of
freedom as possible, the objective is mounted on a mount consisting of an XYZ translation stage
from Newport, an Edmund (66-534) goniometer for pitch control (centered about the objective),and an Edmund rotary stage (66-516) providing roll control.
The objective itself is held in close proximity to the bottom window of the science cell, and isconcentric with the primary coils surrounding the science cell. The objective consists of five parts:
two lenses, two spacer rings, and the holder those four elements are mounted within. The holder
and the two rings are constructed out of eddy-current-free Ultem to avoid stray magnetic fields
which could likely complicate the magnetic field quenches. The original primary lens is a Thorlabs

52

Figure 3.6: Top imaging system used for in-situ imaging.

CNC-polished asphere (AL4532-B) which sets the numerical aperature of the system. The aspheresits below a much thinner plano-convex lens (LA1978-B), and the two are separated by a spacer
ring. The eyepiece used is a 1500 mm plano-convex lens from CVI optics (PLCX-50.8-77.6-C). The
specifications of the imaging system using the original primary lens were determined by previous
work to be the following: Effective focal length of 31.03 mm; NA (defined by the acceptance cone) of
0.384; exit pupil diameter (defined as the image of the stop) of 26.15 mm; limiting aperture (limitedby the spacer ring between the two pieces) of 22 mm; Airy disk radius at 767 nm of 1.13 µm [131].
However, the primary lens was replaced in September 2021 due to unexpected aberrations.

53

Figure 3.7: The objective for the high magnification imaging system, taken with permission fromXin Xie’s thesis [131]. The objective is comprised of a +750 mm plano convex lens and a +30 mm
asphere. Additional information on the original design and testing process may be found in Xin
Xie’s thesis [131]. The new objective has the same holder, retaining rings, and secondary lens but
a different asphere with a similar focal length which did not require turning down to fit the holder.

3.3.1.1

Top imaging aberrations

Unfortunately, we encountered significant problems when using the above design to image
our atoms. The problem first became evident when the lattice was turned on and the fringes of
the atomic density corrugations were located at a dramatically different location on the CCD thanthe bulk cloud itself. A large amount of time, several months, were spent trying to understand thenature and cause of this phenomenon. It was possible to bring these fringes onto the atomic cloud,
but the cloud as a whole was inevitably defocused. It seemed as though the objective was exhibiting
dramatic spherical aberration, well in excess of what anyone in the group expected or though likely.
Only after exhausting every idea we could test with the objective in the experiment did we decideto remove the objective for additional bench testing. Our bench tests confirmed the presence of
substantial spherical aberration, and we switched the high NA asphere from the one machined here
at JILA to the one Thorlabs machined down to size for us. This seemed to help some, but the
improvements were minor. We tested and falsified a hypothesis that the optical coating was to be
blamed when another JILA machined lens exhibited similar aberration. Eventually, a new stock

54

Figure 3.8: Left, out of focus, but fringes on top of cloud. middle, cloud size vs objective Z. rightfocused bulk cloud, but fringes far away.

asphere compatible with the objective holder geometry was found from Edmund Optics (37-434)
which did not require custom machining, and immediately provided the least aberrated images todate. Our conclusion is that the machining of these precision aspheric lenses was the cause of thisaberration. This was surprising as three different lenses exhibited this problem despite the careful
efforts of machinists with experience in machining precision glass. Enough time had been wasted
on this puzzling problem that we were happy to proceed with the new Edmund asphere.
However, our imaging system had been carefully designed with the at the time unquestioned
assumption that the lenses were optically sound. Several months of modeling the optical systemusing Zemax and the objective holder and mount with solidworks had promised a diffraction limited
resolution of 1.1 micrometers. At this point we were behind schedule, and low on manpower, so theoptical layout with the new lens could not be studied with the same level of detail. Fortunately,
the optical lattices provided a shortcut for verifying the magnification and measuring the actual
resolution, which turned out to be 2.5 micrometers, more than a factor of two worse than originallyplanned, and placing the wavenumber of the Y-lattice used for many of the experiments describedin Chapter 5 at the very edge of the resolution.

55
3.3.2

Top camera calibration

The high magnification, in-situ images used to measure the static structure factor in thiswork are all taken using the top imaging system. Immediately below the science chamber sits thehigh NA objective, perched carefully on a mount with five degrees of freedom. Because of thelarge magnification M = 50 of the top imaging system, the probe intensity at the CCD is a factor
M 2 = 2500 lower than at the atoms. This prevents imaging at probe intensities at or below thesaturation intensity because such a weak signal on our camera is lost in the noise. Additionally, theextremely high OD of the in situ clouds, which can be on the order of ODB ≈ 300 would effectively
block all resonant light. For this reason, we operate at I/Isat ≈ 150. The probe intensity wascalibrated in the same manner as described in section 3.2.1.

3.4

Fringe removal with PCA
The OD of the cloud is, as already mentioned, calculated using three frames, one with theshadow of the atoms, one without, and a dark frame to account for common mode noise. If theprobe intensity profile changes between the shadow and light frames, this difference will manifest
as structured noise. Our CCD is limited in speed such that our shadow and light frames aretaken 400 ms apart. Air currents, dust, thermal effects, reflections, etaloning, and more may causesuch deviations. The top imaging system is particularly prone to substantial interference fringes,
which can make the atomic density difficult to discern. As the primary aim of this work is tocharacterize density fluctuations in BEC, fluctuations or fringes from imaging are also of great
concern should they appear at the same spatial frequencies as the fluctuations to be measured.
Fortunately, these fringes rarely appear at frequencies of concern. However, even a low frequency
fringe can have substantial impact on the analysis of that image, generally making it difficult todetermine the size shape and position of the cloud and overall making life harder. Fortunately,
established techniques exist to mitigate this problem, chief among them being principal component
analysis (PCA) [82, 133, 4, 90, 112].

56
Principal component analysis allows one to reconstruct a new light frame for a given shadow
frame given a set of reference frames. A few tens of reference light frames are decomposed into acharacteristic basis set capturing the most important variations from the average. A portion of the
shadow image, away from the atoms, is then decomposed into this basis, and a new light framewith the same noise profile away from the atoms is constructed. Although only pixels without
influence from the atoms are used to determine the relative weights of the reconstructed shadow
image, the characteristic basis set is computed from a reference set with no atoms at all, and so
the reconstructed image fills in the region with atoms as though there were non. The end result ismuch closer to what the shadow frame would look like without the shadow. This is exactly what
an ideal light frame is, and yields dramatically reduced fringes. A total of M images are taken, so
M light, and M shadow (and M dark) which form both our reference set and the data set, as you
can use all your light frames as the reference pool for each shadow frame. It is essential that theimages be generally similar, e.g. if you change the size of the probe beam or alter the magnification,an accurate reconstruction will be difficult, PCA is not magic, we are trying to recreate what the
shadow frame would look like without the atoms. Each image is j × k pixels, typically about 5002
pixels total. Before reconstructing ideal light frames to match our shadow frames, we compute thebasis we will use for said reconstruction [112, 133, 82]. Let the reference images (the light frames)
be Li for i = 1...M and flatten these into column vectors Li (jk × 1)let
M

hLi =

1 X
Li
M
i=1

be the average and let dLi = Li − hLi be the deviation of the ith reference vector Li . So then let
A = [L1 , L2 , ..., LM ]be a matrix. A is thus jk rows by M columns as described above. We now compute the matrix
C = AAT ,Which is a huge square matrix, letting p = jk be the number of pixels in each frame, C is p×p totalelements. Each image on the top camera has p = jk = 4962 pixels so C therefore has 4964 > 60

57
billion total elements. After converting to double precision, each element is eight bytes, so C is anover 400 gigabyte matrix. The eigenvectors of C are the so-called eigenfaces of the images [82], thatis, the eigenvectors of C, when reshaped into a j × k matrix are the typical noise patterns in the
M reference light frames, and those eigenfaces with the largest eigenvalues are the most prominent
noise patterns. So the eigenvectors of C provide the perfect basis from which to reconstruct our
new light frame. However, C is so large MATLAB won’t let me compute it. This makes me sad.
However, it turns out that the eigenvectors of C = AAT and C 0 = AT A are related:
C 0 Vi = AT AVi = λVi
AAT AVi = AλV i(3.25)

C(AVi ) = (AAT )AVi = λ(AVi ) = λUi
Therefore, we can compute these eignevectors via Ui = AVi where Vi are the eigenvectors of
C 0 = AT A which is much smaller; C 0 has size M × M . After normalizing each eigenface, we have an
orthonormal basis with which to characterize the noise patterns in a shadow frame. Normalization
is not a necessary step, but makes compensating for the masked out region of the shadow moreconvenient.
Consider one shadow frame S1 and the ensemble of of M reference images Li (some of which
are a closer match to S1 than it’s original light frame L1 , and so we could get some improvement
by using a different light frame which produces fewer fringes in the OD calculation, but we can dobetter). We then compute the deviation of S1 from the average Si = S1 − hSi. To ensure we do
not let the atoms influence the reconstruction we mask the atom region in S1 . The inner products(3.26)

S1 · Ui = cigive the decomposition coefficients ci of the shadow frame deviations from the average into thecharacteristic noise patterns Ui . Then the reconstructed light frame can be computed
L1 = hLi + hDi +

M
X
i=1

ci Ui .

(3.27)

(b)

(e)

(c)

(d)

Figure 3.9: Principle component analysis applied to an in-situ image. The original OD image computed from the using the original lightframe is shown in real space in (a) and in Fourier space in (b). After applying PCA and reconstructing a better matched light frame, thefringes are reduced, through not completely eliminated. The first ten eigenfaces/PCs used to reconstruct the light frame are shown in (e).

(a)

58

59
We employ an additional minimization routine wherein the OD is computed for the new lightframe as we vary the coefficients ci and the noise away from the atoms is minimized as a functionof each ci . This fine-tunes the coefficients and provides an even more well matched reconstruction
L1 .
3.4.1

Nature of fringes

In this section we discuss how fringes may appear in absorption imaging. By fringe we mean
a ripple pattern on either one of the raw images used to compute the optical depth or on the
OD image itself. A fringe which appears with perfect repeatability between two raw images will
subtract out and the resulting OD image will not exhibit that fringe. The relationship between
fringe amplitude and optical depth is not in general a linear one, as fringes can be caused in a
number of different ways, and may even be dependent on the atomic density itself.
Suppose there is a piece of dust on the edge of the probe after it has passed through theatoms which scatters light into the camera, including (but not only) the same region of the camera
as the shadow of the atoms appears. If this piece of dust is stuck on a mirror and moves slightly
from shot-to-shot then it will produce fringes of a consistent and characteristic pattern and willappear on light frames taken with no atoms present in a similar fashion to shadow frames where
the atoms cast their shadow. In this case, the constructed light frame will accurately reconstruct
the fringe and it will subtract out nicely.
If the piece of dust is directly in the shadow of the atoms (on the mirror with the dust,
not the camera), and if the imaging were otherwise perfectly repeatable, then the fringes will only
appear on light frames, but always subtract between two light frames. In fact, all OD images takenwith no atoms would produce perfectly flat, featureless OD profiles since we’ve asserted that theimaging in the absence of atoms are identical. In the limit of infinite OD, the shadow frames would
not possess the fringe from this piece of dust as it is in the shadow of the atoms. In this case, abasis set of reference frames with no atoms would not be usable to reconstruct a light frame tomatch the shadow frame closely because the shadow frame lacks the fringe but all light frames used60
as the basis have the fringe.

3.5

Fitting the Images
The density distribution as measured by imaging the cloud is the observable at the core of all

our analysis. Generally, the cloud is only partially condensed, and accurate characterization of thecondensate vs thermal components of the cloud is critical in determining the density, temperatureand other properties vital for thorough understanding. This section describes the functional forms
we use to fit our imaged density distributions and how we subsequently extract information about
the atom number of both the condensate and the thermal cloud and the temperature of the system
as a whole. The results of this section are specific to the case of interacting Bose gases in athree-dimensional harmonic potential. We will begin by presuming that the thermal cloud is welldescribed by the Bose distribution and that the condensate is in the so-called Thomas-Fermi regime[102]. In this case the three-dimensional density profile for the entire cloud ntot (r) is [113, 10]:
g3/2 (z(r))
x2y2z2
ntot (r) =+ nc,0 1 − 2 − 2 − 2(3.28)
Rx Ry
Rz
λ3T
p
where nc,0 is the peak value of the condensate density, λT = 2π~2 /mkB T is the thermal wavelength, and z(r) = exp (µ − V (r))/kB T is the fugacity in terms of the chemical potential µ. The
Bose function gν is defined by the equationgν (z) =

∞
X
zn
n=1

nν

.

(3.29)

The form of Eq. 3.28 neglects the mean-field effects of the condensate on the thermal atoms (andvice versa) which will perturb the potential experiences by the thermal cloud from the external
harmonic potential. The parabolic form of the condensate density is invalid near the edge of thecloud where the density experiences a sudden bend, but describes the bulk density quite well. For
the thermal cloud, we assume z(r) = 1 since at the transition temperature µ = 0. Therefore thetotal density is given by






nth,0,3D
1 x2
y2
z2
x2
y2
z2
ntot (r) =
g
exp −
+
+
+ nc,0 1 − 2 − 2 − 2
g3/2 (1) 3/2
2 σx2 σy2 σz2
Rx Ry
Rz

(3.30)

61
where nth,0,3D is the thermal (subscript th) peak (subscript 0 indicating r = 0) density for thethree-dimensional cloud (subscript 3D). The widths σi are related to the temperature by
1
σi =
ωi

r

2kB T
.
m

(3.31)In Appendix C we discuss the density distribution in more detail and consider the ballistic expansionof the gas. In the absence of collisional or interaction effects upon the thermal component, thethermal component will retain the form of Eq. 3.30 with the size changing in time asσi (t) = σi0 (1 + ωi2 t2 )(3.32)

which through Eq. 3.31 allows the temperature to be measured via
Ti =mωi2 σi2 (t).
kB (1 + ωi2 t2 )(3.33)

When the cloud is in thermal equilibrium the temperature is obviously the same in each direction,
so Eq. 3.33 is a useful check to ensure that the cloud is indeed in thermal equilibrium. It is possible
for a sufficiently rapid change in conditions, such as dramatically reducing the trap depth along
the z axis while maintaining strong confinement along x and y, to require a significant amountof time to reestablish an isotropic momentum distribution and thus reestablish a well defined
temperature. Significant in this context meaning significant in an experimental context and maychange depending on the goals of the experiment. Generically, a few (∼3) collisions must occur on
average for cross-dimensional equilibrium to occur, and a few trap periods to ensure this occurs
globally [5, 89]. We should note that the σi are not RMS widths of the thermal cloud because theprofile is not a Gaussian, but rather a sum of Gaussians according to Eq. 3.29. If only the firstterm in the sum is kept then the profile is purely Gaussian.
3.5.1

Column densities

Of course the imaging procedure yields a two dimensional image as the probe beam actually
measures the column density along the imaging axis. The fitting functions we actually use to fit the62data are thus obtained by integrating Eq. 3.30 along one dimension. The integral of the thermalpart isZ∞
nth,0,3D X 1
mx2 my 2 mz 2
exp − 2 −−dz
g3/2 (1)
2σx
2σy2
2σz2
m3/2
m=1


∞ Z
nth,0,3D X
1mx2 my 2 mz 2=
exp − 2 −−dz
3/2
g3/2 (1)
2σx
2σy2
2σz2
m
m=1


∞ Z
nth,0,3D X
1mx2 my 2
z 02 dz 0=
exp − 2 −− 2 √g3/2 (1)
m2
2σx
2σy2
2σz
2
m=1


∞
X
nth,0,3D √1mx2 my 2
=√
2π σz
exp
−
−
m2
2σx2
2σy2
2 g3/2 (1)
m=1



nth,0,2D
mx2 my 2
ntherm,0,2D (r) =g2 exp − 2 −g2 (1)
2σx2σy2

nth,0,3D
g (z) =
g3/2 (1) 3/2

Z(3.34)

where in the last line we simply defined the peak column density to be ñth,0,3D with normalization.
The thermal column density is therefore fit with three terms nth,0,2d , σx , σy . The condensate term
is quite straightforward yieldingnc,2D (x, y) = nc,0,2Dx2
y2
1− 2 − 2
Rx Ry

3/2

(3.35)

with the three fit parameters nc,0,2D , Rx , Ry for the peak density and size of the condensate. The
total column density is the sum of these two components
ncol (x, y) = nth,2D (x, y) + nc,2D (x, y)(3.36)

In time of flight it may be shown (see Appendix C)that both the thermal component and thecondensate maintain their functional form, with only the length scales changing [73, 26, 40, 113].

3.6

Potential problems
Due to the time constraints following the receipt of the new primary lens used to gather thedata we shall present in Chapter 6, the high-resolution top imaging system used in the implosion
experiments may have been insufficiently or inaccurately characterized. An enormous amount oftime was spent on pathological imaging defects in a series of custom cut lenses, and when the lenswith which the data discussed in this chapter arrived, we proceeded with all deliberate haste to63data collection. Additionally, the condensate itself is so optically dense that we used extremely
high probe intensities to compensate, and the condensate always appeared aberrated in the insitu images as indicated by the presence of rings and ripples around the condensate. The odd
appearance of the Y-lattice peak in Fourier space (see Fig. 5.13a) may be a symptom of a poorly
understood point-spread function.
Fundamentally, we are limited by our incomplete knowledge of the point spread function. The
measurement contrast suppression factor we shall discuss in section 5.5 is essentially our attempt
to measure the point spread function of our imaging system for two wave numbers only. If theresults of low vs high density X-lattice calibrations are accurate then our imaging system displays
a density dependent point spread function. In Figure 3.10 the Fourier profile of the condensate
is compared with that of the thermal cloud by controlled which portion of the density profile isunmasked in computing the FFT. The resolved disk for Fig. 3.10a appears noticeably smaller thanin Fig. 3.10b, suggesting that the resolution of the imaging system is different for different partsof the cloud. If the density distribution itself is a major contributor to the point spread functionof the imaging system, then a full understanding of the imaging system requires a sophisticatedmodel outside the scope of this thesis.
Any imaging system will additionally possess noise. The fringes discussed in Section 3.4 are
one type of noise we encountered, and the saturation intensity correction discussed in Section 3.1.3
may also introduce noise. Since we are in the business of measuring density fluctuations, which one
might consider “noise” in other contexts, we are faced with the problem of distinguishing genuine
density variations on the cloud from apparent variations which in reality are imaging noise and do
not correspond to real density fluctuations. The influence of imaging noise on our results will be
discussed in more detail in Section 6.1.

64

(a)

(b)

Figure 3.10: Average Fourier space profiles of the high-temperature low-density tImp = 0. In (a),only the region of position space containing the condensate was unmasked, the typical treatment forthe data in this thesis. In (b), an annulus in position space was selected to correspond to thermal
atoms, with the central condensate masked. Both (a) and (b) are derived from the same set ofimages, yet the edge of the plateau, which represents the resolution of our imaging, is at a larger
wavenumber for the thermal atoms. The features at the top left and bottom right corner of each
image, and the arrow-like feature on either side of (b) are imaging artifacts which may be ignored.

Chapter 4
Fluctuation Theory of BEC

The goals of this chapter are twofold. The overarching conceptual goal is to understand thenature of fluctuations, and in particular why fluctuations in three dimensions are experimentallychallenging even in the absence of technical limitations. The other goal is to provide the quantitative
predictions for the fluctuations of a 3D bulk BEC in terms of the experimental data; the structure
factor provides the connection between experiment and theory.

4.1

Criterion for BEC
The standard definition of a BEC is a system with a macroscopic occupation of the single

particle ground state. This definition is clear and unambiguous in non-interacting systems, butwhen interactions are present the single particle states are modified by interactions with all theparticles present and this definition, while still fundamentally valid, becomes more complicatedto define. Here we formalize the concept of Bose-Einstein condensation in terms of the one-bodydensity matrix. We shall see that the BEC corresponds to the existence of off-diagonal long rangeorder in the one-body density matrix and, equivalently, the macroscopic occupation of one of thesingle particle states. The treatment here resembles that given in standard texts on BEC andsuperfluidity [102, 98]. The one-body density matrix is defined:n(1) (r, r0 ) = hψ̂ † (r)ψ̂(r0 )i(4.1)

The field operators ψ̂ † (r) and ψ̂(r0 ) are the creation (at r) and annihilation (at r0 ) operators inposition space. The meaning of the average h...i is not identical for all situations. If the state is

66
pure state then the average is the standard inner product of the state over its Hilbert space. If
instead ψ̂ represents a statistical mixture of possible states, then the average is an ensemble average
according to the relative probabilities of each state. It is in the second sense in which the density
matrix will generally be used in this thesis.
The diagonal of the one-body density matrix is simply the density of the system n(1) (r, r) =Rn(r) and obeys the normalization n(r)dr = N . The off-diagonal components represent the correlations between different parts of the the system. We may define the momentum-space representationof the one-body density in exactly the same way as in position space(4.2)

n(1) (p, p0 ) = hψ̂ † (p0 )ψ̂(p)i

where the momentum-space field operator is the Fourier transform of the coordinate-space field
operatorZψ̂(p) =

ψ(r)exp (−ip · r)dr.
(2π~)3/2

Evaluating Eq. 4.2 for p = p0 yields the momentum distribution n(p)1n(p) =(2π~)3Zss  −p·s/~n(1) R + , R −
e
dRds
2
2

and where s = r − r0 . The normalization is the same as in coordinate space
4.1.1

(4.3)
R

n(p)dp = N .

Uniform systems

We first consider a uniform ideal system of N particles in a volume V . In such a system theeigenstates are simply plane waves:
eip·r/~
φ(p) = √
.
V

(4.4)

The volume V and the number N will generally drop out in the thermodynamic limit wherein those
two extensive quantities N , V → ∞ while their ratio N/V = n is held constant. We can express
the density matrix as a sum over momenta:n(1) (r, r0 ) =1 X
0Np eip·(r−r )/~
V p

(4.5)

67
where Np are the occupation numbers for each momentum state. At larger separations |r0 − r| only
the p = 0 mode contributes while all other terms interfere to zero and we havelim

|r0 −r|→∞n(1) (r, r0 ) =

N0
.
V

(4.6)

Where N0 is the occupation number of the zero-momentum state. This motivates using non-zero
asymptotic behavior of the one-body density matrix as the criterion for Bose-Einstein condensation.
4.1.2

Trapped systems

Of course an alternative formulation is needed for spatially finite systems encountered in thelab since the infinite limit is obviously problematic. The customary approach is to expand thedensity matrix in terms of its eigenfunctions φj (r) satisfying the equation:Zn(1) (r, r0 )φj (r0 )dr = λj φj (r).

(4.7)

The density matrix is positive definite and Hermitian, ensuring its eigenvalues λj are positive asRthey must be since, under the normalization convention φ∗i φj dr = δij , we identify λj = nj as thePoccupation of the single particle wavefunction φj (r); additionally i ni = N . We may express thedensity matrix in the diagonalized formn(1) (r, r0 ) =Xnj φ∗j (r0 )φj (r).

(4.8)

j

At zero temperature, the ideal gas possesses exactly one non-zero eigenvalue, that of the lowest
energy single-particle state which has eigenvalue N indicating all the particles are condensed intothat state. Therefore, for an interacting system, we require that the occupation of one of the states
be of order N and all others finite and of order unity in the thermodynamic limit N → ∞ forBose-Einstein condensation to occur. In our lab, N0 approaches N for T → 0, although even at
T = 0 interactions will deplete the zero-momentum state such that N0 /N never quite reaches unity.
For superfluid Helium, N0 /N ≈ 0.1, a sufficiently large portion that it’s reasonable to regard thatsystem as being at least partially condensed.

68
One may reasonably ask “if we’ve just argued that long-range off-diagonal behavior of thedensity matrix is our criterion, but this doesn’t apply to trapped system, what was the point ofconsidering the long-range off diagonal behavior of the one-body density matrix?”. A fair question,
and we shall return to it later in more detail when the connection between the phase correlation
function and the density matrix is established. However, while the infinite limit may be irrelevant
to trapped gases, the meaning of long-range is still perfectly coherent in a spatially finite system.
It will turn out that the opposite regime from BEC the density matrix drops off exponentially overa microscopic distance and with an asymptotic limit of zero. It’s true that in a finite system thedensity matrix necessarily goes to zero as |r − r0 | → ∞ since there are no more atoms beyond acertain distance, but the density matrix may drop off non-exponentially or over an intermediate
distance, so the criterion still has value.

4.2

Presence and absence of BEC in 1, 2 and 3 dimensional ideal gases
This section aims to to provide a cohesive argument for the presence or absence of Bose-

Einstein condensation in one, two, and three dimensions in terms of the behavior of the chemical
potential. This approach compliments arguments regarding the existence of long-range order alluded to in the previous sections. We will return to the question of long-range off-diagonal orderand the connection to phase fluctuations in later sections, but here we understand BEC as a state
where one of the low-lying states is microscopically occupied. Here we consider ideal gases. We
begin in section 4.2.1 with a presentation of results about Bosonic systems which are true as aconsequence of Bosonic physics (e.g. the Bose-Einstein distribution function). In section 4.2.2 weshall introduce results, themes, and definitions which will be common to all systems discussed inthis section, and we do so for the purpose of convenience and foreshadowing. We will then discuss
uniform ideal systems in Sec. 4.2.3 and then consider the influence of harmonic external potential
in Sec. 4.2.4.

69
4.2.1

General facts of Bose systems

The total number of particles in a Bose gas is given in general by [12]
N=

X
i

1
e(i −µ)/kB T − 1(4.9)

where the sum is over each individual eigenstate i . Converting the sum to an integral we haveZ ∞N = Vd0dd ke(−µ)/kB T − 1 (2π)d1(4.10)

where Vd is the system volume in d dimensions. These integrals may be directly evaluated over thewavenumber as in Eq. 4.10, but an alternative approach in terms of the density of states D() can
also be illuminating. In terms and integral over the density of states the total particle number is:Z ∞N = Vd0D()de(−µ)/kB T(4.11)The density of states in one, two, and three dimensions are given in Appendix A for the uniform
gas and for harmonic confinement.
4.2.2

Foreshadowing

In the following sections, the chemical potential µ will be evaluated for D = 3, D = 2 and
D = 1 dimensions according to Eq. 4.10 and Eq. 4.11. The behavior of µ in the highly degenerate
regime will be determined and the occupancy of the ground state N0 thus estimated according to
N0 ≈ kB T /|µ|. We will then evaluate whether or not N0 ≈ N is possible in the thermodynamic
limit. In anticipation of the forthcoming results it is useful to have certain facts in mind from theoutset. The chemical potential µ is a central parameter for these systems, and for a classical orBose gas µ ≤ 0, as long as we define E = 0 to be the energy of the system’s single-particle groundstate. At zero temperature for all non-interacting cases described here the chemical potential will
vanish µ(T = 0) = 0 − O(1/N ). The chemical potential µ is often expressed through the fugacityz = exp(µ/kB T ) where z ∈ (0, 1) since µ ∈ (−∞, 0). Throughout my thesis, a statement that µ = 0
is, unless otherwise stated, ignoring the 1/N fractional corrections to µ. As a consequence of the70form of the Bose-Einstein distribution a particular function, referred to in this thesis as the Bosefunction1 gν (z) will be a recurring feature [12]. The Bose function gν (z) is defined bygν (z) =∞X(4.12)

z n /nν .

n=1

Finally, the degeneracy parameter (also called the phase-space density) of the gas is defined forq
2π~2
a d-dimensional gas with particle density nd as the product nd λdT where λT =mkB T is thethermal wavelength. When nd λdT  1 the particle wavefunctions are all overlapping and the gas
is considered deeply degenerate. In the opposite limit, nd λdT  1, the particles lose their wavelike
qualities and behave as classical free particles well separated from each other and thus quantum
effects are less important.
4.2.3

Uniform ideal system

With these facts in mind, we turn to the specific case of uniform, ideal, Bosonic systems.
We shall begin with three dimensional systems before turning to two and finally one dimensional
systems. The quantitative results for each case will be consolidated and listed in summary at theend for each case.
4.2.3.1

3D uniform potential

In three dimensions equation 4.11 becomes problematic as in 3D the density of states is√
D3U () = 2m2  /(2π 2 ~3 ) which assigns zero weight to the ground state  = 0. Hence a correction
is required and the total number is given byz
V m3/2N=+√
1−z
2 π 3 ~3

Z ∞
0

√



e(−µ)/kB Td(4.13)

The correction term z/(1 − z) is identically equal to N0 according to Eq.4.9 since the ground statehas  = 0. A detailed derivation of the correction term z/(1 − z) can be found in [12] and is notincluded here. A naive implementation of Eq. 4.10 converges perfectly fine but does not accurately
1

The Bose function g(z) is also known as the polylogarithm, Liν (z) and sometimes the Bose-Einstein integral or
the Fermi-dirac integral when expressed in terms of integrals over those distributions

71
account for the ground state, although this failure is less clear in the momentum integral of Eq. 4.10
compared to the energy integral of Eq. 4.11. Therefore, we use the replacement N → N − N0 andproceed, finding
g3/2 (e−µ/kB T )g3/2 (z)N − N0
=
=
3
V3
λT
λ3T

(4.14)

where N0 is the number of particles in the ground state, is the thermal wavelength and is the Bosefunction. The case2 of g3/2 (1) = ζ(3/2) ≈ 2.61 shows how in three dimensions the population of
the excited state is bounded from above to a finite value. The chemical potential is bounded fromabove by zero, and so the fugacity z ≤ 1. Since g3/2 (z) converges at z = 1, the number of excited
particles is limited by Eq. 4.14 and any additional particles must therefore occupy the groundstate, which has  = 0 and hence the chemical potential µ is zero when N exceeds the value given
by Eq. 4.14.
What temperature will this saturation of excited states take place for a given atom number
N ? We determine T3D by setting N0 = 0 and µ = 0 which are approximately true at the onset of
Bose-Einstein condensation, and solving Eq. 4.14 for T . The result is
2π~2
T3D =
mkB



2/3
N
ζ(3/2).
V(4.15)

We can describe the degeneracy of the cloud according the the phase space density n3D λ3T
which is well below unity for non-degenerate systems but exceeds it in degenerate systems. For
temperatures well above the transition temperature, n3D λ3T  1 and g3/2 (z) ≈ z so that thechemical potential may be approximated:µ(T  T3D ) = kB T ln(n3D λ3T ).

(4.16)The behavior of µ changes dramatically at T = T3D . As T → T3D from above, the chemical
potential increases monotonically. Once the temperature reaches the transition temperature, µbecomes zero and remains zero for all lower temperatures. This behavior differs from that seen inlower dimensions as we shall see.
2

The Riemann Zeta function ζ(s) obeys gν (1) = ζ(ν)

72
3D uniform potential results

For the 3D uniform gas we have the following results:√
Z ∞V m3/2

z
N=√
d
+
1−z
2 π 3 ~3 0 e(−µ)/kB T
g3/2 (z)
Nex =
λ3T

2/3
2π~2 NT3D =ζ(3/2)
mkB V

Total number (with correction)(4.17)

Occupation of excited states(4.18)

Transition Temperature to BEC(4.19)µ(T ≤ T3D ) = 0

Chemical potential below T3D(4.20)µ(T  T3D ) = kB T ln(n3D λ3T ).

Chemical potential well above T3D(4.21)

4.2.3.2

2D uniform potential

In two dimensions the density of states D2U () = m/(2π~2 ) does not treat the ground stateas pathologically as the density of states in three dimensions D3U . Therefore, we shall not treat
the ground state separately and simply evaluate Eq. 4.10. The result isN =−V2 mkB T
ln(1 − z)
2π~2(4.22)

Solving for the chemical potential µ yieldsµ = kB T ln(1 − e−T2D /T ) < 0(4.23)

where
T2D =

2π~2 n2D
mkB(4.24)

is the characteristic temperature for the uniform 2D system, n2D = N/V2 is the number density
and V2 is the area (2D volume) of the system. The chemical potential thus increases monotonicallywith T . For T  T2D , the exponential in the logarithm approaches unity from below. Expandingln(1 − e−x ) about x = 0 and substituting x = T2D /T yields
" 




 #
T2D
1 T2D
1 T2D 2
µ = kB T ln
−
+
...
T
2
T
24
T

(4.25)

which for high temperatures reduces to simply the first term. So for T  T2D , µ ≈ −kB T ln(T /T2D )(this agrees with [98, 99]). For T  T2D the exponential et2D /T approaches zero from above.

73
Expanding ln(1 − y) about y = 0 and substituting y = e−T2D /T yields"#
e−T2D /Te−2T2D /T
e−3T2D /Tµ = −kB T
+
+
...
1
2
3

(4.26)

which reduces to simply the first term. Hence T  T2D implies the chemical potential approaches
zero from below as µ ≈ −kB T e−T2D /T . The occupation of the ground state can be approximated1
N0 ≈ kT /|µ|. Then for N0 ≈ N , the condition for BEC, requires that T /T2D ≈ ln(N) whichprecludes BEC of an ideal gas in the thermodynamic limit for two dimensional systems (except,
perhaps, in the theoretical case of exactly T = 0). Since for any finite temperature we do no have amacroscopic occupation of any of the low lying states, we are justified in our decision not to treat
the ground state separately as was required in three dimensions. We can express Eq. 4.23 in termsof the degeneracy parameter n2 λT2µ = kB T ln(1 − e−n2 λT ).
2D uniform potential results

AmkB T −|µ|/kB T
ln
1
−
e
2π~2
2π~2 n2D
T2D =
mkB


µ = kB T ln 1 − e−T2D /T(4.27)In 2D we have the following results:
Total number(4.28)

Transition temperature(4.29)

chemical potential(4.30)µ ≈ −kB T e−T2D /TT2Dµ ≈ kB T ln
T
T2D
n2D λ2T =Tµ = kB T ln 1 − e−n2D λT

T  T2D(4.31)T  T2D(4.32)

degeneracy parameter(4.33)

chemical potential(4.34)kB T≈ eT2D /T|µ|T1
N0 ≈ N =⇒
≈
T2D
ln N

ground state occupation(4.35)

Precludes BEC in thermodynamic limit(4.36)

N =−

N0 ≈

74
4.2.3.3

1D Uniform potential√

In one dimension the density of states diverges at zero momentum since D1U () = √2mπ~ √1 .
Therefore, we need not worry about the integral neglecting a macroscopic occupation of the groundstate if one should occur (which it will not for finite temperature). In one dimension the integral
in Eq. 4.10 becomes
1
N = V1
2

r

mkB T
L
g1/2 (e−|µ|/kB T ) =g (e−|µ|/kB T )2
2π~
2λT 1/2(4.37)

where V1 is the length of the system. Defining the 1D degeneracy parameter n1D λT and squaring
yields(n1D λT )2 =2π~2 n21D
T1D
=
mkB T
T

which motivates the identification
T1D =2π~2 n21D
mkB

as the characteristic temperature in 1D. For T  T1D (equivalently n1D λT  1) we can use thebehavior near zero of gν (z) ≈ z sor
2n1D λT = 2

T1D= g1/2 (z) ≈ eµ/kB T .
T

Therefore, The high temperature behavior for µ is given by!
r
T1Dµ = kB T ln 2
= kB T ln(2n1D λT ).
T(4.38)In the opposite limit, where n1D λT  1 we can make use of the limit as the argument of g1/2
approaches unity from below[130]lim gν (ew ) = Γ(1 − ν)(−w)s−1 .

|w|→0

So
µ/kB T

2n1D λT = g1/2 (e√ √−2 π kB Tp)≈|µ|

The low temperature behavior of the chemical potential is thus
µ≈−

πkB T
πkB T 2=−.
(n1D λT )2
T1D

75
The occupation of the ground state is thus
N0 ≈ kB T /|µ| ≈

T1D
.
πT

Again we see that large ground state occupation N0 ≈ N requires T /T1D ≈ 1/N
1D uniform potential results:
L
N=
2

r

In 1D uniform systems we have the following results:



mkB T
−|µ|/kB Tge
2π~2 1/2
2π~2 n21D
T1D =
mkB
r
T1D
n1D λT =T!
r
T1Dµ = kB T ln 2
= kB T ln(2n1D λT )T
πkB T
πkB T 2=−
(n1D λT )2
T1DT1
N0 ≈ N =⇒
≈
T1D
πN
µ=−

4.2.4

Total number(4.39)

Characteristic temperature(4.40)

Degeneracy parameter(4.41)T  T1D(4.42)T  T1D(4.43)

Precludes BEC in thermodynamic limit(4.44)

Harmonically trapped ideal systems

In uniform ideal systems, we’ve seen how the dimensionality precludes true BEC in thethermodynamic limit for one and two dimensional systems. This was done by determining the behavior of the chemical potential µ in the limit of highly degenerate systems where the d-dimensional
degeneracy parameter nλdT  1. For d = 1, 2, µ was a monotonically-increasing function for all temperatures, and so the fugacity z = eµ/kB T < 1. The occupation of the ground state N0 = z/(1 − z)then remains negligible for all finite temperatures. In this section we consider harmonically trappedsystems. We shall see that BEC is present in two and three dimensions in the thermodynamic limit
below a well defined transition temperature. The one dimensional system will be an edge case forwhich the transition temperature approaches zero in the thermodynamic limit, but BEC may still
occur in lab settings [76].

76
4.2.4.1

3D Harmonic confinement

The density of states for the three dimensional harmonically trapped gas is D3H () = 2 /(2~3 ωx ωy ωz )and so as in the uniform case, we must treat the ground state with  = 0 separately. The total
number of excited particles is:
Z ∞1
2
Nex = 3 3d
2~ ω̄ − e(−µ)/kB TkB T 3=g3 (z)~ω̄

(4.45)

The critical temperature occurs for z = 1 as in the uniform case and solving for the criticaltemperature yields
~ω̄
T3D =
kB



N
ζ(3)

1/3

(4.46)

We can use this temperature to rewrite the number of excited atoms as
Nex
=
N



kB T
~ω̄

3

ζ(3)
=
N



T
T3D

3

which is a very familiar looking expression. 3D harmonic trap results(4.47)We have the following

results for for a harmonically trapped ideal Bose gas in three dimensions.
Z ∞
1
2
d
2~3 ω̄ 3 − e(−µ)/kB TkB T 3=g3 (z)
~ω̄


~ω̄
N 1/3
T3D =
kB ζ(3)" #T 3
N0 = N 1 −T3D(4.48)

Nex =4.2.4.2

Number of excited atoms(4.49)

Transition temperature(4.50)

Condensate fraction(4.51)

2D Harmonic confinement

The density of states for the two dimensional harmonically trapped gas is D2H () = /(~2 ωx ωy ).
In the uniform 2D system we did not need to separate the ground state out of the integral, butsince the density of states erroneously gives the ground state a weight of zero for the trapped gas,

77
we must separate the ground state:
1
N = N0 + 2 2
~ ω̄

Z ∞



0

e(−µ)/kB Td.

(4.52)

This is suggestive that there will be BEC in two dimensions. The number of excited particles is
thus:

Nex =

kB T
~ω̄

2
g2 (z).

(4.53)

the two dimensional trapped system seems much more analogous to the three dimensional trappedsystem than was the case for uniform systems. The transition temperature is analogously derived:


~ω̄
N 1/2
T2D =
(4.54)kB ζ(2)The condensate fraction follows as in 3D.
2D harmonic trap results:
ideal Bose gas in two dimensions.
Z ∞1Nex = 2 2d(−µ)/kBT~ ω̄ − ekB T 2=g2 (z)
~ω̄


~ω̄
N 1/2
T2D =
kB ζ(2)π2ζ(2) =6" #T 2
N0 = N 1 −T2D
4.2.4.3

We have the following results for for a harmonically trapped(4.55)
Number of excited atoms(4.56)

Transition temperature(4.57)

ζ(2) unlike ζ(3) has a convenient (and famous) value(4.58)

Condensate fraction(4.59)

1D Harmonic confinement

The density of states for the one dimensional harmonically trapped gas is D1H () = 1/~ωand is thus simply a constant. Therefore the integral
Z ∞
1
1
Nex =
d
(−µ)/kBT~ω̄ 0 ediverges. Let us return to the sum
N=

X


1
e(−µ)/kB T − 178which is necessarily true. It turns out that one can proceed by separating the sum into threecomponents and exploiting some obscure mathematical identies. This is is the approach of [99].
An alternative derivation is given by [76] and yields



kB T
~ω
N − N0 = −ln 1 − z exp −
~ω
2kB T

(4.60)

This expression relies on the approximation ~ω  kB T . The transition temperature is determinedfrom the following transcendental equation
kB T1DN=ln~ω

2kB T1D~ω(4.61)

which has the solution in terms of the Lambert W function W0 (x). In the thermodynamic limit
N → ∞ with N ω → C for some constant C the transition temperature goes to zero according to
T1D =

~ωN
kB W0 (2N )(4.62)

since the denominator blows up but the numerator remains fixed. However, if the confining potential
is at all tighter than harmonic, so that V (x) ∝ kx2−ε for ε > 0 then T1D > 0 in the thermodynamic
limit [9]. Of course in finite harmonic systems BEC will still occur, an the condensate fraction goes
as
N0
=1−
N
T

T ln

1D ln



2kB T
~ω





2kB T1D~ω

(4.63)The ideal one dimensional gas in a harmonic potential thus will generally undergo BEC in laboratory
conditions.
1D harmonic trap results:

We have the following results for for a harmonically trapped

ideal Bose gas in one dimension.



kB T
~ω
Nex = −
ln 1 − z exp −
~ω
2kB T


kB T1D
2kB T1D
N=
ln
~ω
~ω
~ωN
T1D =
kB W0 (2N )

 

BT
T ln 2k~ω
N0 = N 1 −T1D ln 2kB~ωT1D

Number of excited atoms(4.64)

Equation for transition temperature(4.65)

Transition temperature using Lambert W(4.66)

Condensate fraction(4.67)

79

4.3

Weakly interacting Bose gases
The previous sections concerned ideal gases only. When interactions are included many

properties of the system may be changed dramatically. The theory of an almost ideal Bose gas
was first developed by Bogoliobov in 1947 and his approach remains popular today. Bogoliobov
theory for weakly interacting dilute Bose gases may be found in textbooks on BEC and superfluidity[102, 98] and we shall not go into as much detail here, but a thesis concerning the measurement offluctuations in interacting BEC would be incomplete without a presentation of the basic theoreticalframework. Bogoliobov theory provides the dispersion relation used in the Introduction of thisthesis to motivate the implosion technique used in this thesis. The prerequisites for applicabilityof Bogoliobov theory will be introduced in the context of interactions and their effects before
outlining the formalism of the Bogoliobov approximation and transformation. We shall conclude
with a presentation of the Gross-Pitaevskii equation.
4.3.1

Dilute and cold gases

In the context of Bose gases dilute means that the interparticle distance is much larger than
the interaction range of interparticle forces. A cautionary note is warranted; the s-wave scattering
length a which parameterizes low energy scattering is the characteristic length ultimately used to√√determine the diluteness condition na3  (15 π /128)2 ≈ 1/25 where (15 π /128) is the LHY√
coefficient to the quantity na3 in the beyond mean-field correction to the (many-body) groundstate energy [102]. However, the scattering length ultimately depends upon the actual interatomic
interactions which may be characterized by interaction range r0 . The diluteness condition is then
 
N
−1/3
r0  d = n=.
(4.68)
V
We generally treat the interaction as vanishing beyond r0 since the interparticle potentials U (r)generally tend towards zero for r  r0 . The Fourier transform of the potential U (r) yields thescattering amplitudeZ
U (p) =

U (r)e−ip·r/~ dr,(4.69)

80
from which we can see that for p  ~/r0 , where U (r) = 0, the scattering amplitude U (p) rolls
off since the exponent will oscillate the integral to zero. However in the opposite limit then the
exponent becomes unity and V (p) is independent of the momentum p. The range of the potential
prompts an obvious characteristic momentum pc = ~/r0 . Therefore, by cold we mean that the√thermal momentum pt ∼ 2mkb T is much smaller than pc . Hence, the coldness condition may be
expressed as
T 

~2
.
2mkb r02(4.70)

From this limitation on the momentum, we see that the low p behavior of Eq. 4.69 will dominateRand hence we can approximate U0 = U (r)dr with the result
U0 = g =

4π~2
a,
m

(4.71)awhere a is the s-wave scattering length we all know and love. The interaction term U0 = g = 4π~m2

may be derived using the Born approximation for the scattering amplitude between two particles,
details may be found in standard texts [102, 98]. We may express the Hamiltonian of the systemin terms of the field operators ψ̂:
Z 
Ĥ =


Z
~2
1
†
∇ψ̂ (r)∇ψ̂(r) dr +ψ̂ † (r)ψ̂ † (r0 )U (r0 − r)ψ̂(r0 )ψ̂(r)dr0 dr
2m2(4.72)

where U (r0 − r) is the two-body scattering potential. In the case of a uniform gas occupying a
volume V the field operator can be expanded as a superposition of plane waves :1 X
ψ̂(r) = √âp eip·r/~ ,V p(4.73)

where âp is the annihilation operator for a single particle state of a plane wave with momentum p
and the commutators for âp are [âp , âq ] = δpq and [â†p , âq ] = 0. Substituting of this ansatz into
Eq. 4.72 yields:
Ĥ =

X p2
p

Here Vq =R

2mâ†p âp +1 X
Vq â†p1 +q â†p2 −q âp1 âp2 .
2V p ,p ,q1(4.74)

2

U (r) exp[−iq · r/~]dr is the Fourier transform of the two-body scattering potential.

To get Eq. 4.74 from Eq. 4.72 one must use Eq. 4.73 with p → [p1 + q, p2 − q, p2 , p1 ] for81[Ψ† (r)Ψ† (r0 ), Ψ(r0 ), Ψ(r)] respectively in Eq. 4.72. Then move the q exponents next to the U (r0 −r)
to get the Vq term via Eq. 4.69. So the Hamiltonian after these substitutions and some rearranging
is expressed
X Z †


1
†
0
0
Ĥ = (...) +
â
â
U
(r
−
r)
exp
iq
·
(r
−
r)/~
p
+q
p
−q
1
2
2V 2 p ,p ,q12(4.75)

âp1 âp2 exp i(−p1 · r − p2 · r0 + p1 · r + p2 · r0 )/~ dr0 dr
The other exponentials cancel and since we have two spatial integrals, one of which is used up for
Vq and the other is then empty (all creation and annihilation operators are independent of r so canbe moved out of the integral). That second integral will also give us a needed factor of V .
In the name of coldness and diluteness we are permitted to simply replace the integral overU (r0 − r) with U0 . So we would simply replace Vq with U0 in Eq. 4.74 and only keep q = 0 terms.
This yields the new Hamiltonian:
Ĥ =

X p2
p

2m

â†p âp +U0 X † †â â âp âp .
2V p ,p p1 p2 1 2
1(4.76)

2

Remember, we are considering a BEC here, so the lowest energy state is occupied with a macroscopic
number of particles in the condensate. This means that one can use the Bogoliobov approximation
which is to replace the operator â0 with a c-number.
â0 ≡

p
√
N0 ≈ N

(4.77)In a cold and dilute BEC, the occupation number of particles in excited momentum states is very
low. Therefore in the lowest-order approximation, we keep only the p = 0 term and the Hamiltonian
yields the energy of the ground stateE0 =U0 2 U0
N ≈
Nn
2V 0
2

(4.78)

We can determine the pressure of the gas and its compressibility.
P ≡−

∂E0
g
= n2
∂V
2

∂n
1
=
∂P
gn

(4.79)

82
The hydrodynamic equation yields the speed of sound c:
r

1
∂n
=
2
mc
∂P

c=

gn
m

(4.80)

The compressibility is infinite for an ideal gas. In an interacting system the compressibility is
positive only when a > 0. In fact, the instability of a uniform Bose gas for attractive interactions is
manifest in these results; negative a gives negative pressure and an imaginary speed of sound, bothof which indicate such a system is not thermodynamically stable. Finally, the chemical potential
is given byµ=∂E0= gn = mc2∂N(4.81)

which is always positive when a > 0. This is the energy required to add one particle into thecondensate, or the energy gained by subtracting one particle from the condensate.
4.3.2

Higher order approximations: Bogoliobov approximation and excitation
spectrum

In computing the lowest order approximation to the ground state energy in Eq. 4.78 we only
considered the zero-momentum terms of the Hamiltonian 4.76. We are ultimately interested influctuations and excitations and so require a higher order approximation to the Hamiltonian. Let’s
examine the intermediate Hamiltonian Eq. 4.75 we used in the previous section
X Z


1
Ĥ = (...) +â†p1 +q â†p2 −q U (r0 − r) exp iq · (r0 − r)/~2
2V p ,p ,q12(4.75)

âp1 âp2 exp i(−p1 · r − p2 · r0 + p1 · r + p2 · r0 )/~ dr0 dr
Even though we go to higher orders in p here, p remains small. This means q is remaining smallenough that Vq = U0 no matter what. So this yields Eq. 4.76
Ĥ =

X p2
p

2m

â†p âp +U0 X †âp1 +q â†p2 −q âp1 âp2 .
2V p ,p ,q1(4.76)

2

Let’s consider the physical meaning of Eq. 4.76 when applied to our state. As written here,
we are annihilating a single p2 particle, then annihilating a single p1 particle, then creating a

83p2 − q particle, and lastly creating a p1 + q particle. But for that to be physically valid, we know
momentum must be conserved, so physically we can deduce certain restraints. Let us consider all
the possible ways of expressing that second term to the next nonvanishing order, which by our
analysis above, will be quadratic in p since linear terms would violate conservation of momentum.
Of course we still have the p1 , p2 , q = 0 term, which we will separate out. We are left with avariety of terms with nonvanishing pi . Suppose p1 and p2 are each zero. Then we can see thatwith q = p the term ∝ â†p â†−p â0 â0 (creation of two particles of opposite momenta) is non-vanishing.
This suggests there is a term corresponding to the complimentary physical process, the annihilationof two particles of opposite momenta, a term ∝ â†0 â†0 âp â−p , for which we have q = −p. There areno other quadratic terms where p1 and p2 have the same magnitude.
If p1 = 0 and p2 = p then our annihilation operators destroy a single particle with p and thecreation operators must fix to preserve momentum conservation. This can be done in two ways,
with q = 0 which gives the term ∝ â†0 â†p â0 âp and with q = p which gives the term ∝ â†p â†0 â0 âp .
Similarly is p1 = p and p2 = 0 then q = 0, −p yields terms ∝ â†p â†0 âp â0 and ∝ â†0 â†p âp â0 . But
the commutation relations allow use to write the four terms described in this paragraph together
in the form â†0 â†p â0 âp . Hence, to quadratic order in p our Hamiltonian may be written:
Ĥ =

X p2
p

2m

â†p âp +U0 † †â â â0 â0 +2V 0 0
U0 X  † †4â0 âp â0 âp + â†p â†−p â0 â0 + â†0 â†0 âp â−p2V

(4.82)

p6=0

Where we haven’t touched the first term. We continue with the Bogoliobov approximation for√â0 ≈ N in the third term of Eq. 4.82, but we need to be more precise in the second term. To
P
do this we use the normalization relation â†0 â0 + p6=0 â†p âp = N , which is true by inspection. We
can hit both sides of that to getâ†0 â†0 â0 â0 = â†0 N â0 − â†0 Xâ†p âp  â0 .

p6=0

Now we use the normalization relation to substitute again for â†0 â0 , which gives the N 2 term andPa single term of N . We use the commutation relations [âp , âq ] = δpq and [â†p , âq ] = 0 to get the84other term like N

P(and dropping the higher order terms). We are left withâ†0 â†0 â0 â0 = N 2 − 2NX

(4.83)â†p âpp6=0

We also require the scattering length a beyond the lowest order Born approximation. Without
deriving them here, we find the results [102]:


2
X
m 
U
m
a=
U0 − 0
2
4π~
2V
p2


U0 = g 1 +

p6=0


g Xm
V

p6=0

p2



(4.84)awhere as always g = 4π~
m . Doing the various substitutions outlined we get the following expression:2

Ĥ =gN 2 X p2 †
1 X †
mgn
+
âp âp + gn
2âp âp + â†p â†−p + âp â−p + 2 .
2V
2m
2
p
p

(4.85)

p6=0

4.3.2.1

Bogoliobov Transformation

We are ready to do our Bogoliobov transformation. Equation 4.85 can be diagonalized via
the Bogoliobov transformation:âp = up b̂p + v−p b̂†−p ,â†p = up b̂†p + v−p b̂−p .

(4.86)The two parameters up and v−p are determined by the following constraints: the new operators b̂†pand b̂p are assumed to obey the same Bosonic commutation relation as the real particles operatorsâp and â†p :hib̂p , b̂†p0 = δpp0

(4.87)

This commutator definition thus constrains the two parameters up and v−p . Some math follows
which is not particularly physically illuminating, see section 4.2 of [102] for the details. The explicit
form of up , v−p is 2
1/2
p /2m + gn 1
up =
+
(p)
2

 2
1/2
p /2m + gn 1
v−p = −
−
(p)
2

where
"

gn 2
(p) =
p +
m



p2
2m

2 #1/2(4.88)(4.89)

85
is the Bogoliobov dispersion law for the elementary excitations of the system. We can now diagonalize the Hamiltonian 4.85 into the simple form
Ĥ = E0 +X

(p)b̂†p b̂p ,(4.90)

where


1
1X
p2
m(gn)2E0 = gnN +(p) − gn −
+
,
2
2
2m
p2

(4.91)

p6=0

is the ground state energy calculated to the higher-order approximation. Looking at the Hamiltonian 4.90 we see it resembles that of a system of non-interacting quasiparticles each with energy
(p). Consider how the dispersion relation depends on p; in the limit of small and large p we haver
p→0 ≈ p

gn
= pc,
m

p→∞ ≈

p2
+ gn.
2m

(4.92)

Looking at the explicit form of u, v substituting the low momentum dispersion relation will yield
up→0 , vp→0 ≈ ±

mc 1
±
p
2

1/2

(4.93)


so |up | ≈ |v−p |  1 and âp ∼ up b̂p − b̂†−p . At low momenta, the real particles resemble aroughly equal number of quasiparticles of positive and negative momenta. For the higher momenta
limit, we know that  ∼ p2 /2m + gn and so |up | ' 1 and |v−p | ' 0 and the quasiparticle b̂p is
indistinguishable from the real particle so âp ∼ b̂p . The fact that the high momentum excitations
resemble real particles is also suggested by the high momentum behavior of Eq. 4.92 which resembles
that of a free particle. The ground state of the weakly interacting system at T = 0 is the vacuum
state of the Bogoliobov quasiparticles :b̂p6=0 |0i = 0(4.94)The transition between the phonon-like and particle like regimes occurs around p2 /2m ∼ gn = mc2 .
Setting p2 /2m = gn and letting p = ~/ξ defines the healing length:
s
~2
1 ~
ξ=
=√
.
2mgn
2 mc

(4.95)

86
The occupation number Np of the quasi-particles is given by the normal Bose expression
Np = hb̂†p b̂p i =1
exp[β(p)] − 1(4.96)

This must not be confused with the average particle occupation number hâ†p âp i. The chemical
potential of a gas of quasi-particles is zero by construction since their number is not fixed but
is determined by the condition of thermodynamic equilibrium, in analogy with Planck’s black
body law for photons. The actual particle occupation number hâ†p âp i is, under the Bogoliobov
transformation,
np = hâ†p âp i = |v−p |2 + |up |2 hb̂†p b̂p i + |v−p |2 hb̂†−p b̂−p i .

(4.97)

holding only for p 6= 0. This is important to note since if we are interested in the number of atoms
in the condensate (for which p = 0) we have
N0 = N −X
p6=0

V
np = N −(2π~)3Z 

|up |2 + |v−p |2|v−p | +d3 pexp[β(p)] − 1
2(4.98)We see how the |v−p |2 term indicates that the interactions cause there to be some particles with
non-zero momenta, even in the limit of T = 0 sincelim

T →0,p6=0

(p)/(kB T ) = ∞ =⇒ hb̂†p b̂p i = 0

which is obvious from the form of Eq. 4.89. Here, we see the difference between the quasi-particles
occupation number with momentum p 6= 0 and the real particles occupation number for such p.
Even though we calculate zero quasiparticles at absolute zero, we still have a non-zero number ofactual atoms in higher p states.
4.3.3

LHY corrections

The ground state energy Eq. 4.91 in the higher-order approximation can be easily calculated
by replacing the sum with an integral in momentum space resulting in


gnN
128
3 1/2
E0 =
1 + √ (na )
,
2
15 π

(4.99)

87
while the chemical potential µ = ∂E0 /∂N is thus easily shown to be



32
3 1/2
µ = gn 1 + √ (na ).
3 π(4.100)

These expressions give the LHY correction to the mean-field values of the ground state energy andchemical potential for a weakly interacting Bose gas. They are included for completion and are not
relevant for the work in this thesis.
4.3.4

The GP equation

The dynamics of the weakly interacting Bose gas is well described by the Gross-Piteavskii
equation, which suitable not only for uniform systems for which the Bogoliobov theory above wasdeveloped but also for non-uniform systems. The GP equation describes the evolution of the
wavefunction Ψ(r, t) according to


~2 2
∂
2
−
∇ + Vext (r, t) + g|Ψ(r, t)| Ψ(r, t) = i~ Ψ(r, t).
2m∂t
The mean field term, g =R(4.101)

Vef f (r)d3 r can be expressed in terms of the scattering length yieldingg = 4π~2 a/m. This equation resembles the Schrodinger equation with an additional term describing
the interactions which renders the GP equation nonlinear. We are interested in fluctuations fromthe mean and so it is natural to expand the wavefunction asΨ(r, t) = [Ψ0 (r) + O(r, t)] e−iµt/~(4.102)where O(r, t) describes small oscillations about the mean equilibrium value Ψ0 (r). ExpandingO(r, t) as
O(r, t) =X[ui (r)e−iωi t + vi∗ (r)e+iωi t ](4.103)iand substituting into the GP equation leads to the coupled differential equations~ωi ui (r) = (Ĥ0 − µ + 2gn(r))ui (r) + g(Ψ0 (r))2 vi (r),−~ωi vi (r) = (Ĥ0 − µ + 2gn(r))vi (r) + g(Ψ∗0 (r))2 ui (r),(4.104)

88∇
where Ĥ0 = − ~2m+ Vext (r). So far these representations are entirely classical, reflecting the fact
2

2that the condensate behaves as a classical field. Quantization is achieved by rewriting these interms of field operators consisting of the still-classical condensate order parameter Ψ0 (r) and an
operator describing the fluctuations Ô(r, t)hi
Ψ̂(r, t) = Ψ0 (r) + Ô(r, t) e−iµt/~ .

(4.105)

The fluctuation operator is defined analogously to Eq. 4.86 asÔ(r, t) =X[ui (r)b̂†i e−iωi t + vi∗ (r)b̂i e+iωi t ](4.106)iwhere b̂†i and b̂i are the creation and annihilation operators for the i-th elementary excitation
which has frequency ωi . The ui , vi and ωi are the solutions to the Bogoliobov equations of Eqs.
4.104. In this way, the Bogoliobov theory of the uniform case is extended to the trapped case. My
understanding is that Eli solves for those things.

4.4

BEC, Superfluidity, and Fluctuations
In the first section of this chapter the criterion for BEC was determined to be a macroscopic

occupation of one of the low-lying single-particle ground states, and this criterion was shown tobe equivalent to the behavior of the off-diagonal behavior of the one-body density correlation
function. We then developed the theory of BEC in ideal gases for dimensions d = 1, 2, 3 forboth uniform systems and harmonically confined systems. The role of dimensionality and theexternal potential was shown to be crucial to the existence, or lack thereof, of BEC at finite orzero temperature. In determining whether or not those systems exhibited BEC we emphasized themacroscopic occupation criterion, but did not explore the existence of BEC from the perspective of
the one-body density matrix. The role of fluctuations has thus been relegated to the backgroundso far in this chapter. Moving forward, we wish to emphasize the role of fluctuations. The physics
of BEC and fluctuations thereof more or less requires Bogoliobov’s theory and its extension to
non-uniform systems and so was introduced in Section 4.3.

89With the vocabulary and formalism established for these systems, we may now focus of therole of fluctuations in BEC. The interplay between the phenomena of Bose-Einstein condensationon the one hand and superfluidity on the other may be understood in the framework of quantum
hydrodynamics. The role of the wavefunction’s phase turns out to be central to this interplay, more
central than its modulus. In the following section we highlight the role of phase and dimensionality
in the behavior of the off diagonal portion of the one-body density matrix, providing a fluctuation
focused perspective on the phenomenon of BEC to compliment the energetics based arguments ofthis chapter’s first section. This framework provides, at least in my opinion, a view which is more
physically intuitive in understanding how fluctuations manifest in d-dimensional system. We shallbegin by discussing a criterion for superfluidity in terms of the spontaneous creation of elementary
excitations. We then highlight the importance of the phase and describe BEC in the language ofquantum hydrodynamics. The nature of phase fluctuations in determining the presence or absenceof long-range order in various dimensions will finally be discussed.
4.4.1

Landaus criterion for superfluidity

Consider a (uniform) fluid with energy E and momentum P in a given reference frame K andsuppose K 0 is a reference frame moving with velocity V with respect to K. If we ignore relativity
and proceed with Galilean transformations we may express the energy of the fluid in K 0 by
1
E0 = E + M V 2 − P · V
2(4.107)0

P = P − MV
where M is the total fluid mass. Suppose one of the frames (K) is the fluid’s rest frame and consider
the creation of an excitation with momentum p and energy (p) within the fluid which is initiallyin the ground state with energy E0 . The total energy of the fluid is then E0 + (p) and the total
momentum is P = p by construction. Let’s give our fluid some context by further supposing thefluid was initially flowing within a tube at constant velocity v and denote the tubes rest frame as

90
K 0 . Then Eq.4.107 indicates:
1
E 0 = E0 + (p) + M v 2 − p · v
2(4.108)0

P = p − Mv
In the tube’s rest frame, the energy of the system has changed by ∆E = (p) − p · v. Therefore we
identify (p) − p · v as the energy of the excitation in the tube frame and conclude that excitations
may appear if
(p) − p · v < 0,

(4.109)and the fluid may transfer momentum to the tube only under the necessary condition v > (p)/|p|.
The creation of such an excitation and the subsequent transfer of the momentum to the environment
represents a dissipative mechanism. Since the system’s excitations follow the dispersion relation
(p), there will be a momentum p minimizing (p)/|p|. We may therefore define a critical velocity
vc of the fluid through the tube for which
vc = min
p

(p)
.
|p|

(4.110)If the flow through the tube is less than vc then it is impossible for the creation of an elementary excitation to be energetically favorable and the flow will not dissipate. The criterion for superfluidity,
which is defined as flow without dissipation, is thus:v < vc = min
p

(p)
.
|p|

(4.111)This is Landau’s criterion of superfluidity [79]. In the case of an ideal Bose gas, we have(p) =

p22m

and so vc = 0 indicating that the idea Bose gas (in an infinite uniform system) is not actuallya superfluid! At a finite but small temperatures the system will necessarily occupy some of theelementary excitations. In the Bogoliobov picture we may view these excitations as non-interacting
quasiparticles and assume the system behaves like a non-interacting gas of these quasiparticles.
Thermalization with the environment takes place via interactions between the quasiparticles and91
the wall of the tube, and so excitations may be created and destroyed through this process. The
thermally excited excitations of course contribute to the flow of particles through the tube, but
in their interaction with the tube, the flow associated with these excitations does dissipate! The
finite temperature fluid thus possesses a superfluid component which flows without dissipation
and a normal component which does dissipate, at least with respect to the walls of the tube. The
normal component will therefore equilibrate such that its velocity vn is zero with respect to thetube wall. The occupation of the elementary excitations is thus
Np =1
exp[((p) + p · (vs − vn ))/kB T ] − 1(4.112)where vs is the superfluid velocity and vn is the normal-fluid velocity.
4.4.2

The role of the condensate phase

One striking result of the previous section was the implication that a BEC in an ideal gas is
not a superfluid. Therefore, the identification of the condensate density |Ψ0 | with the superfluid
density ns is not a general truth. Of course, BEC and superfluidity are still related and we shedsome light on the nature of that relationship in this section. Under a Galilean transformation thefield operator Ψ̂(r, t) transforms according to:
 

i
1
2
Ψ̂ (r, t) = Ψ̂ (r − vt, t) exp
mv · r − mv t
~
2
0

0

(4.113)

In a uniform fluid we may write the condensate part of the wavefunction in the coordinate system
where the sample is in equilibrium [102] asΨ0 =√n0 e−iµt/~ .

(4.114)In the frame where the fluid is moving with velocity vs the order parameter may be written√

n0 eiΦ 
11
Φ(r, t) =mv · r −mv 2 + µ t
~
2
Ψ0 =

(4.115)

92
where Φ is the phase in the moving frame. From Eq. 4.115 we extract the important result
vs =

~
∇Φ
m

(4.116)

relating the superfluid velocity to the phase of the condensate order parameter. The appearanceof the order parameter’s phase Φ but not its density n0 is indicative of the more important role ofthe condensate phase to superfluidity compared to the density.
4.4.3

Quantum Hydrodynamics and the phase correlation function

The formalism of hydrodynamics is a useful framework for understanding supfluidity. For an
normal fluid we have the total energy of the system H:
Z
H=

(∇φ)2

ρ
+ e(ρ)dr,2(4.117)

where φ is the so-called velocity potential obeying vs = ∇φ, ρ is the density of the fluid ande(ρ) is the internal energy per unit volume. The variation δH of the fluid energy appears in thehydrodynamic equations:
∂ρ
δH
=
,
∂t
δφ

∂φ
δH
=−
.
∂t
δρ

(4.118)Thus ρ and φ are the congucate variables of the system and the corresponding quantum mechanicaloperators obey the commutation relation[φ̂(r), ρ̂(r0 )] = −i~δ(r − r0 ).

(4.119)We may write the Hamiltonian for a quantum fluid as
Z
Ĥ =

ρ̂
∇φ̂ ∇φ̂ + e(ρ̂)dr,r(4.120)~where φ̂ = m
Φ̂. In terms of the operators b̂† and b̂ for the phonons which obey the commutation

relation b̂k0 b̂†k − b̂†k0 b̂k = δkk0 we may write the density deviation and phase operators as
1 X ρ̄~k 1/2  ik·r
δ ρ̂ = √b̂k e+ b̂†k e−ik·rc2V k(4.121)

93


1/2 

µt
1 X
m2 c
Φ̂ = − − √ib̂k eik·r − b̂†k e−ik·r(4.122)~ρ~k
2V k 2 1/2c
There is a typo in Equation (6.66) of [102], where the m in the mquantity is erroneously a
ρ~k 1/2mc
linear power (giving ρ~k). The form here is correct. The one-body density matrix n(1) (r, r0 ) =hΨ̂† (r)Ψ̂(r0 )i may be computed in terms of the field operator
Ψ̂(r, t) =√n0 eiΦ̂ ,(4.123)where n0 is the condensate density and Φ̂ is the phase operator. The one-body density matrix
obeys
D
E
1
0
0 2
n(1) (r, r0 ) = n0 e−i[Φ̂(r)−Φ̂(r )] = n0 e− 2 h[Φ̂(r)−Φ̂(r )] i = n0 e−(χ(0)−χ(s)) ,(4.124)whereχ(s) = hΦ(r)Φ(r0 )i(4.125)is the phase correlation function and s = |r − r0 | [74, 102]. This is a key relationship, as it directly
relates the long range order to the phase properties of the system. If regions of the system are highly
correlated, i.e. if the phase coherence of the system is large, then the phase correlation function
χ(0) ≈ χ(|r − r0 |) and the off diagonal part of n(1) (r, r0 ) will not drop off. We demonstrated in Sec.
4.1 that this is defines a BEC. Using Eq. 4.122 we have:χ(s) = hΦ(r)Φ(r0 )i*+
m2 c X −1  ik·r
00
† −ik·r†χ(s) = −kb̂k e− b̂k eb̂k eik·r − b̂k e−ik·r
2V ρ~
k
*
+

2
m c X −1  † ik·(r−r0 )†−ik·(r−r0 )χ(s) =
kb̂k b̂k e+ b̂k b̂k e
2V ρ~
k
*
+


2
X
m c
†
†
−1
ik·(r−r0 )−ik·(r−r0 )χ(s) =k(1 + b̂k b̂k )e+ b̂k b̂k e
2V ρ~
k

 ik·(r−r0 )m2 c X 1e†χ(s) =+ hb̂k b̂k i
V ρ~
2
k

(4.126)k

In the above, we have only one sum over k since all momentum cross terms vanish. We havealso dropped the terms containing b̂k b̂k and b̂†k b̂†k as they vanish (or contribute to higher order?).
From the fourth line to the last we have switched the sum over k in the last term (as is allowed by

94
symmetry) to produce like exponentials. Expressing the integral in terms of momentum p = ~kwe havem2 c
χ(s) =
ρZ 
1 eip·s/~ dd pNp +2p (2π~)d

(4.127)

where Np = hb̂†p b̂p i = [ecp/kB T − 1]−1 is the thermal distribution of phonons. In three dimensions
this is identical to Eq. (6.71) of [102]. In the high temperature limit kB T  cp we have Np ≈kB T /cp  1/2 and som2 kB Tχ(s) =
ρ

eip·s/~ dd p
p2 (2π~)dZ(4.128)

Expanding the exponential according to Euler’s formula and dropping the sine term since it is an
odd function yields:
m2 kB Tχ(s) =
ρcos (p· s/~) dd p.
p2(2π~)dZ(4.129)The density matrix thus decays according to n(1) (s) = n0 e−(χ(0)−χ(s)) with
χ(0) − χ(s) =m2 kB TρZ

1 − cos (p· s/~) dd p.
p2(2π~)d

(4.130)We can evaluate the integral for d = 1, 2, 3 finding:
m2 kB Th[Φ̂(r) − Φ̂(r )] i =ρZ ∞m2 kB T
ρ

Z pco

0 2

1 − cos (p· s/~) dp
p2
(2π~)
−∞

=

m2 kB T|r − r0 |2ρs ~2

1 Dimension(4.131)h[Φ̂(r) − Φ̂(r0 )]2 i =m2 kB Th[Φ̂(r) − Φ̂(r )] i =
ρ
0 2

0

Z pco
0

1 − cos (p· s/~) (2π)dpp(2π~)2≈m2 kB T 
0
γ
+
ln
2p
|r
−
r
|
co
2πρs ~2(4π)dp(1 − cos (p· s/~))(2π~)3

2 dimensions(4.132)
p |r−r0 |
co

~ sin
~
kB T 
≈ 2
pco −2
2π ρs ~
2|r − r0 |


(4.133)where γ ≈ 0.577 is Euler’s and constant we have introduced an ultraviolet cutoff pco because thehydrodynamic formalism is only valid for long wavelengths, and the integrals in d = 2 and d = 3
dimensions do not converge otherwise. We omit the cutoff in the one dimensional integral because
it does converge and otherwise the result is ugly, but strictly speaking, it too should be limited in

3 dimensions

95
its range. The 2π and 4π are simply the angular parts for two and three dimensions. The resultsof Eqs. 4.131, 4.132 and 4.133 were first given by Kane and Kadanoff, although the 2D and 3D
results differ slightly in form [74]. Other sources also have slight differences [98, 60] in 2D and 3D,
but they do not affect the limiting behavior. In 1D the one-body density matrix tends towards zero1
exponentially: n(1) (r − r0 ) ∼ e−|r−r | . In 2D it decays algebraically: n(1) (r − r0 ) ∼ |r−r
0 | . In 3D, it
0settles to a constant non-zero value. Therefore, in 1D and 2D we do not expect BEC in uniform
systems. Note that quantum hydrodynamics requires an interacting fluid, so these results do not
map perfectly onto the demonstration in Sec. 4.2.3, which were ideal systems. Of course, theseresults do not apply to T = 0 exactly, so we should not place too much confidence in them for zerotemperature systems.
Let us also compute the density-density correlation function hδ ρ̂(r)δ ρ̂(r0 )i. The one-bodydensity matrix was given in Eq. 4.124 only in terms of the phase correlation function, because thecontribution of the density is negligible in comparison. To understand why, we compute analogously
to Eq. 4.126. The operator math is identical so we can simply write down the answer%(s) = hδ ρ̂(r)δ ρ̂(r0 )i*+ρ̄~ X  ik·r
† −ik·r† −ik·r0
ik·r0=
k b̂k e+ b̂k eb̂k e
+ b̂k e2V ckρ̄~ X 10†=
+ hb̂k b̂k i keik·(r−r )Vc2(4.134)k

For kB T  pc we haveρ̄~%(s) =cZ
cos (p· s/~)dd p(2π~)d

(4.135)for the density correlations in d dimensions.

4.5

Discussion
The low momentum behavior for the density correlation function %(s) = hδ ρ̂(r)δ ρ̂(r0 )i andthe phase correlation function χ(s) = hΦ(r)Φ(r0 )i are very different. The phase correlations go as
k −2 in the integrand while k drops out of the integral entirely; although naturally the integrals

96will introduce a factor of k d−1 due to the volume element, and it is this which differentiates the
behavior of the one-body density matrix in different dimensions. The stronger dependence on thephase Φ̂ of the order parameter justifies the omission of density effects in Eq. 4.124.
The quantum hydrodynamic picture has provided us with quantitative results on the natureof fluctuations in multiple dimensions. We have seen how the phase fluctuations destroy long rangeorder in 1D and 2D for finite temperature, a complimentary view to arguments from the chemical
potential for ideal gases in earlier sections. We have also seen how the density fluctuations are muchless important in determining the off-diagonal long range behavior of the system. We have also seen
how in low dimensions the fluctuations of both the phase and the density are greater than in higher
dimensions. We have even seen some interesting intermediate behavior in the two-dimensional casethat has not yet been explored in our discussion, the algebraic, as opposed to exponential (as in
1D), decay of the one-body density matrix to zero. This suggests a sort of intermediate behavior
for these systems which turns out to be related to a transition to a superfluid state which does notexhibit BEC: the Berezinskii-Kosterlitz-Thouless (BKT) transition. However, the conceptual goalof this chapter, which was to explain how and why dimensionality affects fluctuations, has been
fairly thoroughly addressed.

4.6

The Static Structure Factor
With the conceptual understanding of fluctuations having been extensively explored so far,

it is time to shift gears and bring the theory back to the lab.
4.6.1

Static Structure factor for a > 0

The static structure factor is the quantity of primary interest in this thesis, and here we
provide its definition. The static structure factor is related to the dynamic structure factor S(p, ω),which characterizes inelastic scattering by a probe which transfers momentum ~p and energy ~ω tothe system [102]. One approach to understanding the static structure factor is from the perspectiveof linear response theory [102]. However, the conceptual task mentioned at the beginning of this

97
chapter was to develop the tools necessary to understand the nature of fluctuations in different
dimensions in order to motivate the technique used in this thesis, and it is this author’s hope thatthat task has been achieved in the previous sections. As such, we shall narrow our focus from thispoint forward to the question of what measurements will reveal the fluctuations in our bulk three
dimensional BEC, and what we expect to measure.
The treatment here is due to our theory collaborator Eli Halperin, and his forthcoming thesis
will be a much more thorough development of these ideas [62]. We begin with the correlation
function.
ν(r1 − r2 ) =1hδn(r1 )δn(r2 )in̄(4.136)

where δn(r) = n̄ − n(r) is the density deviation from the average. The static structure factor S(k)is the Fourier transform of the static correlation function
Z
S(k) =

ν(r)e−ik·r dr(4.137)

S(k) =

1h|δn(k)|2 i
N(4.138)δn(r)e−ik·r dr(4.139)We may rearrange to express this aswhereZ
δn(k) =

is the Fourier transform of the density deviation. The static structure factor is the ensemble averageof the square of the Fourier transform of the density deviation. In the Bogoliobov approximation
the static structure factor may be written
S(k) = hâ†k âk i + hâ−k â†−k i + hâ†k â†−k i + hâ−k âk i(4.140)where the âk are the annihilation operators for real particles with wavenumber k. We perform the
Bogoliobov transformationâk = uk b̂k + v−k b̂†−kâ†k = uk b̂†k + v−k b̂−k(4.141)

98
where
s
uk , v−k = ±and~2 k 2 /2m + gn 1
±
2(k)
2

(4.142)

 2 2 2~ k.
2m

(4.143)s
(k) =~2 k 22gn+
2m

The structure factor may now be expressed as
S(k) =


1 ~2 k 2  †hb̂k b̂k i + hb̂−k b̂†−k i + hb̂†k b̂†−k i + hb̂−k b̂k i
(k) 2m(4.144)

We proceed as usual withhb̂†k b̂k i = hb̂†−k b̂−k i =1e(k)/kB T(4.145)hb̂†k b̂†−k i = hb̂k b̂−k i = 0(4.146)[b̂k , b̂†k ] = 1(4.147)

which reduces to
1 ~2 k 2
S(k) =coth(k) 2m



(k))
2kB T


.

(4.148)This is the structure factor for a uniform BEC. The use of the Bogoliobov approximation â0 =√â†0 ≈ N in deriving Eq. 4.140 assumes a condensate contains most of the atoms.
4.6.2

Evolution of structure factor at a < 0

The details of the computation for the time evolution of the structure factor following a quench
to a < 0 is beyond the scope of this thesis, which, despite the topics covered in the last thirty orso pages, is an experimentally focused work. Details may be found in the forthcoming thesis of
our colleague Eli Halperin who performed the numerical simulations and theoretical work much ofthis thesis relies on [62]. The validity of Bogoliobov theory for systems with attractive interactions
is not obvious a priori. It turns out that the Bogoliobov procedure is incapable of diagonalizing
the Hamiltonian and producing Bosonic quasiparticles for all values of k [62]. Progress can still be99
made within the Bogoliobov framework and the evolution of the structure factor may be expressed
as





~2 k 2 n(gi + |gf |)f (k)t
1 ~2 k 2
(k)
coth
1+
sin2
(k) 2m
2kB T
m|f (k)|
~



2
2
~ k n(gi + |gf |)f (k)tS(k, t) = S(k, 0) 1 +sin2(gf < 0)2
m|f (k)|~S(k, t) =(4.149)This is the prediction for the evolution of the structure factor following a quench to attractive
interactions. The structure factor may be measured with in-situ density measurements of thecloud following a quench.

4.7

Related work in 2D
In two dimensions, fluctuations have been the subject of extensive experimental efforts. The

group of Cheng Chin as the University of Chicago and more recently the group of Chen-Lung Hung3
at Purdue University have explored interacting BECs in 2D and measured the static structure factor
in a number of experiments over more than a decade [72, 59, 71, 104, 29, 30]. In particular the workin [72] concerned the direct measurement of the structure factor in 2D for interacting condensatesat various interactions strengths, all with a > 0. They place great emphasis on a comprehensive
understanding of the imaging systems optical transfer function to accurately extract the density
correlations from in-situ images of their 2D BEC. They used images of a non-interacting thermal
cloud for which the expected structure factor is flat to characterize the aberrations in their system.
This was the motivation for analyzing the structure factor of the thermal component away fromthe condensate in Section 3.6 and comparing the resulting profile with the corresponding analysis
we use to determine the structure factor of the condensate.
In a later work, the oscillations predicted by Eq. 4.149 were observed following a quench
between two positive values of a [71]. In the experiment most analogous to our own, the cloud
was quenched from gi = 0.25 to gf = 0.079 where the interaction strength for the 2D superfluid√is g2D = 8π a/lz where a is the scattering length and lz is the harmonic oscillator length in thez-direction which is frozen out the dynamics to make the system effectively 2D. Growth of the3

Dr. Hung was a graduate student in Dr. Chin’s lab

100
structure factor was observed immediately following the quench, which was performed by abruptly
changing the magnetic field over 300 µs, which was sufficiently slow compared to all other relevanttimescales. One should note that the form of Eq. 4.149 presumes the gf < 0, we may rewrite theevolution for the structure factor following a quench between two values of a as
~2 k 2 n(gi − gf )2 f (k)tS(k, t) = S(k, 0) 1 +sin.
mf (k)2~(4.150)This expression is valid for all quenches, and correctly predicts no dynamics for gi = gf , whichEq, 4.149 as written, does not as it presumes gf < 0. The oscillations of of the structure factor
after the quench-down gi > gf measured in [71] matched the predictions of Eq. 4.150. In our case,for which gf < 0 (and g = g3D = 4π~2 a/m) the (k) may be imaginary and since sin(ix) = sinh(x),we do not expect oscillations for unstable mode but rather runaway growth. In the quench up
experiments of [71], for which gf > gi > 0, the initial trajectory of the structure factor for a given k
is downward, and the observed oscillation frequencies were 65% greater than predicted. Subsequent
theoretical work [104] considered the influence on a thermal bath of atoms on the evolution of thestatic structure factor. Most recently the Hung group at Purdue quenched to attractive interactions
and observed the growth of unstable modes followed by the fragmentation of the gas into 2D solitonsin [29]. In a subsequent paper, the interactions were quenched to attractive interactions and held
for up to 2 ms, and the growth of the structure factor was observed for stable and unstable modes,
as predicted by Eq. 4.150, see Fig. 2 of [30]. However, the amplitude of the growth was lower thanthat predicted by Eq. 4.150, and they fit their data toS(k, δτ ) = e−Γk ∆τ Scoh (k) + Sinc (k)(4.151)

where Scoh is given by Eq. 4.150 and
" "
#
#







2
Γ2k
1
∆τ
∆τ
f
Sinc (k) =η− 2 2
1 − e−Γk ∆τ cos 2f−
sin 2f+ η+ 1 − e−Γk ∆τ
2
~
~Γk
~
~ (Γk + 42f )(4.152)
with η± = 1 ± 2i /2f . Sinc (k) is an incoherent additive background which does not grow over time4 .
The coupling between the the coherent and incoherent contribution is governed by the k-dependent4

It is also the most annoying equation I’ve ever typeset in my life

101
damping rate Γk , which is a fit parameter of their model for which they find Γk ∼ 0.5i /~. They alsoemphasize the used of a high temperature thermal gas to calibrate their imaging transfer function.

Chapter 5
Controlled Demolition

5.1

Optical Lattice overview
The use of optical lattices to provide a well understood control case for implosions is an

important stepping stone to understanding how intrinsic fluctuations are amplified. This chapter
shall describe the geometry of the optical lattices used to write our imposed density perturbations,
how they are calibrated, and how those density perturbations evolve following a quench to negative
scattering length along with a comparison to numerical simulations of such a quench. We shallfinish with a calibration of our resolution at the two wavevectors of the X-lattice and Y-lattice.

5.2

Lattice geometries
The optical lattices used as seed potentials in this work are formed by the interference of

two running beams intersecting at shallow angles. Each lattice is referred to by the direction of its
wavevector. Only the X-lattice and Y-lattice are used in this thesis, so discussion of the Z-lattice
is omitted.
5.2.1

X-lattice geometry

The X-lattice geometry is shown in Fig. 5.1a. The X-lattice beams flank the V-beam of the
ODT and enter the science cell from above; the +125 mm lens which focuses the VODT beam steers
the X-lattice beams toward the atoms. Due to the small waist, and large corresponding divergence,of the X-lattice beams, each beam’s individual focus is longitudinally shifted down-beam from the103(a) X-lattice(b) Y-lattice

Figure 5.1: Geometry of X-lattice (a) and Y-lattice (b). Prior to the two-inch mirror which focuses
the V-beam of the ODT, the X-lattice beams are parallel to the V-beam (not shown). The V ODTlens steers the X-lattice beams inward such they intersect at the atoms with a half-angle of ∼ 8.5
degrees. The Y-lattice beams enter on either side of the H2-beam (not shown). The Y-lattice
beams intersect at ∼12 degrees.

atoms’ position and thus the lattice beams do not provide a noticeable perturbation to the trapping
potential for typical intensities used for seeding the density perturbation. The maximum attainable
power of the X-lattice beams is sufficient to shift the cloud’s equilibrium position, but this power
is only used in the trap-frequency measurements described in 2.6.1. For the typical purpose of
imposing a small seed potential the X-lattice beams do not perturb the cloud.
5.2.2

Y-lattice geometry

The Y-lattice beams enter the science cell on either side of the H2 ODT beam which provides
the vertical confinement for the atoms. These beams, like the H2 beam, are directed to the sciencecell by mirrors approximately 45 cm away from the science cell to provide room for the cart whichtransfers the atoms from MOT2 into the science cell. This severely constrains the space available for
the Y-lattice periscopes. As a result, proper Y-lattice alignment is much more difficult to achieve

104
since alignment must be performed using hex keys socketed into the adjustment knobs of the 1/2inch kinematic mirror mounts atop the Y-lattice periscopes. This is in contrast to the alignment
optics of the X-lattice which possess differential micrometers. Fortunately, the alignment of eachlattices is quite stable, so once the optimum alignment is achieved the lattice intensity does notvary day to day. As mentioned in Chapter 2, the Y-lattice light is fiber coupled and routed from theupper optical breadboard of Figure 2.10 to the section of the optical table beneath the H1/H2 andside-imaging breadboards. The fiber limits the maximum power of the Y-lattice to substantiallybelow that of the X-lattice. Unlike the X-lattice, the beams are not focused by any lens near thescience cell, and this further reduces the maximum Y-lattice depth compared to the X-lattice, thebeams of which are somewhat focused by the VODT lens, although the focus is down-beam of theatoms.

5.3

Optical lattice characterization
This section shall describe how the lattice’s are characterized. We shall first describe thephysics used for absolute measurements of the lattice depth using a short, powerful lattice pulse.
In this thesis, a short laser pulse means 0.5 − 500µs timescale control of a CW laser as controlled
by an AOM not the vastly faster timescales of pulsed lasers. A lattice pulse is simply both lattice
beams being pulsed at the same time. The lattice pulse calibration method will be shown to beinvalid for weak lattice of the sort used as seed potentials in this work. We shall then describe how
the shallow-lattice regime is calibrated with respect to the absolute lattice depth measurements ofdeep lattices.
5.3.1

Absolute lattice depth measurements of deep lattices

Consider a uniform BEC with uniform phase at t = 0 subject to a lattice potential for theduration of a short laser pulse. We write the potential as
U (x) = −U0
(1 − cos(2kL x)).
2

(5.1)

105
If pulse is on for a short enough time that none of the atoms will have moved much, then each portionof the wavefunction will acquire a phase given by U (x) in the normal manner ψ(t) = ψ(0)eiEt/~[110, 109]. If we suppose ψ(x, t = 0) = ψ0 thenψ(x, t > 0) = ψ0 eiU (x)t/~i= ψ0 eU0 t cos(2kL x)
2~

(5.2)

where we have ignored the global portion of the phase irrelevant to the dynamics. Now if we takea look at Appendix B we have the identity Eq. B.9 which allows us to write∞X
U0 t i2nkL xnψ(x, t > 0) = ψ0
i Jn
e
.
2~
n=−∞

(5.3)

where Jn is the Bessel Function of the first kind [1]. The pulse transforms the original wavefunction
ψ0 containing only a single component with momentum p = 0 into a superposition of components
with p = 2n~kL (n∃Z) each with probability2|ψn | = |Jn



U0 t
2~



|2 .

(5.4)

In the trapped case, the argument is identical differing only in the minor detail that ψ0 has a finite
momentum spread because of the finite spatial extent. The evolution of a trapped BEC subjected
to a suitably short pulse and immediately released from the trap is characterized by a symmetricpattern of clouds exploding away from the central cloud with relative populations given by Eq. 5.4.
We shall refer to the argument of the Bessel function Jn as the pulse parameter.
This effect, which we refer to as Bragg diffraction, provides a convenient method for characterizing the lattice depth U0 when the pulse is sufficiently short. If the amplitude evolves substantiallyduring the lattice pulse, Eq. 5.3 becomes invalid. Therefore, the pulse duration must be short compared to the timescales of the atomic motion. An obvious upper bound to the pulse duration istherefore the trap period of the ODT. The trap period of the ODT is on the order of 10 ms, andindeed the atoms can not be considered stationary for long pulses on this order. However, thelattice itself provides a timescale. In an optical lattice the atoms experience a sinusoidal potential
due to their off-resonant interaction with the 1064 nm ODT beams:
U (y) = U0 cos2 (kL x).

(5.5)

106
Where U0 is proportional to the laser intensity and the lattice wavevector is kL = π/d. Expanding
2 x2 we see thatand noting that in the harmonic approximation U (y) = 12 mωho
2
ωho=

2U0 kL2
4U0 EL
=
m
~2

(5.6)

π~
.
U0 EL

(5.7)

and easily rearrange to show that [70]Tho = 2π/ωho = √The oscillation period here provides a characteristic timescale of the lattice relevant for the calibration method described above, and satisfies the physicist’s biological requirement to express physicalsystems in terms of harmonic oscillators whenever possible. The pulse duration must be only afraction of Tho for the atoms to be considered stationary. In the lab the validity of Eq. 5.4 isevaluated using the linearity of the pulse parameter with respect to the pulse duration as describedin the next section.
5.3.2

Deep Lattice calibration

The results of the previous section form the basis for absolute calibrations of the latticepotential. A priori consideration of the beam geometries provides an estimate of the lattice depth,
but imperfect alignment of the beams with respect to each other and the atoms, uncertainties in the
beam size and shape at the atoms, and non-ideal polarization make an entirely a priori evaluationof the lattice potential unwise. The diffraction of non-zero momentum sidebands due to a shortpulse is immune to these uncertainties since the evolution according to Eq. 5.4 provides a measureof exactly how strong any potential of the form of Eq. 5.1 is at the location of the cloud. In otherwords, Bragg diffraction only measures the amplitude of the Fourier component of the potential
corresponding to that of the lattice, and does so accurately and precisely where the atoms are, nocalculation of beam powers required.
To characterize the lattices a small BEC is produced, then exposed to a short lattice pulse
at the same moment the main ODT beams (H2 and V) are turned off. The cloud, or perhaps

107

Figure 5.2: Time-of-flight absorption images of a condensate taken after variable-duration latticepulses for fixed lattice power. Pulse duration for top cloud is zero (no lattice pulse) and increasing
by ten microseconds each row. Bottom row pulse duration is 50 microseconds. Each row is thedata from a single experimental cycle.

clouds, are imaged from the side after >20 ms of time-of-flight. The entire density distribution (ofmultiple clouds corresponding to different multiples of the lattice recoil momentum) is then fit witha single global condensate size and shape given by a Thomas-Fermi profile where the amplitude ofeach mode is given by Eq. 5.4. The separation between adjacent modes is simply ∆x = 2n~ktT OF
assuming purely ballistic motion. Several images of clouds for different values of the lattice pulse
parameterα(U0 , t) =

U0 t
2~

(5.8)

are shown in Fig. 5.2, where the pulse parameter was measured for various times at a fixed lattice

108

(a)

(b)

Figure 5.3: Lattice calibration using Bragg diffraction. In (a) the entire multi-component cloud is fit
according to a single Thomas-Fermi profile where the amplitude of the several clouds is constrained
according to Eq. 5.4 with a single pulse parameter. In (b) the populations of the first nine diffracted
order from lattice pulses versus pulse duration. Each curve is an independent fit to the given mode
from which the pulse parameter may be extracted. The pulse parameters for each order agree tobetter than 1%. We use the global pulse parameter from (a) to characterize the lattices, (b) is
included for reference.

power. When the pulse duration becomes too long for Eq. 5.4 to be valid, the density distribution
will no longer resemble the precise prediction of Eq. 5.4. This can be seen in Fig. 5.3a where thepulse parameter saturates beyond a certain pulse duration, indicating where the measurement is
not providing an accurate measurement of the lattice depth. By checking that we remain in thelinear regime of Fig. 5.3a we ensure accurate characterization of the lattice depth.
5.3.3

Shallow lattice calibration

The purpose of the lattice is to write a well understood density fluctuation on an in-situ
cloud as an analogy to intrinsic fluctuations. Therefore, in-situ images of perturbed clouds mustbe translated into absolute units. Diffraction experiments in situ are certainly possible, but are illsuited for characterizing the lattice since presence of the trapping potential spoils the free expansionof the various momentum modes. Furthermore, the Bessel function distribution of Eq. 5.3 is difficult
to fit until momentum populations are well separated, but the trap means they may never separate

109

(a)

(b)

Figure 5.4: Averaged, squared Fourier transforms of the density with and without the X-lattice.
The color scale is the same for both images, and the natural log of the squared Fourier transformis shown. Prior to computing the Fourier transform the average density profile in the absence ofthe lattice is subtracted from each individual cloud. On the left fifteen clouds without lattices wereimaged, and no structure is apparent at kx = 1.73µm−1 . In contrast, the presence of the lattice isapparent in (b). Only four density profiles are averaged in (b), making it slightly noisier.

at all, and even if they do separate, the field of view of the in-situ imaging system is far too small dueto the large magnification. The only signature of the lattice available to us for high magnification,
in-situ imaging is the density distribution itself; it is in the Fourier transformation of the data
where the seed lattice perturbations are evident, as in Fig. 5.4. Extracting the lattice depth fromthe density profiles themselves turns out to be complicated by imperfect imaging as well, making a
careful calibration chain between in-situ images and time-of-flight Bragg pulse measurements of thelattice depth indispensable. The first link of the calibration chain is the absolute measurement of
the lattice depth possible for deep enough lattices, as described in the last section. A photodetector
in the beamline of the unused X-lattice light (see Figure 2.10) simultaneously records the residual
lattice power. This photodetector is directly proportional to the depth of the X-lattice and Ylattice, with the proportionality constant dependent on attenuation and alignment. It is criticalthat the alignment of the lattice optics for both the X and Y lattices be robust enough that the

110proportionality constant between the lattice PD and the lattice depth at the atoms be stable.
Substantial work on the optomechanical stability of the lattice optics was required to achievesufficient stability. Alignment and optimization of the lattices is performed by maximizing thelattice pulse parameter for fixed power while adjusting each of the four beams, two for each lattice.
When alignment is completed several neutral density (ND) filters are placed in the latticebeam path to attenuate the beam. Electronically, the lattice power may only be reliably controlled
over a factor of about 20 using the H1 AOM which is insufficient dynamic range. The positionof these filters is important as ND filters may be damaged and cause intensity dependent effects
on the beam profile, pointing, and attenuation. The lattice PD conversion factors are measured
using the most powerful lattice possible given the overall lattice attenuation, and so in practice
the attenuation is empirically chosen that maximum lattice power is just sufficient to allow Bragg
diffraction measurements to be validly performed, this is the reason that only one lattice depth
is used to calibrate the lattice PD in Fig. 5.5a; if the lattice were weaker the diffraction method
could not be used. The power must also be low enough that additional ND filters will operate in
the linear regime.
Once the system is well characterized and robust such that one may be confident in inferring
the true lattice depth at the atoms from the lattice PD, the in-situ signature of the lattice may
be calibrated. In real space, the lattice is quite subtle, and relatively deep lattices are needed todirectly see the perturbation in the density. In Fourier space the lattice is quite visible, as shownin Fig. 5.4, where on the left a BEC with no lattice applied exhibited no obvious feature at the
X-lattice wavevector of kx = 2π/(3.63µm), and on the right a very clear peak is prominent. The
region around the peak is integrated and plotted versus the known power for various lattice depths,
as in Fig. 5.5b. The inegrated quantity is then fit to the form
F (U0 ) = F0 + AU0p

(5.9)

where the exponent p should be approximately two since we are squaring the FFT of the density.
Due to noise and imperfect imaging p is a floating parameter of the fit with the result p = 1.90(11).

111

(a)

(b)

Figure 5.5: Calibration of the in-situ X-lattice. In (a) the lattice depth is determined in time-offlight imaging following a short pulse and compared to a photodetector in the X-lattice beamline(see Fig. 2.10). The slope of the linear fit in (a) is the conversion factor which determines thelattice depth in (b) based on the measured PD voltage. For (b) the clouds are imaged in situ forvarious known lattice depths and the square of the Fourier transform of the density is computed
and subsequently integrated in the region of Fourier space corresponding to the X-lattice.

If p = 2 is fixed, the intermediate lattice depths of 10-40 nK are systematically biased. Several
images with no lattice at all are taken and used to fix the background signal F0 so the fit of theintegrated lattice peak only has the amplitude A and the exponent p as floating parameters. This
provides the calibration of lattice depth vs in-situ measurement we require.

5.4

Lattice implosion experiments
In this section we discuss implosion experiments in the presence of a seed potential. For theseexperiments the cloud is cooled to sufficiently low temperatures that no thermal component was visible in time-of-flight measurements. The atoms undergo all-optical evaporation in the |1, −1i state.
The density is determined by what value with which the scattering length the cloud equilibrates;
larger values of a produce a larger, less-dense cloud. The lattice potential is then adiabaticallyapplied over 200 ms and held for the same amount of time to ensure the system is in equilibrium
at the time of the quench. The quench is performed by using a resonant microwave π-pulse to112transfer the cloud from the |1, −1i state to the |2, −2i state; the duration of the microwave pulse is∼7 microseconds. Following the quench the cloud evolves for a variable length of time we shall callthe implosion time before the cloud is imaged in situ using the |2, −2i → |3, −3i cycling transition.
5.4.1

Lattice sanity checks

We performed some experiments with the lattice as sanity checks to ensure we had an accurate
understanding of how the density perturbation behaves after the lattice potential is removed. The
first such experiment was performed by turning the lattice potential off a variable amount of time
before the imaging and observing the relaxation and revival of the lattice signal, as show in Fig.
5.6a. The quench is always performed since |2, −2i is both the imaging state and the state withattractive interactions, but the time between quench and imaging is 5 µs, negligible in comparison
to the dynamics. Only the first few oscillations of the density perturbation are observable. From
the timescale of the oscillation evident in Fig. 5.6a we performed implosion experiments where
the lattice was turned off at a moment where the trajectory of the density oscillation was either
positive, negative, or flat, see Fig. 5.6b. If the lattice remains on, the initial trajectory is flat but
the perturbation begins to grow rapidly. If the lattice is turned off at the moment of quench thetendency for the perturbation to relax is in conflict with the instability of the mode.
5.4.2

Persistent lattice

The first implosion experiments with our control lattices were performed without removing
the lattice seed-potential following the quench. A comparison between numerical simulation andexperiment is given in Fig. 5.7 for an initial lattice depth of ∼8 nK [62]. The numerical results
have been adjusted with an overall scaling factor which accounts for imperfect contrast but which
does not affect the time dependence. These data were taken with a cloud of 4.44e5(2.3%) atoms
and trap frequencies [ωx , ωy , ωz ] = 2π × [74.8, 74.7, 72.7] giving a peak density 10.3e13 cm−3 andchemical potential 95 nK. The agreement between the simulation and the prediction is excellent,and even reflects the drop off in the data at long times. The drop-off observed in Fig. 5.7 is caused

113

(a)

(b)

Figure 5.6: Perturbation response to turning off the lattice. In (a) the lattice is turned off avariable time before imaging and the perturbation relaxes before reviving with opposite phase.
If the system were an infinite uniform (except the perturbation) superfluid we would expect theoscillation of the perturbation amplitude to last indefinitely. In (b) we see the influence of thistransient response depicted in (a). When the lattice is turned off 100 µs before the quench (green)
we see the perturbation continuing to shrink. The initial trajectory of the amplitude is downward.
When the lattice is turned off 300 µs before the quench (blue), the perturbation has just begun
to grow again, and we see the perturbation grow from a smaller initial amplitude and continue togrow. The initial trajectory here is upward. If the lattice is left on during implosions, we see thebehavior depicted by the red squares. The initial trajectory is flat in this case. If the lattice isturned off at the moment of quench, the initial trajectory remains flat (black circle, gray stars) but
the lattice does not grow very much.

by both number loss through three-body recombination and the increasing deviation of the density
perturbation from a perfect sinusoid.
We fit the initial growth also to the functionñkx (t) = ñkx (0) cosht
τ

(5.10)

where ñkx (0) is the initial size of the perturbation and τ is an exponential growth constant. The
growth time constant τ does depend on the time range to be included in the fit, so to capture theearly times a cutoff time for the fit domain is chosen. We select 300 µs to be our cutoff time andfind that the perturbation growth time constant is 176(7) microseconds. Fitting the simulation
results to Eq. 5.10 over the same domain as the experimental data gives consistent results, which

114

Figure 5.7: Growth of the X-lattice seed perturbation at kx = 1.73 for the 3.63µm X-lattice. The
solid line indicates numerical simulations using the GPE. Open circles correspond to individual
experimental cycles. The simulation is fit to the experimental data with a single overall scaling
factor.

is unsurprising given the good agreement apparent in Fig. 5.7.
In the introduction to this thesis the amplification of fluctuations at a given wavevector was
motivated through the dispersion relation"
 2 2 2 #1/2
gn 2 2~ k
(k) =
~ k +
,
m
2m

(5.11)

which becomes imaginary for negative scattering lengths. We argued that the timescale for the
implosions would grow with a time constant τ = ~/(k), how has this prediction been borne
out? In Fig. 5.8a the dispersion relation for the peak and average density of the cloud is given
and in Fig. 5.8b are the expected growth timescales ~/. The timescales predicted by τ = ~/fare significantly slower than the observed timescales and the simulation results. The observed
timescales were independent of the initial lattice depth as indicated by Fig. 5.9, as expected.
5.4.3

Vanishing lattice

To better reflect implosions seeded by intrinsic fluctuations of the BEC, the seed potential
should not persist during the implosion process. After all, the fluctuations we wish to magnify are

115

(a)

(b)

Figure 5.8: Dispersion relation (k) and corresponding timescales τ = ~/(k) expected for the
peak and average density of the condensate. The wavevector of the X-lattice is 1.73 µm−1 and isindicated with a vertical black line. The timescale predicted by the dispersion relation alone aresubstantially slower than indicated by Fig. 5.7. Note the change in the y-axis labels of (a) aboveand below the  = 0 reference line, where below the line the energy is imaginary. The blue (orange)
densities indicated are the initial peak (average) densities of the data of Fig. 5.7.

Figure 5.9: Growth of the X-lattice seed perturbation at kx = 1.73 for the 3.63µm X-lattice vsinitial lattice depth. For deeper lattices, the local density may be sufficiently high as to change thegrowth rate, but no dependence is observed for shallow lattices.

DatasetName

λL (nK)kL(µm−1 )

YHD
YLD
XHD
XLD

2.45
2.45
3.63
3.63

2.56
2.56
1.73
1.73

Seed
Depth
(nK)
11.3
11.3
6.2
11.3

Peak
Density(1013 cm−3 )
10.3
5.42
8.98
4.81

Average
Density(1013 cm−3 )
5.89
3.10
5.13
2.75

ChemicalPotential
(nK)
95.3
118.3
83.1
105.1

116

Table 5.1: Conditions for the four experiments performed with the lattice potential removed at the
moment of the quench.

not reinforced by any external potential. Therefore, we performed several implosion experimentswith a seed potential which vanishes at the time of the quench. Implosions seeded by either the Xlattice or Y-lattice alone were performed for high density and low density clouds. The temperatureof the cloud was low enough that a precise measurement could not be made as the condensate
did not show a thermal component which could be fit, suggesting the temperature to be below
80 nK, as temperatures above this contain a fitable thermal component. The conditions for theseexperiments are given in Table 5.1.
The evolution of the ky = 2.5 µm−1 perturbations of the Y-lattice is shown in Fig. 5.10
for both high and low densities. Results for the X-lattice with kx = 1.73 µm−1 are shown in Fig.
5.11. Numerical simulations for each scenario are plotted alongside the experimental data. The
Y-lattice data do not agree with simulations well. The reasons for this are complex and perhaps not
fully understood. One possibility is a systematic error in determining the density. This possibility
is explored for The Y-lattice implosions in Fig. 5.12 by comparing simulations performed for arange of densities. The high density implosion data can be reasonably well accounted for by a 13%
underestimation of the density, but the low density data cannot be explained in this way.
A separate contributing factor is the noise floor of the lattice depth measurement. The Ylattice wavenumber is near the edge of our resolution, and the seed depth of 11.3 nK was selected to
be deep enough to be detectable, but shallow enough to remain a small perturbation on the chemical
potential, which is on the order of 100 nK. Therefore, the relaxation predicted by simulations will
not be visible as the density perturbation shrinks below the noise floor. This is certain to impact

117

(a)

(b)

Figure 5.10: Perturbation depth vs implosion time for clouds seeded with a 11.3 nK perturbationfrom the Y-lattice where ky = 2.5 µm−1 . The initial density-weighted density is 5.89 × 1013 cm−3 in(a) and 3.10 × 1013 cm−3 in (b). Solid lines are numerical simulations of the first 1000 microseconds
after the quench [62]. In the high density case (a) the ky mode is unstable only near the center ofthe cloud, and we see initial decrease of the perturbation, before it begins to grow rapidly at latertimes. In the low density case (b) the ky mode is stable throughout the cloud

118

(a)

(b)

Figure 5.11: Perturbation depth vs implosion time for clouds seeded with a perturbation from theX-lattice where kx = 1.73 µm−1 . Solid lines are numerical simulations of the first 1000 microseconds
after the quench [62]. The initial density-weighted density is 5.13 × 1013 cm−3 in (a) and 2.75 × 1013
cm−3 in (b). In the high density case (a) the seed perturbation is 6.2 nK and the kx mode isunstable throughout most of the cloud and immediate growth follows. In the low density case (b)the seed perturbation is 11 nK and kx mode is unstable only near the center of the cloud, andthe pertuirbation was observed to grow more slowly than in the higher density cloud. Numericalsimulations indicate an initial decrease of the perturbations amplitude in (b) prior to growth atlater times; this is not evident in the experimental data

119

(a)

(b)

Figure 5.12: Perturbation depth vs implosion time for clouds seeded with a 11.3 nK perturbationfrom the Y-lattice where ky = 2.5 µm−1 . To illustrate the possibility of a density calibration error
numerical simulations for a range of densities is plotted with the initial depth fixed to 11.3 nK (solid
curves) [62]. The solid curves have no free parameters. In (a), we see that if the actual density is
6.66e13 cm−3 the simulation better reflects the observed growth. This would require a systematic
error of 13%, a realistic possibility. However, the experiments with lower density shown in (b)cannot be explained through density miscalibrations alone. One factor guaranteed to be present is
the noise floor near 10 nK; this prevents the initial decrease of the perturbation amplitude frombeing clearly observed. The legend gives the density-weighted density in cm−3 for the simulation
results.

120

(a)

(b)

Figure 5.13: Fourier space images of in-situ data. In (a) the peak it sue to the Y-lattice, in (b) it isdue to the phantom lattice. A magenta circle indicates the rgion used to integrate the peaks. The
Y-lattice is very near the phantom in Fourier space, and since the phantom cannot be removed, it
effectively behaves similar to the persistent X-lattice of Fig. 5.7. This can contaminate the Y-latticesignal.

the low density Y-lattice data of Fig. 5.10b and 5.12b, since for any plausible density error, theperturbation is predicted to relax in the first few hundred microseconds.
A third contributing factor to the tension between numerics and observations exemplified inFig. 5.10 is due to what we call the phantom lattice, see Fig. 5.13b. This is a lattice presumablycaused by reflections of the H2 ODT beam which interfere to create a subtle lattice potential whichcannot be removed without turning off the H2 beam, which is required to hold the atoms against
gravity. The phantom lattice is very near the Y-lattice in Fourier space, as shown in Fig. 5.13. If
the contribution from the phantom is not sufficiently masked out in Fourier space, effects like theearlier than expected growth of the Y-lattice seen in Fig. 5.10a may occur.
.

121

Figure 5.14: Influence of the phantom lattice on the integrated Y-lattice signal. Each curve corresponds to a different size mask (in Fourier space) applied to mask the phantom lattice prioir
to integrating the Y-lattice Fourier peak. Unlike Figs. 5.10b and 5.12b, the lattice signal is notconverted to perturbation amplitude. It is clear that as the mask on the phantom lattice peak in
Fourier space decreases, more of the phantom lattice signal leaks into the Y-lattice data. For Figs.
5.10b and 5.12b, a mask of 0.4 µm−1 was used.

5.5

Contrast correction
As mentioned at the close of Chapter 3, one of the major challenges faced in this experimentis limited knowledge of the transfer function of our in-situ imaging. In particular, we are concernedthat the fidelity with which we image the density fluctuations is worse than expected. Even for
extremely deep lattices for which the cloud is divided into several non-overlapping local condensates,
we were unable to directly see the density drop to zero between lattice sites. The impact of thiscontrast suppression on the data analysis of unseeded implosions will be discussed in Section 6.4.
To attempt to compensate for this problem we returned to our lattice experiments to evaluate

122
for kx = 1.73 and ky = 2.5 our efficiency in detecting density structure at these two values ofk. To this end, we created synthetic density distributions numerically by defining a 3D trap
potential V0 (r) with similar trap frequencies to our experiments, defined a value of the chemical
potential, defined a lattice potential Vx (r), and filled in the density according to the Thomas-Fermi
prescription n(r) = Θ(µ − V0 (r) − Vx (r)) where Θ is the Heaviside function. These 3D density
profiles were summed along the z-direction to mimic the profiles observed in the experiment. The
square modulus of the two-dimensional FFT was calculated for each density profile. In Fig. 5.15
the density profiles and cross section in real space and frequency space are shown for Vx = 0 while
in Fig. 5.16 the lattice depth is 20%. The lattice potentials Vx (r) are strictly nonnegative, so thedensity was uniformly renormalized to ensure constant N between clouds of differing lattice depths.

To quantify the strength of the lattice for different fractional lattice depths, the integrated
volume of the region indicated in magenta in Fig. 5.16 is divided by the integrated volume of thecentral peak indicated by the green circle of the same figure. In Fig. 5.17 this ratio is shown forvarious lattice depths as a fraction of µ. Let us call Vx /µ = V˜x and the ratios between the lattice
peak and the bulk peak RN (V˜X ) where the subscript N indicates it is computed numerically. We
then reanalyzed our in-situ lattice calibration curves and calculated the corresponding ratio weshall dub RE (V˜x ) where the E indicates experimental data. The experimental data were fit to a
curve of the forms
RE (V˜x ) =

c20 +



Vx 2
A
µ

(5.12)

where the constant c0 is fixed by reference shots, see Fig. 5.18. We then compute the ratio
RN (V˜X )/RE (V˜x )|c0 =0 to determine by what factor the density fluctuation is suppressed in our
measurement. The quantity RE (V˜x )|c0 =0 indicates we are only comparing the lattice contributionof our calibration curve, i.e. we set c0 = 0 but do not modify A from its fitted value. We see from
Fig. 5.19 that the suppression factor we infer from this procedure depends on what normalized
lattice depth V˜x is chosen as the reference. For the data in Fig. 5.18 the unperturbed chemical

123

(a)

(b)

(c)

(d)

Figure 5.15: Synthetic Thomas-Fermi profiles without any lattice potential in (a) and its crosssection in (c). the square modulus of the FFT is shown in (b), and its cross section in (d). thecolor scale for (b) is chosen to match Fig. 5.16

124

(a)

(b)

(c)

(d)

Figure 5.16: Synthetic Thomas-Fermi profiles with a lattice potential equal to 20% of the chemical
potential in (a) and its cross section in (c). the square modulus of the FFT is shown in (b), andits cross section in (d).

125

Figure 5.17: Ratio between the volume of the lattice peak and central peak in Fourier space.

potential was 82 nK, and we can see that at a lattice depth of 16 nK is sufficiently deep to havea clearly resolved lattice signal compared to the unseeded cloud. We therefore select 20% as ourreference depth to determine the suppression factor from Fig. 5.19, with the conclusion that at√
kx = 1.73 our sensitivity is suppressed by a factor of 27 compared to a perfect imaging system.
The data of Fig. 5.19 are derived from the squared magnitude of the FFT, so the suppression factor
in terms of the density is the square root of the quantity plotted.
Repeating this procedure for the other three data sets of Table 5.1 yielded interesting results.
√
The suppression factor for the low-density X-lattice data set was found to be 43 , over 50% greaterthan the high-density X lattice case. The Y-lattice high-density suppression factor was found to√√
be 155 while the low-density suppression factor was 135 . The Y-lattice wavenumber lies atthe edge of our resolution, so the suppression factor was expected to be much larger than for the
X-lattice wavenumber. The density dependence at kx = 1.73 µm−1 was not expected, but may be
the result of a slight misalignment of the X-lattice along the line of sight and the fact that the lowdensity cloud is larger than the high-density cloud. Imaging through a larger cloud with a slightly
misaligned density corrugation will result in lower contrast. That this density dependence was much

126

Figure 5.18: Calibration data for the high density implosion experiments with a X-lattice seedpotential discussed in the previous chapter.

Figure 5.19: Ratio of ideal sensitivity to observed sensitivity for the X-lattice Fourier component.
Depending on which normalized lattice depth fraction V˜x is chosen as the reference, our inferredvalue of the suppression factor for our observations varies from 23 to 37. We select a 20% lattice
depth (V˜x = 0.2) to balance between sufficient signal to noise with the desire to investigate thesmall perturbation limit (see Fig. 5.18).

127
weaker for the Y-lattice may be because the Y-lattice slices are better aligned to the imaging axis.
However, the fact that the Y-lattice suppression factor is actually lower for the low density caseseems inconsistent with this hypothesis, although the closer agreement in the Y-lattice suppression
factors may indicate that they are essentially the same, and the difference in suppression factors is
simply a statistical fluctuation.

Chapter 6
Intrinsic Fluctuations

This chapter concerns our use of implosions to amplify and measure intrinsic fluctuations inbulk 3D BECs. We sought to test the theoretical predictions of the structure factor S(k, t) andits evolution following the quench to negative scattering length. We were unable to clearly observe
many of the effects predicted. Based on the forthcoming discussions and those of the previous
chapters, particularly in regard to the challenges in imaging a bulk 3D BEC, we shall approach ouranalysis assuming as correct the time dependence for the growth of the structure factor predicted by
the theoretical model introduced at the end of Chapter 4 in order to test the model’s prediction forthe pre-quench structure factor. We shall begin with a discussion of the challenges we faced in ourdata collection and analysis and potential sources of error. We shall then proceed to a descriptionof the data collection process followed by a review of the results for four data sets covering both
high and low regimes of temperature and density. We shall focus on the wavenumbers for which
we have calibration data.

6.1

Discussion of sources of error
We begin our discussion by highlighting the difficulties we encountered in the experiment.

Before we begin, it is useful to be explicit about what sources of noise exist in the system. First andforemost, we are seeking to measure the density fluctuations of BEC, and density fluctuations may
arise in several ways. An ideal BEC in an perfect harmonic potential, completely isolated from the129
surrounding environment, including a measurement apparatus,1 will exhibit density noise even atzero temperature due to the discrete nature of atoms and the associated shot-noise. An interacting
BEC will exhibit interaction effects which deplete the single particle ground state and modify thedensity fluctuation spectrum. The effects of the interactions is not necessarily to increase thedensity fluctuations for all wavelengths, and indeed the repulsive interactions used in this thesis
have been shown to suppress the structure factor at small k [72]. Figure 6.1 shows the predictedsuppression for our densities at zero temperature. Let us call the total density noise which exists at
T = 0 quantum density noise. The zero temperature BEC corresponds to the quasiparticle vacuum
state in the Bogoliobov picture, as discussed in Chapter 4. At any finite temperature T > 0 someof these quasiparticle states will be excited and thermal density fluctuations will exist. Let us callthese density fluctuations thermal density noise. The contributions of both quantum density noise
and thermal density noise are included in the expression for the structure factor which, for a BEC
in equilibrium at positive scattering length, is given by
S0 (k) =1 ~2 k 2coth(k) 2m



(k))
2kB T


.

(6.1)

The cloud may also have density noise which we shall call ‘superthermal’. If one taps their knuckle
on a table, the table will vibrate with some frequency spectrum, but these excitations are notthermally excited. If one were to compute the temperature required to thermally populate thevibrational states the table occupies following the knuckle tap the resulting temperature would bemany orders of magnitude greater than the actual temperature. Eventually, dissipative processes
like turbulence and friction will distribute the energy as heat and the noise spectrum would become
thermal (to the extent that one may neglect the quantum noise of classical object like a table).
Superthermal density noise in our condensate may arise in several ways, such as sloshing or stirring.
Uncharacterized contributions to the confining potential may also cause superthermal density noise,such as the phantom lattice mentioned in Chapter 5. We shall include quantum density noise,
thermal density noise, and superthermal density noise in the broader group of real density noise,1

We are ignoring the measurement problem; possibly the concept of a physical system in the absence of a measuringapparatus is incoherent, but we aren’t making metaphysical claims

130

(a)

(b)

Figure 6.1: Structure factor S0th vs scattering length for our data sets. The longer wavelength
fluctuations are more strongly suppressed than those at shorter wavelengths. The densities used
are hni = 5.5e13 cm−3 for the high-density data and hni = 3.2e13 cm−3 for the low-density data.
The fact that S0th for various densities and wavenumbers all converge to 1.0 in the limit of vanishinginteractions clarifies the normalization of S(k). S(k) = 1.0 is the result for naïve, poisonnian shotnoise in the density.

because they reflect the actual density configuration of the condensate.
Additionally we have the effects of imperfect resolution and knowledge of the imaging point
spread function. These include the contrast fidelity discussed at the end of the last chapter andmay stem from aberrations in the imaging system. We shall refer to these as imaging calibration
errors.
Finally we will have noise contributions from the measurement process. These include photon
shot noise on the detector, imaging artifacts like interference fringes, and anything which does not
represent actual density structure of the atomic cloud. Any contribution to the measured density
fluctuations that does not arise from the real density noise or imaging calibration errors we shall callimaging noise. In other words, imaging noise refers to sources of inferred density structure which are
not due to the actual density distribution of the atoms or are unrelated to the imaging point-spreadfunction. Imaging noise itself may be divided into two categories, cloud-related imaging noise andcloud-unrelated imaging noise. Cloud-related imaging noise contributions are those which only
appear in the presence of the atomic cloud, such as certain types of fringes and probe saturation
related noise. In the following discussion we shall consider in turn imaging calibration error, thenimaging noise, and finally real density noise. First, we shall elucidate the strategy we take in131analyzing the data in the presence of these sources of noise.
6.1.0.1

Strategy for coping with noise

In Chapter 4 we discussed how the static structure factor is predicted to evolve following aquench to negative a. The structure factor evolves according to




~2 k 2 n(gi + |gf |)1 ~2 k 2(k)2 f (k)tS(k, t) =
coth
1+
sin
(k) 2m
2kB T
m|f (k)|2
~

(6.2)

S(k, t) = S0 (k)(1 + Srel (k, t))If we had perfect quantitative understanding of both our imaging transfer function (resolution
fidelity as a function of k) and the imaging noise contributions to our images (such as photon shotnoise, or unsubtracted probe beam structure), and we were confident that superthermal density
noise is a negligible contribution to the real density noise of the condensate, we would be able
to robustly test the predictions for both the structure factor S0 of the pre-quench cloud and thefractional growth Srel (k, t). Unfortunately, our understanding of our resolution is limited as is thatof the imaging noise introduced by the imaging procedure. We have quantitative understanding of
our resolution only at the two lattice wavenumbers, and have very little quantitative understandingof the imaging noise. Therefore, we adopt the strategy of assuming the accuracy of the relative
growth factor Srel (k, t) and use this to to extract a value for the noise background and the initialstructure factor S0 (k). This value of S0 (k) may then be compared with theory. We therefore fit
the data to
Sraw (k, t) = c0 (k) + Sexp (k, t = 0)(1 + Srel (k, t))(6.3)

where c0 (k) and S0exp (k) are the fitting parameters. Time-independent imaging noise is characterized by c0 (k) and S0exp (k) is the experimentally determined initial structure factor. The ‘raw’
signal Sraw (k, t) is2 the squared Fourier amplitude normalized by the atom number, as discussed in
Section 4.6; the details of computing Sraw (k, t) from our in-situ images will be discussed in Section
6.3. In terms of real density noise, imaging noise, and the resolution correction discussed in Section2

Since the raw signal Sraw (k, t) is comprised of both real density noise contributions and imaging noise, it is nottechnically correct to call Sraw (k, t) the ‘raw structure factor,’ although this is an easy detail to forget.

132
5.5, the raw signal is
Sraw (k, t) ∝RealDensityN oise(k, t)+ ImagingN oiseCloud (k) + ImagingN oiseN oCloud (k). (6.4)
Resolution(k)6.1.0.2

Imaging calibration error

The top imaging system was designed and intended to yield a diffraction limited resolution of
1.1 µm, but we were unable to achieve this performance and estimate the actual resolution to be 2.5µm. Additionally, the fidelity of our imaging structure well within the resolution limit seemed to becompromised. We calculated the contrast suppression factor for our two lattice wavevectors in theprevious chapter. In Section 6.4 we will discuss in more detail how we implement this correctionin our analysis and its influence on our results.
6.1.0.3

Imaging noise

The imaging noise is dominated by a few effects which may influence the inferred density
profile for reasons unrelated to the calibration of the optical transfer function of our in-situ imaging,which we addressed in the last section, or due to the real density noise to be discussed in the next
section. Sources of imaging noise include fringes and photon shot noise on the detector. In Chapter
3 we discussed fringe removal from absorption images. This procedure works well when the fringes
themselves do not depend on the atomic cloud. A fringe may appear due to a speck of dust drifting
through the side of the probe beam as it propagates from the atoms to the camera, and typicallythis will not be problematic. However, in some situations the fringes depend strongly on the atomic
cloud. In my early years at JILA, while I was working with rotating gases and Bose polarons, we
had a particularly pathological fringe which plagued our imaging for several days, but only appeared
when the atoms were present. The culprit turned out to be a speck of dust which sat down-beamfrom the atoms and lay within the shadow of a large atomic cloud such that when atoms were there,
the speck did not scatter any light as it was in the shadow; but for the light frame, taken with noatoms, the speck scattered a great deal of light and produced rings on the light frame which wereabsent in the shadow frame.

133
It was only when our laser became unlocked and several images were accidentally takenwith no atoms did we notice the fringe didn’t appear in the OD analysis, and the rings fromthe scattered light suddenly appeared in the “shadow” image as well, allowing the fringe to be
perfectly subtracted away in the computation of the OD, but only when the atoms were absent.
That frustrating experience taught me the valuable lesson that imaging noise can manifest in
unexpected ways, and the value that absorption images with no atoms present may provide theexperimenter in understanding noise. In the limit of small pixel-by-pixel differences between thelight and shadow frame the effect on the optical depth is linear as structure subtracts nicely, butwhen the shadow frame deviates substantially (fractionally speaking) the effect of noise becomesincreasingly nonlinear.
We analyzed the Fourier spectrum of a set of images taken in the absence of atoms during
the experimental runs for the main data. Occasionally during the data scans a laser would become
unlocked or the cart would become stuck and the experiment would proceed for hours without
producing any atoms, but with each OD image taken faithfully by the machine. One such set ofover 600 images were taken during collection of the lowest temperature data set. We compute thesquare of the FFT for the same region of real-space as we do for images with atoms. In Fig. 6.2 weanalyze these images both with and without subtracting the mean OD profile for all 600 images, andthe fractional impact on the amplitude (squared) in Fourier space is shown in Fig. 6.2c. Structurefrom the probe beam is apparent, but this structure is most readily observed in the presence ofatoms. In Fig. 6.3 the 51 images taken for the tImp = 0 time step of the same data set are analyzed
in the same manner as Fig. 6.2. Structure from the probe beam is written onto the atoms when
computing the OD through the saturation intensity correction as described in Section 3.1.3, see
Eq. 3.23. This structure is subtracted out when the mean density profile is subtracted from eachindividual image, as shown by Fig. 6.3b which does not show any apparent structure for k > 1.
It remains possible that in subtracted the mean OD profile in computing the structure factor
we fail to fully eliminate the noise written onto the atoms by the saturation intensity correction. Indeed, if we attempt to compensate for shot-to-shot changes in the cloud position by post-processing

134

(a)

(b)

(c)

Figure 6.2: Squared FFT in the absence of atoms, x and y axes are wavenumbers k = 2π/λ ininverse microns. The color axis in (a,b) is ln(S0raw ). The same region of real space is analyzed foreach panel. In (a,b) no atoms were present and the logarithm is display for increased visibility. In(a) the OD image for each individual experimental cycle was analyzed with no subtraction of theaverage. In (b) we do perform the subtraction of the mean OD prior to computing the FFT. In (c)the effect of subtracting the average is shown for the central region of Fourier space. To compute(c), the array values (without taking the log) in (b) are subtracted from (a) (also without taking
the log) and divided (pixel-wise) by the values shown in (a). Therefore, (c) indicates the fractional
increase in noise for each pixel in Fourier space if the average is not subtracted; despite the visual
similarity of (a) and (b), subtracting the mean OD suppresses noise on the order of 25% for k > 1.
In (a,b) the entire Fourier profile is shown; features at k > 3 are imaging artifacts well beyond ourresolution.

the density profiles to align them to a common center, these shifts, typically only a few pixels, reintroduce a large amount of noise to the structure factor calculation.
For each data set in Table 6.1 we evaluated a number of images taken during data collection
for which there were no atoms in the same manner as to be described in detail in Section 6.3
without normalizing by the observed atom number. This provides a quantitative estimate of how
much noise the images exhibit from noise contributions present even when atoms are not present.
In Figure 6.4 the unnormalized structure factor evaluated for images both with an without atoms
present is shown for one of the high-density data sets in Fig. 6.4a and for one of the low-density
data sets in Fig. 6.4b. The data indicate that from k = 1 µm−1 to k = 2.3 µm−1 cloud-unrelated
imaging noise accounts for nearly all of the signal in the low density data, while it accounts forabout half of the noise in the high-density data. This trend was consistent across all high density
and low density data. Near the edge of the resolution limit kres = 2.5 µm−1 cloud-unrelated noise

135

(a)

(b)

(c)

Figure 6.3: Squared FFT in the presence of atoms, x and y axes are wavenumbers k = 2π/λ ininverse microns. The color axis in (a,b) is ln(S0raw ). The same region of real space is analyzed
for each panel. In (a) the OD image for each individual experimental cycle (51 total images) was
analyzed with no subtraction of the average. In (b) we do perform the subtraction of the mean
OD prior to computing the FFT. In (c) the effect of subtracting the average is shown for thecentral region of Fourier space. To compute (c), the array values (without taking the log) in (b)are subtracted from (a) (also without taking the log) and divided (pixel-wise) by the values shownin (b). Therefore, (c) indicates the fractional increase in noise for each pixel in Fourier space if theaverage is not subtracted. Contrary to the case where no atoms were present shown in Fig. 6.2,
subtraction of the mean density profile is required, as otherwise noise from the probe beam iswritten onto the atoms. In some regions of k-space the noise from probe beam may be more than
10x the atomic signal. Unlike in Fig. 6.2, we only show the resolved region, indicated clearly by
the plateau reaching out to nearly k = 3 µm−1 in (b).

begins to dominate for all cases.
6.1.0.4

Real density noise

Of course, some of the density noise we see is due to actual structure on the cloud. The
fluctuations we wish to amplify and measure by imploding the BEC are the quantum density noise,
which is present even a T = 0, and the thermal density noise. The sum of these two contributions
is given by Eq. 6.2. However, density structure may exist on the cloud for other reasons and these
density fluctuations are not included in Eq. 6.2. We refer to these fluctuations as superthermal
density noise. As we saw in the previous chapter, there appeared to be what we described as a’phantom’ lattice, which is a type of superthermal density noise which is well correlated shot to shot,a ‘static’ source of superthermal density noise. This lattice seemed to be caused by self interferenceof the H2 ODT beam as it remained visible when the atoms were confined by only that beam and

136

(a)

(b)

Figure 6.4: The unnormalized structure factor evaluated for images both with an without atoms
present. In (a) the low-temperature, high-density images with atoms yield nearly twice the signal
obtained in the absence of atoms. In contrast, (b) indicates the low-temperature, low-density
images with atoms yield very similar Sraw as images without atoms. In the language of Section
6.1 the orange circles represent the cloud-unrelated imaging noise, while the difference between theblue and orange markers is due to a combination of cloud-related imaging noise and real density
noise, which itself is comprised of quantum and thermal density noise and superthermal density
noise. We shall see that imaging noise dominates over real density noise.

all others were physically blocked from leaking into the science cell. We were able to mitigate thephantom lattice with a carefully constructed aperture, but other phantoms may exist of which weare unaware. That the phantom is visible is already a piece of luck, since if the phantom wave
vector possessed a large enough z component, then the density corrugation would be washed out
by the geometry of our imaging. If such a phantom exists due to back-reflection of the V ODTbeam, it would be entirely invisible to the top imaging system. Locally, real density structure may
exist from phantom lattices or stray reflections and thus perturb the density from our expectation.
These perturbations to the confining potential cannot be turned off as the calibration lattices could
be, and it was shown in the previous chapter how influential such a persistent structure may be on
the implosion process.
If there are unknown density features in the confining potential, then the impact of any
initial velocity the condensate possesses at the time of quench may be to generate waves as the137
condensate sloshes over the perturbation. We saw in the last chapter how the initial velocity of thelattice ripples had a profound effect on the subsequent evolution, as shown in Figure 5.6b. This is a
plausible explanation for the inconsistent temperature dependence we shall observe in forthcoming
sections. Colder clouds will actually be more sensitive to roughness of the optical potential, as(non-static) superthermal density ripples will take longer to damp away in colder clouds. During
preparation of the initial conditions for each data set, measurements were performed to ensure thecloud is at rest and in equilibrium at the moment of the quench, however we saw in Chapter 4
how fluctuations in the phase are far more important to the behavior of a Bose superfluid, and so
the possibility that minute phase or velocity fluctuations might be magnified and surpass the effectof the density fluctuations we seek to amplify cannot be discounted. Density fluctuations arisingin this way are a form of superthermal density noise which may not be correlated shot to shot as
opposed to the static superthermal noise sources mentioned above.
6.1.1

Shall I compare BEC to summer’s day?

At the beginning of this section (6.1) the example was used of tapping one’s knuckles on atable to describe superthermal density noise. A better analogy would be the surface of a pond on
a pleasant summer’s day. A perfectly calm pond3 will have fluctuations of the local water level dueto thermal excitations. These are not waves or ripples that one can see by eye, but rather the same
fluctuations which cause Brownian motion. They are microscopic, but they do exist. Only at zerotemperature would these fluctuations vanish. If a perfectly stationary boat were floating on theperfectly calm pond, the water level would contain a static dip where the boat displaces the water.
This would be most obviously seen if the pond froze and the boat were then removed; a depression
would exist in the otherwise flat surface of the frozen pond. But suppose the boat and the (liquid)
pond were struck by a summer’s breeze; the boat would begin to rock and the surface of the pond
would develop waves. In our experiment, the phantom lattice is like the boat, and the BEC is3

The pond is classical, we are not considering quantum effects in this analogy, don’t @ me.

138
like the pond.4

As the water flows around the boat, ripples will form but these ripples have no

relationship to the thermally excited fluctuations which are microscopic. Eventually, turbulence
will redistribute all the energy of these macroscopic, superthermal waves as heat and the pond willappear still once more, but still have the static dip in water level where the boat floats.
In our experiment, we know we have at least one boat on our pond, the phantom lattice.
If this phantom lattice (or others) is moving, density and velocity fluctuations will be generated
like waves on the pond. How can we distinguish waves of this sort from the thermal fluctuationsof the pond? Let’s consider a different question: how could we determine the temperature of thewater from a video of a pollen grain, or the molecular mass of an unknown liquid at a known
temperature?5

If one were to perform such an experiment, perhaps in a small dish under amicroscope, held by an undergraduate with too much coffee in their system, the challenge would
be to separate the thermal, Brownian motion of the pollen from the superthermal macromotionof the the sloshing liquid. If one wished to determine the molecular mass of an unknown liquid,
one could record the motion of the pollen at a range of temperatures, and analyze its motion in
frequency space. Based on the principles of atomic theory and the mass of the pollen, one could
fit6 a model of the power spectral density of the pollen containing two free parameters, the molar
mass of the liquid, and a background noise profile reflecting the trembling hands of the student.
We embrace a similar approach, where the validity of Srel (k, t) in Eq. 6.2 is assumed a priori anda fit to Eq. 6.3 is used to isolate the portion of our observed density which evolves with the timedependence predicted by Srel (t) from a noise background which does not.
Of course, if our input parameters to the theoretical predictions are faulty, or if other unaccounted for phenomena are present then the theoretical prediction of Eq. 6.2 cannot be relied
upon. Such errors may arise from the typical types of experimental error in BEC experiments, likedensity calibration errors or the presence of spin impurities. Discussion of these more mundane
sources of error will be deferred to Section 6.6.
4

and background gas collisions are when ‘rough winds do shake the darling buds of May’ into the pond
This is apparently the Comps 2 portion of the thesis.
6
It is left as an exercise to the reader to do this.
5

data set
Name
LTHD
HTHD
LTLD
HTLD

Temperature
(nK)
100
226
91
210

Peak
Average
ChemicalDensityDensity
Potential(1013 cm−3 )(1013 cm−3 ) (nK)
10.3
5.89
93
9.39
5.37
84
5.82
3.32
121
5.54
3.17
115

139

Table 6.1: Conditions for unseeded implosions. The data sets are labeled with HT (LT) indicatinghigh (low) temperature and HD (LD) indicating high (low) density.

6.2

Sample preparation
Clouds of a known density and temperature are prepared for four different combinations oftemperature and density, with the conditions for each data set given in Table 6.1. The system isthen quenched to a final scattering length of −33.4 a0 by way of a 6-8 µs microwave pulse resonantwith the |1, −1i → |2, −2i transition. We then wait up to 1000 µs, in steps of 40 µs, and image
the atoms using the |2, −2i → |3, −3i cycling transition. The imaging is done on resonance in the
I  Isat regime for the probe intensity. In all cases, a five microsecond gap between the quench
and the start of image acquisition is inserted which we shall ignore in our treatment. We shall referto the shortest implosion time, which technically takes place beginning five microseconds after thecompletion of the microwave pulse, as the tImp = 0 time step, and similarly ignore this short time
gap for all other time steps. The high temperature (HT) and low temperature (LT) data for a given
density (high density (HD) and low density (LD)) were taken during the same data run. So 26
experimental cycles were taken with an HTHD cloud, first with zero µs implosion time (tImp = 0),then tImp = 40 µs, and proceeding to the final implosion time of tImp = 1000 µs. The conditions
are then changed to prepare a LTHD cloud and 26 experimental cycles, one for each time step, areperformed. This procedure is repeated back and forth until a few dozen images for each time step
are collected for each temperature. The HTHD and LTHD data sets were taken together this way
over three days in early November 2021. The HTLD and LTLD data sets were taken together inlate November 2021. The primary way the densities of the LD data set are reduced compared to

140
the HD data sets is by adiabatically increasing the scattering length of the cloud to 250 a0 over
several dozen microseconds.

6.3

Computation of the structure factor
In Chapter 4 we introduced the structure factor which quantifies the fluctuation spectrum

as a function of temperature, density, and scattering length. Recall that the structure factor is theensemble average of the squared Fourier transform of the real-space density deviations normalized
by the total number N
1h|δn(k)|2 i.
N

(6.5)

Z

(6.6)

S(k) =The quantity δn(k) is computed via
δn(k) =

δn(r)e−ik·r dr

where δn(r) = n(r) − hn(r)i and hn(r)i is the ensemble average density. Let us go through thedetails of the calculation for a few data points of the HTLD data set by considering Fig. 6.5.
We first compute the two-dimensional density profile of each image for a given time step.
The mean density profile is then computed for the specific tImp , such as the tImp = 0 time step ofFig. 6.5a. For each of the individual images with the implosion time (Fig. 6.5b) a mask is applied
to both the mean (hn(x, y)i) and individual density profiles ni (x, y) and the total atom number forthe mean (individual) density hNW i (NiW ) within the unmasked region is computed by summing
the pixel values and multiplying by the effective pixel area. The subscript on NW stands forwindow to indicate we only care about atoms within the unmasked region of real space. The mask
is constructed to isolate only the portion of the image containing the condensate, and the samemask is used for all the images for a given tImp . The density deviation δni (x, y) for the individual
profile ni (x, y) is computed as δni (x, y) = ni (x, y) − C hn(x, y)i where C = NiW / hNW i is a scaling
factor to account for total number fluctuations. An example δni (x, y) is displayed in Fig. 6.5c. The
structure factor Si (k) = |δni (x, y)|2 /NiW is then computed by performing a two-dimensional fast
Fourier transform (FFT), squaring the magnitude, and dividing by NiW . The resulting structure

141
factor is shown in Fig. 6.5d. The structure factor is symmetric about the origin, but we display theentire plot for aesthetic reasons.
Unlike the last chapter, the condensate is not seeded with a well defined density perturbation,and so we extract the fluctuation amplitudes from an annulus in Fourier space. We divide thestructure factor Si (k) into rings of constant width in k-space as shown in Fig. 6.5d and computethe mean value for each ring. We chose a ring width of 0.2 µm−1 as this reflects the characteristic
size of the features in k-space; if the rings were much wider, we would lose resolution in our
determination of S(k) and if the rings were smaller, neighboring rings would become increasingly
correlated and signal-to-noise would suffer due to the decrease in the number of pixels in k-space. It
is evident in Fig. 6.5d that the innermost rings are influenced by the low-k peak which reflects thebulk size, shape and position of the cloud. Number fluctuations in the absence of size or position
fluctuations would manifest only at k = 0 if the scaling factor C were fixed to unity, but since
the mean density is scaled before computing the density deviation the k = 0 term is near zero byconstruction. It is not identically zero due to machine precision. Additionally the phrase “number
fluctuations in the absence of size fluctuations” does not reflect reality because the size of the cloudis a function of total number. Since we compute the density deviation from the average density
distribution, shot-to-shot changes in the position are actually the largest contributing factor for
k < 1. In principle these can be compensated for by shifting the optical densities of each image
to center them to within a pixel of the common center, but this cannot be performed when thesaturation intensity correction to the absorption imaging is included in the computation of the
optical depth or density profile. If the saturation correction is used, spatial structure from theprobe beam is written onto the density profile; when δn(x, y) = n(x, y) − C hn(x, y)i is computed
without shifting the images, this structure subtracts out very well, to the point of being negligible.
However, if the images are aligned to a common center, the pixel-by-pixel subtraction will actuallymagnify the probe beam’s structure.
As the cloud evolves following the quench, we expect to see features develop in k-space asdensity fluctuations are amplified. We present in Fig. 6.6 a tImp = 800 µs rendition of the tImp = 0

142

(a)

(b)

(c)

(d)

Figure 6.5: Example data for calculation of the structure factor for the tImp = 0. Panels (a,b,c) arein real space, with the x, and y axes given in pixels; the effective pixel size is 0.32 µm, a red scalebar is shown for reference. The color axes of (a,b,c) are column densities in m−2 . The average twodimensional density profile, displayed in (a) is calculated from 33 individual experimental cycles.
In (b) on of those 33 density profiles is displayed. A mask is applied to both density arrays andthe density deviation δn(x, y) = n(x, y) −R C hn(x, y)i is computed where C is the scaling factor
described in the main text to ensure that δn(x, y) = 0 within the unmasked region; the resulting
density deviation array, with mask applied, is shown in (c). The structure factor Si (k) is computedfrom the FFT of (c) and displayed in (d) along with the eleven annuli over which we average Si (k).
The innermost ring extends from k = 0.6 µm−1 to k = 0.8 µm−1 , while the outermost ring extends
from k = 2.6 µm−1 to k = 2.8 µm−1 . The first couple rings are influenced by the bulk structure of
the cloud as a whole, particularly from shot-to-shot changes in cloud position.

143

(a)

(b)

(c)

(d)

Figure 6.6: Example data for calculation of the structure factor for the tImp = 800 µs. Panels(a,b,c) are in real space, with the x, and y axes given in pixels; the effective pixel size is 0.32 µm,
a red scale bar is shown for reference. The color axes of (a,b,c) are column densities in m−2 . The
density distributions at 800 microseconds of implosion time do not look much different from thet = 0 case presented in Fig. 6.5. For the individual experimental cycle shown here in (b), the cloud
position is offset slightly from the mean profile in (a). In Figs. 6.5a and 6.5b the individual profile
was well aligned to the mean. The offset is visible in the density deviation (c) and contributes tothe low-k portion of the structure factor shown in (d). At larger k we see more fluctuations in(d) than was evident in Fig. 6.5d due to the growth of fluctuations over 800 microseconds. Notethe mask at [kx , ky ] = [0.3, −2.15] (mind the flipped y-axis) bounded by the upper magenta circle
which is introduced to mask the phantom lattice, which could not be perfectly eliminated. The
lower magenta circle shows the unmasked signal at the phantom wavevector; the prominent feature
within this region may be the signature of the phantom lattice growing, and so it is omitted incomputing the mean S(k) for the rings over which the phantom region overlaps.

144µs data of Fig. 6.5 which illustrates this growth. The phantom lattice discussed in the last chapter
may be responsible for the feature within the lower magenta circle. For this reason a mask is applied,
indicated by the upper magenta circle, although for computing the ring results both regions aremasked.

6.4

Contrast correction
One of the challenges we were confronted with in our experiment was the poor contrast ofdensity perturbations in our in-situ images. This was the problem considered in Section 5.5, andcan be understood by comparing the predicted growth at a given k compared to the observations.
Figure 6.7 shows the raw signal Sraw (k, t) Figure 6.7 shows how, under the presumption thatthe prediction is accurate, but our detection may be inefficient, the dominant contribution to theobservations at tImp = 0 must be noise which is not growing over time, which is to say it is not real
density noise. If the noise offset is subtracted uniformly from the blue data points in Fig. 6.7, thenthe remaining contribution grows with the same relative behavior as the prediction, but with anamplitude given by the fitting parameter S0exp in Eq. 6.3 such that S0exp /S0th = 0.04, suggestingthat our sensitivity at this wavenumber (k = 1.7) is only about 4%. We multiply the observed
structure factor by a factor of 27, the contrast suppression factor we calculated for the highdensity, X-lattice wavenumber kx = 1.73 µm−1 and repeat the fit to Eq. 6.3, with the result shownin Fig. 6.8. The growing component is well captured by the fit model and our inferred structure
factor of S0exp (1.7) = 4.37 is consistent with the prediction of S0th (1.7) = 4.12. Unfortunately,
applying the contrast suppression correction to the low-density k = 1.7 µm−1 rings as well as therings at k = 2.5 µm−1 wavenumber (high density and low density) did not yield the same agreementwith theory.

6.5

Results
Here we present the results of the implosion experiments investigating the structure factor.

A direct fit to the functional form with a scaling factor and background offset was performed. The

145

Figure 6.7: Observed structure factor vs implosion time for the k = 1.7 ring in the high temperature,high density data set. We fit the data to Eq. 6.3 which has two free parameters, an overall scaling
factor and an offset to account for sources of imaging noise. The prediction is given by the dotted
purple line and the fit by the solid maroon line. Clearly, the predicted growth is much more rapid
than was observed; after 170 us the prediction doubles from 4.12 to 8.32, while the observations
go from 2.28 to 2.45, an increase of only 7%. The blue line gives the value of the imaging noise
fit parameter, 2.12 for these data. If the background is subtracted out, then the fit indicates the
non-noise contribution grows by a factor of (2.45 − 2.12)/(2.28 − 2.12) = 2 as constructed.

experimental data, along with the analytic prediction of Eq. 6.2 and is shown in Figs 6.12, 6.13,
6.14, and 6.15. In those figures the correction has not been applied to any of the rings. In the next
section the contrast suppression correction factor discussed in section 5.5 is applied to the k = 1.7µm−1 and k = 2.5 µm−1 rings, and in section 6.5.2 we give the results without this correction.
6.5.1

Calibrated wavenumbers

We first present results using the contrast correction factor discussed in the previous section.
Only two values of k could be calibrated in this manner. We display the inferred structure factor
along with the fit parameters in Fig. 6.9; where we have included the results for the other valuesof k for which we lack knowledge of the contrast suppression factor to indicate the influence of the146

Figure 6.8: Refit to the model of Eq. 6.3 after applying the contrast correction factor calculated
for the high-density X-lattice wavenumber. While the noise is substantial, the growing component
matches the prediction of Eq. 6.2 very well, with the S0exp /S0th equal to 1.06.

contrast suppression. Only for the high-density k = 1.7 µm−1 ring do we clearly see the expected
temperature trend, and only the HTHD data are quantitatively consistent with the prediction. The
LTHD data lie below the high temperature data, but above and inconsistent with the prediction.
No temperature dependence is evident in the low-density k = 1.7 µm−1 ring, and the inferred
structure factor is much higher than expected. For both densities, the k = 2.5 µm−1 ring wasobserved to be much greater than the expected value. The temperature dependence of the LT andHT high-density sets has the correct trend. With the previous mentioned exception of the HTHD
k = 1.7 µm−1 ring, the ratios S0exp (k)/S0th (k) lie well above one as shown in Figs 6.9c and 6.9d.
The imaging noise, shown in 6.9e and 6.9f is quite large, as expected.

147

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.9: Results with the contrast correction applied according to Section 5.5. The correction
is only applied to the k = 1.7 µm−1 and k = 2.5 µm−1 ring, the other values are displayed with nocorrection to illustrate the influence of the contrast suppression. Figures (a,c,e) are high-density
while (b,d,f) are low density. The measured structure factor (circles) along with the prediction (x’s)
is shown in (a,b), the scaling parameter is shown in (c,d); if the growth were perfectly describedby the model, the scaling amplitude would be one. The noise offset is shown in (e,f).

148
6.5.2

Uncorrected results

In Fig. 6.10 are the results for all rings without the contrast suppression taken into account for
any of the rings. With the exceptions mentioned in the previous paragraph, the data are identicalbetween Figures 6.9 and 6.10, but are plotted again in Fig. 6.10 or visibility. The small ratio
S0exp /S0th of Fig. 6.10c and 6.10d are an order of magnitude below unity, highlighting how little
growth is apparent if the contrast suppression factor is not accounted for. If the simplistic approachof simply evaluating the signal at tImp = 0, without regard to the subsequent evolution and without
attempting to characterize the signal as a growing component on top of a noise component, we dofind temperature behavior which qualitatively resemble the predicted behavior for majority of thedata, as shown in Fig. 6.11.
As we have seen, the observed growth of the fluctuations is difficult to reconcile with theoretical predictions for the majority of our data. Only the high-density data are in qualitative agreementwith the expected temperature dependence, and of those only the high-temperature data are inquantitative agreement. The raw data (without any consideration of the contrast suppression) fortImp = 0 are in qualitative agreement, but the growth of the fluctuations is the signature we rely
upon to distinguish “signal noise” from “noise noise” and only a small fraction of the tImp = 0 signal
appears to grow according to the predictions of Eq. 6.2. The contrast correction initially seemed
quite promising, as indicated by the good agreement between the HTHD data at k = 1.7 µm−1
and shown in Fig. 6.7 and Fig. 6.8, but applying the same method to the other data set for which
we have the calibration data did not provide the same level of agreement. After accounting for thecontrast suppression, the k = 2.5 µm−1 appears to disagree with the prediction by a similar factor
as when the correction is omitted, but in the opposite direction. In this section we discuss various
sources of error beyond those discussed in Section 6.1 which may have affected the experiment andhow we attempted to account for them.

149

(a)

(b)

(c)

(d)

(e)

(f)

Figure 6.10: Results without the correction applied to any of the wavenumbers. With the exceptionof the k = 1.7 µm−1 and k = 2.5 µm−1 ring, these values are identical to Fig. 6.9 but are provided
hear with more visibility.

150

(a)

(b)

Figure 6.11: The initial value of the structure factor signal evaluated from the fit to Eq. 6.3 without
distinguishing between the noise component and the growing component. For most values of k thetemperature dependence is qualitatively in agreement with predictions.

Figure 6.12: Low temperature High density implosions. Resolution corrections have not been applied to the signal Sraw (t). Everywherein this figure Sexp,th (0) indicates tI mp = 0 µs, the k dependence is implied. The bottom left panel’s axis labels apply to all plots, but thevertical scale differs between plots. Conditions are given in the bottom right, as is the fit model, which is identical to Eq. 6.3. The fit
parameter c0 (k) is time independent and represents the imaging noise contribution to Sraw .

151

Figure 6.13: High temperature High density. Resolution corrections have not been applied to the signal Sraw (t). Everywhere in this figure
Sexp,th (0) indicates tI mp = 0 µs, the k dependence is implied. The bottom left panel’s axis labels apply to all plots, but the vertical scalediffers between plots. Conditions are given in the bottom right, as is the fit model, which is identical to Eq. 6.3. The fit parameter c0 (k)is time independent and represents the imaging noise contribution to Sraw .

152

Figure 6.14: Low temperature Low density. Resolution corrections have not been applied to the signal Sraw (t). Everywhere in this figure
Sexp,th (0) indicates tI mp = 0 µs, the k dependence is implied. The bottom left panel’s axis labels apply to all plots, but the vertical scalediffers between plots. Conditions are given in the bottom right, as is the fit model, which is identical to Eq. 6.3. The fit parameter c0 (k)is time independent and represents the imaging noise contribution to Sraw .

153

Figure 6.15: High temperature Low density. Resolution corrections have not been applied to the signal Sraw (t). Everywhere in this figure
Sexp,th (0) indicates tI mp = 0 µs, the k dependence is implied. The bottom left panel’s axis labels apply to all plots, but the vertical scalediffers between plots. Conditions are given in the bottom right, as is the fit model, which is identical to Eq. 6.3. The fit parameter c0 (k)is time independent and represents the imaging noise contribution to Sraw .

154

155

6.6

Other potential sources of error
The factors discussed in Section 6.1 are particularly salient for the experiments of this chapter

because they represent ways the fluctuation signal seek is entangled7 with various sources of noiseor imaging limitations. The theoretical predictions were then leveraged to aid in untangling said
signal, with limited success. The theoretical prediction can only be trusted to the degree that thecalculation is performed with the correct parameters. Furthermore, if unaccounted for phenomena,such as an unexpected impurity atoms causing inelastic collisions, are present, then the theory
cannot be leveraged as described in Section 6.1.1. Such problems may arise if the conditions we
use as the inputs to our theoretical predictions are in error.
6.6.1

Miscalibrated conditions

The predictions for the growth of the structure factor depend on the density of the system,and the value at tImp = 0 depends on temperature. The number and the temperature are extracted
from time-of-flight imaging using the side imaging system, and this was checked before and after
data collection, as were the trapping frequencies and the magnetic field which determines theinitial scattering length. There are thus six quantities which determine the density profile of thecloud (N, a(B), T, ωx , ωy , ωz ). Of these, a and ωxyz are most precisely known and are unlikely tocontribute to miscalculating the density. The number and temperature are more prone to error,
and as we saw in section 5.4.3 a density error of 20% may have a significant impact. It is also
possible that the system was not in equilibrium at the time of the quench. Care was taken to
ensure that the cloud was not breathing or sloshing at the time of the quench, and the cloud wasallow to equilibrate for as long as possible to ensure the condensate was in equilibrium with thethermal component. However, the thermalization time, particularly for the two low-density data
sets, may have been insufficient. Since the low-density conditions were generated by preparing thecloud at a = 250 a0 compared to a = 110 a0 , the three-body recombination was significantly worse
7in the colloquial sense

156
for those two data sets and the equilibration time was limited to 50 ms which may not provide
enough time for superthermal density fluctuations to relax.
6.6.2

Spin purity

The spin purity of the initial or final state may also influence the dynamics. If the initialstate were pure but five percent of the atoms were not transferred to the |2, −2i state, interspecies
interactions between the |1, −1i and |2, −2i state may have had unforeseen consequences. The
pulse frequency and duration were checked before after and during each period of data collection,
but if these parameters changed during data collection, the spin purity of the final state could becompromised. Time-of-flight imaging from the side showed no detectable population remaining
in the |1, −1i state following our π pulses, suggesting the fidelity of the pulse to be >99%. If
inelastic collisions were to take place after the quench, presumably the effects would be limitedby the number of the impurity species. No systematic trend as a function of time in the data
was observed for any of the data sets, although individual experimental cycles were occasional
outliers and were removed from the analysis. Impurities in the initial state would also cause thesecomplications, but are even less likely to go unnoticed while operating the experiment as this ischecked periodically with Stern-Gerlach and spectroscopic measurements and the spin purity is
known to be indistinguishable from 100% in earlier stages of the evaporation.
6.6.3

Revisiting S0exp (k)/S0th (k)Consider Figure 6.9c and Figure 6.9d which represent a bottom-line proportional comparisonof experiment versus theory for the pre-quench static structure function S0 (k) . It is worth notingthat S(k) goes as the square of the amplitude of a density fluctuation. If we consider only calibrated
wavenumbers k = 1.6 µm−1 and k = 2.6 µm−1 of Figure 6.9c and Figure 6.9d, there are eight
distinct experiment-theory comparisons, for high/low temperature, high/low density and high/low
wavenumber. Six of those ratios are in the range 1-7, so that the agreement in density amplitude
is good to about a factor of 2.6. Only on the low-n, high-k points are the amplitudes off by as157much as a factor of size. In the ratio plots shown in Figure 6.9c and Figure 6.9d, the temperaturedependencies seem in some cases to be very far from prediction, far in excess of error bars. In
Figure 6.16a and Figure 6.16b we replot the bottom-line comparison but: (i) presenting only thecalibrated wavelengths for clarity (ii) putting the high and low density on the same vertical scale,and (iii) making the comparison on an absolute scale, not as ratios. For clarity, in four places weinclude a black bar showing the theoretically expected difference (not ratio) between high and low
temperatures, translated vertically to allow comparison with the experimental data. In this picture
the deviation from observed temperature dependence is within error bars in two cases and not sostarkly out of range for the other two. The overwhelming discrepancy is not in the temperaturedependence but that the experimental points are displaced upwards from theory by an amount

(a)

(b)

Figure 6.16: Direct comparison between experiment and theory for high-density (a) and low-density(b) conditions. The two plots share the same vertical axis. A black bar with representing thedifference between the high and low density predictions has been placed alongside the experimental
observations as a visual aid to comparing between experiment and theory.

158
which depends on the initial density of the cloud and on the wavenumber observed.
Plotting the theory-experiment comparison on an absolute scale makes more sense if we wantto understand the discrepancy as a possible additive offset. We suspect most of the difference isdue to superthermal density noise of various kinds, whether due to laser speckle in the confining
optical potential, or to as-yet not thermalized ripples leftover from the preparation. It particularly
makes sense that the low density points show greater additional noise than high density points,
because, as discussed in Section 6.6.1, the low-density points undergo and additional “last-minute”
preparation step, wherein the scattering length is made more positive so as to expand the cloud.
As the cloud expands it moves over hypothetical imperfections in the trapping potential at speeds
approaching the speed of sound. As for the k-dependence of the discrepancy, this may reflect thestructure of trap imperfections. Certainly the imperfection we are explicitly aware of has a k-vectorof magnitude 2.3 µm−1 , quite similar to the magnitude of our“high-k” ring.

Chapter 7
Summary and Outlook

The explosion of research that has followed in the wake of physicists’ ability to cool andtrap atoms has been truly remarkable, and the work done here at JILA as been a big part ofthat story. In writing this thesis and in working in the lab over the last several years, I havebeen impressed, and somewhat intimidated, by the astonishing amount of knowledge that has beenacquired by the scientific community and by how much remains unknown. Over the years, many
perplexing phenomena in ultracold atomic physics, particularly those involving strongly interacting
systems, have been explained, or sometimes explained away, by appealing to the mysterious roleof fluctuations. I am reminded of the role turbulence plays in classical fluid dynamics, incredibly
important, but extremely difficult to understand. It is the goal of our work to tackle the difficult
problem of testing our models of quantum fluctuations in bulk 3D gases to shed light on this murky
topic. It is our hope that these efforts will provide insight into these phenomena and prompt future
work which may improve upon our own.
In this thesis, we presented a vision for a fluctuation microscope which uses the instabilityof a Bose-Einstein condensate quenched to attractive interactions and the subsequent implosion
process to amplify fluctuations which are otherwise too small and subtle to see. The challenges anduncertainties which come with in-situ imaging a bulk 3D cloud, particularly when resolving finedensity structure is the goal, prevented the thorough, independent testing of the both fluctuation
spectrum of the initial cloud and the growth characteristics following the quench. Fundamentally,
we lack the precise quantitative understanding of the imaging resolution, technical noise, and160
unexpected contributions to the actual density profile to provide a thorough comparison to theory.
However, it was always understood by our group that these issues would be challenging, andwe developed a strategy to calibrate our measurements using a well known density perturbation
supplied by an optical lattice. This turned out to be incredibly valuable, as the imaging difficulties
were numerous and extremely time consuming. Only with the optical lattice data was a quantitative
test to the theoretical predictions of the structure factor possible.
The lattice experiments demonstrated quantitatively that our detection of density fluctuations was significantly less efficient than anticipated, and accounting for this phenomenon for thetwo wavenumbers for which our lattice data permitted a retroactive calibration brought clarity to
our understanding. At these wavenumbers, k = 1.7 µm−1 and k = 2.5 µm−1 , we were able to see
an effect of temperature on the static structure factor which qualitatively agrees with the theoretical prediction that structure factor grows with temperature for our high density experiments. To
arrive at this conclusion, we leveraged the assumption for the magnitude of the relative growth to
extract the initial structure factor as a free parameter. While the result agreed within error only
for the high density, high temperature data, the dependence on temperature seen for both calibrated wavenumbers provides support for the predictions of the density fluctuations of the initialsample. Without these results, one may be tempted to concluded that our observations falsified
the predictions for the structure factor.
The utility of the lattice experiments in the calibration of our imaging supports the general
notion that fluctuation measurements in-situ of a bulk 3D condensate require a thorough, quantitative understanding of one’s imaging transfer function and sources of technical noise. Subsequent
experimental work must emphasize this in order to place the experimenter in the position to robustly test each prediction for the static structure factor. In a hypothetical alternate reality whereanother year of lab work were possible, we would likely perform additional lattice experiments andattempt to develop a quantitative model of imaging resolution and noise, and certainly perform thesystematic characterization of the high resolution objective which time constraints prevented.
Theoretically, our understanding has matured from the vision presented in the introduction.

161
The fact that the structure factor is predicted to grow for all wavenumbers, even those for which
the mode is stable, was surprising to our group. The importance of the velocity/phase of theorder parameter at the time of the quench was also not appreciated until after data collection
was concluded. The phase fluctuations are larger than density fluctuations, as we determined in
Chapter 4. The role of phase fluctuations was not definitively experimentally established by thiswork. It is my hope that future work will shed light on this and related questions, and that thiswork will help others in their efforts studying these subtle quantum systems.

Bibliography[1] Milton Abramowitz, Irene A. Stegun, and Robert H. Romer. Handbook of Mathematical
Functions with Formulas, Graphs, and Mathematical Tables. American Journal of Physics,
56(10):958–958, oct 1988.
[2] U. Al Khawaja, H. T.C. Stoof, R. G. Hulet, K. E. Strecker, and G. B. Partridge. Bright
Soliton Trains of Trapped Bose-Einstein Condensates. Physical Review Letters, 89(20), 2002.
[3] J. F. Allen and A. D. Misener.
141(3558):75–75, jan 1938.

Flow of Liquid Helium II.

Nature 1938 141:3558,[4] Emine Altuntaş and I. B. Spielman. Self-Bayesian aberration removal via constraints forultracold atom microscopy. Physical Review Research, 3(4):043087, dec 2021.
[5] M. Anderlini and D. Guéry-Odelin. Thermalization in mixtures of ultracold gases. PhysicalReview A - Atomic, Molecular, and Optical Physics, 73(3):73, 2006.
[6] 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–201,
1995.
[7] 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–201,
jul 1995.
[8] A. Aspect, C. Cohen-Tannoudji, E. Arimondo, N. Vansteenkiste, and R. Kaiser. Laser cooling below the one-photon recoil energy by velocity-selective coherent population trapping:theoretical analysis. JOSA B, Vol. 6, Issue 11, pp. 2112-2124, 6(11):2112–2124, nov 1989.
[9] Vanderlei Bagnato and Daniel Kleppner. Bose-Einstein condensation in low-dimensional
traps. Physical Review A, 44(11):7439–7441, dec 1991.
[10] Vanderlei Bagnato, David E. Pritchard, and Daniel Kleppner. Bose-Einstein condensation in
an external potential. Physical Review A, 35(10):4354–4358, may 1987.
[11] J. Bardeen, L. N. Cooper, and J. R. Schrieffer. Theory of superconductivity. Physical Review,
108(5):1175–1204, dec 1957.
[12] Paul D. Beale and R. K. Pathria. Statistical Mechanics. Elsevier Science, 2011.

163[13] V L Berezinskii. Destruction of long-range order in one-dimensional and two-dimensional
systems having a continuous symmetry group I. classical systems. Sov. Phys. JETP, 32(3):493,
1971.
[14] T. P. Billam, S. L. Cornish, and S. A. Gardiner. Realizing bright-matter-wave-soliton collisions with controlled relative phase. Physical Review A - Atomic, Molecular, and Optical
Physics, 83(4):41602, 2011.
[15] T. P. Billam, S. A. Wrathmall, and S. A. Gardiner. Variational determination of approximate bright matter-wave soliton solutions in anisotropic traps. Physical Review A - Atomic,
Molecular, and Optical Physics, 85(1):13627, 2012.
[16] Sebastian Blatt. Ultracold Collisions and Fundamental Physics with Strontium. PhD thesis,
University of Colorado at Boulder, 2011.
[17] Nikolay Bogoliobov. ON THE THEORY OF SUPERFLUIDITY *. Journal of Physics,
11(1):23–32, 1947.
[18] Bose. Plancks Gesetz und Lichtquantenhypothese. Zeitschrift für Physik, 26(1):178–181, dec
1924.
[19] Eric Braaten and H. W. Hammer. Three-Body recombination into deep bound states in abose gas with large scattering length. Physical Review Letters, 87(16), 2001.
[20] C. C. Bradley, C. A. Sackett, and R. G. Hulet. Bose-einstein condensation of lithium: Observation of limited condensate number. Physical Review Letters, 78(6):985–989, 1997.
[21] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet. Evidence of Bose-Einsteincondensation in an atomic gas with attractive interactions. Physical Review Letters,
75(9):1687–1690, 1995.
[22] L D Carr and J Brand. Spontaneous Soliton Formation and Modulational Instability in
Bose-Einstein Condensates. Physical Review Letters, 92(4):040401, jan 2004.
[23] L. D. Carr, Charles W. Clark, and W. P. Reinhardt. Stationary solutions of the onedimensional nonlinear Schrodinger equation. I. Case of repulsive nonlinearity. Physical Review
A - Atomic, Molecular, and Optical Physics, 62(6):063610–063611, 2000.
[24] L. D. Carr, Charles W. Clark, and W. P. Reinhardt. Stationary solutions of the onedimensional nonlinear Schrodinger equation. II. Case of attractive nonlinearity. PhysicalReview A - Atomic, Molecular, and Optical Physics, 62(6):063611–063611, 2000.
[25] L. D. Carr, J. N. Kutz, and W. P. Reinhardt. Stability of stationary states in the cubic
nonlinear Schrödinger equation: Applications to the Bose-Einstein condensate. PhysicalReview E - Statistical Physics, Plasmas, Fluids, and Related Interdisciplinary Topics, 63(6),2001.
[26] Y. Castin and R. Dum. Bose-Einstein Condensates in Time Dependent Traps. PhysicalReview Letters, 77(27):5315–5319, dec 1996.
[27] Roman Chapurin. Precise Measurements of Few-Body Physics in Ultracold. PhD thesis, University of Colorado at Boulder, 2019.

164[28] Roman Chapurin, Xin Xie, Michael J. Van de Graaff, Jared S Popowski, José P. D’Incao,
Paul S Julienne, Jun Ye, and Eric A Cornell. Precision Test of the Limits to Universality in
Few-Body Physics. Physical Review Letters, 123(23):233402, dec 2019.
[29] Cheng An Chen and Chen Lung Hung. Observation of Universal Quench Dynamics andTownes Soliton Formation from Modulational Instability in Two-Dimensional Bose Gases.
Physical Review Letters, 125(25):250401, dec 2020.
[30] Cheng An Chen, Sergei Khlebnikov, and Chen Lung Hung. Observation of Quasiparticle
Pair Production and Quantum Entanglement in Atomic Quantum Gases Quenched to an
Attractive Interaction. Physical Review Letters, 127(6):060404, aug 2021.
[31] Qian Yong Chen, P. G. Kevrekidis, and Boris A. Malomed. Quenched dynamics of twodimensional solitary waves and vortices in the Gross-Pitaevskii equation. Journal of Optics(United Kingdom), 15(4):044012, apr 2013.
[32] Cheng Chin, Rudolf Grimm, Paul Julienne, and Eite Tiesinga. Feshbach resonances in ultracold gases. Reviews of Modern Physics, 82(2):1225–1286, apr 2010.
[33] Gui Shu Chong, Wen Hua Hai, and Qiong Tao Xie. Breathing Bright Solitons in a BoseEinstein Condensate. Chinese Physics Letters, 20(12):2098–2101, dec 2003.
[34] Steven Chu, L. Hollberg, J. E. Bjorkholm, Alex Cable, and A. Ashkin. Three-dimensional
viscous confinement and cooling of atoms by resonance radiation pressure. Physical Review
Letters, 55(1):48, jul 1985.
[35] C. Cohen-Tannoudji and J. Dalibard. Laser cooling below the Doppler limit by polarization gradients: simple theoretical models. JOSA B, Vol. 6, Issue 11, pp. 2023-2045,
6(11):2023–2045, nov 1989.
[36] Alan Corney. Atomic and Laser Spectroscopy, volume 9780199211. Oxford University Press,
oct 1977.
[37] Simon L Cornish, Sarah T Thompson, and Carl E Wieman. Formation of Bright MatterWave Solitons during the Collapse of Attractive Bose-Einstein Condensates. Physical Review
Letters, 96(17):170401, may 2006.
[38] R. Côté, A. Dalgarno, and M. J. Jamieson. Elastic scattering of two Li7 atoms. PhysicalReview A, 50(1):399–404, 1994.
[39] David Landhu Dale G. Fried, Thomas C. Killian, Lorenz Willmann. Bose-Einstein condensation of atomic hydrogen. Applied Physics B: Lasers and Optics, 69(5):357–361, 1999.
[40] F. Dalfovo, C. Minniti, S. Stringari, and L. Pitaevskii. Nonlinear dynamics of a Bose condensed gas. Physics Letters, Section A: General, Atomic and Solid State Physics, 227(34):259–264, mar 1997.
[41] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. Van Druten, D. S. Durfee, D. M. Kurn,and W. Ketterle. Bose-Einstein Condensation in a Gas of Sodium Atoms. Physical Review
Letters, 75(22):3969, nov 1995.

165[42] B. J. Da�browska-Wüster, S. Wüster, and M. J. Davis. Dynamical formation and interaction ofbright solitary waves and solitons in the collapse of Bose-Einstein condensates with attractive
interactions. New Journal of Physics, 11(5):053017, may 2009.
[43] Dirk Van Delft. History and significance of the discovery of superconductivity by Kamerlingh
Onnes in 1911. Physica C: Superconductivity, 479:30–35, sep 2012.
[44] Elizabeth A. Donley, Neil R. Claussen, Simon L. Cornish, Jacob L. Roberts, Eric A. Cornell,and Carl E. Wieman. Dynamics of collapsing and exploding Bose–Einstein condensates.
Nature 2001 412:6844, 412(6844):295–299, jul 2001.
[45] V. N. Efimov. Low-energy properties of three resonantly interacting particles. Soviet Journalof Nuclear Physics, 29(4):546, 1979.
[46] A Einstein. Quantum Theory of a Monoatomic Ideal Gas A translation of Quantentheorie
des einatomigen idealen Gases ( Einstein , 1924 ). Sitzungsberichte der Preuss. Akad. derWissenschaften, 1(3):1–5, 1924.
[47] By A Einstein. Quantum theory of the monoatomic ideal gas Second treatise. In Proceedingsof the Prussian Academy, volume 1, 1925.
[48] Jason Remington Ensher. The First Experiments with Bose-Einstein Condensation of 87 Rb.
PhD thesis, Univeristy of Colorado at Boulder, 1996.
[49] U. Fano. Effects of Configuration Interaction on Intensities and Phase Shifts. Physical Review,
124(6):1866, dec 1961.
[50] Ugo Fano. Sullo spettro di assorbimento dei gas nobili presso il limite dello spettro d’arco. Il
Nuovo Cimento, 12(3):154–161, mar 1935.
[51] Ugo Fano, Guido Pupillo, Alberto Zannoni, and Charles W. Clark. On the Absorption
Spectrum of Noble Gases at the Arc Spectrum Limit. Journal of Research of the National
Institute of Standards and Technology, 110(6):583, 2005.
[52] Herman Feshbach. Unified theory of nuclear reactions. Annals of Physics, 5(4):357–390, dec
1958.
[53] Herman Feshbach. A unified theory of nuclear reactions. II. Annals of Physics, 19(2):287–313,
aug 1962.
[54] R. P. Feynman.
Atomic theory of the λ transition in helium.
91(6):1291–1301, 1953.

Physical Review,[55] C. J. Foot. Atomic physics. Oxford University Press, 2005.
[56] Kōstas. Gavroglou. Fritz London : a scientific biography. Cambridge University Press, 1995.
[57] Rudolf Grimm, Matthias Weidemüller, and Yurii B. Ovchinnikov. Optical Dipole Traps for
Neutral Atoms. In Advances in Atomic, Molecular and Optical Physics, volume 42, pages
95–170. 2000.
[58] Walter Grimus. 100th anniversary of the Sackur-Tetrode equation. Annalen der Physik,525(3):A32–A35, mar 2013.

166[59] Li Chung Ha, Chen Lung Hung, Xibo Zhang, Ulrich Eismann, Shih Kuang Tung, andCheng Chin. Strongly interacting two-dimensional Bose gases. Physical Review Letters,
110(14):145302, apr 2013.
[60] Z. Hadzibabic and J. Dalibard. Two-dimensional Bose fluids: An atomic physics perspective.
Rivista del Nuovo Cimento, 34(6):389–434, 2011.
[61] Wenhua Hai, Chaohong Lee, and Guishu Chong. Propagation and breathing of matter-wavepacket trains. Physical Review A - Atomic, Molecular, and Optical Physics, 70(5 B), 2004.
[62] Eli Halperin. Fluctuations from in situ density measurements in a quenched BEC. (October):1–10, 2021.
[63] T. W. Hänsch, A. L. Schawlow, T. W. Hansch, A. L. Schawlow, T. W. Hansch, and A. L.
Schawlow. Cooling of gases by laser radiation. Optics Communications, 13(1):68–69, jan
1975.
[64] Charles E. Hecht. The possible superfluid behaviour of hydrogen atom gases and liquids, jan
1959.
[65] Eugene Hecht. Optics: International Edition. Pearson, 4 edition, 2003.
[66] W. R. (William Russell) Hindmarsh. Atomic spectra. 1967.
[67] P. C. Hohenberg. Existence of long-range order in one and two dimensions. Physical Review,158(2):383–386, jun 1967.
[68] P. C. Hohenberg and B. I. Halperin. Theory of dynamic critical phenomena. Reviews ofModern Physics, 49(3):435–479, jul 1977.
[69] Ming Guang Hu, Michael J. Van De Graaff, Dhruv Kedar, John P. Corson, Eric A. Cornell,and Deborah S. Jin. Bose Polarons in the Strongly Interacting Regime. Physical Review
Letters, 117(5), jul 2016.
[70] J. H. Huckans, I. B. Spielman, B. Laburthe Tolra, W. D. Phillips, and J. V. Porto. Quantum
and classical dynamics of a Bose-Einstein condensate in a large-period optical lattice. PhysicalReview A - Atomic, Molecular, and Optical Physics, 80(4):043609, oct 2009.
[71] Chen Lung Hung, Victor Gurarie, and Cheng Chin. From cosmology to cold atoms: Observation of Sakharov oscillations in a quenched atomic superfluid. Science, 341(6151):1213–1215,
sep 2013.
[72] Chen Lung Hung, Xibo Zhang, Li Chung Ha, Shih Kuang Tung, Nathan Gemelke, and Cheng
Chin. Extracting density-density correlations from in situ images of atomic quantum gases.
New Journal of Physics, 13(7):075019, jul 2011.
[73] Yu. Kagan, E. L. Surkov, and G. V. Shlyapnikov. Evolution of a Bose-condensed gas undervariations of the confining potential. Physical Review A, 54(3):R1753–R1756, sep 1996.
[74] J. W. Kane and L. P. Kadanoff. Long-range order in superfluid helium. Physical Review,155(1):80–83, mar 1967.

167[75] P. Kapitza. Viscosity of Liquid Helium below the λ-Point.
141(3558):74–74, jan 1938.

Nature 1938 141:3558,[76] Wolfgang Ketterle and N. J. van Druten. Bose-Einstein condensation of a finite number ofparticles trapped in one or three dimensions. Physical Review A - Atomic, Molecular, andOptical Physics, 54(1):656–660, 1996.
[77] L. Khaykovich, F. Schreck, G. Ferrari, T. Bourdel, J. Cubizolles, L. D. Carr, Y. Castin, and
C. Salomon. Formation of a matter-wave bright soliton. Science, 296(5571):1290–1293, may
2002.
[78] J M Kosterlitz and D J Thouless. Ordering, metastability and phase transitions in twodimensional systems. Journal of Physics C: Solid State Physics, 6(7):1181–1203, apr 1973.
[79] L Landau. Theory of the superfluidity of helium II. Physical Review, 60(4):356–358, 1941.
[80] L. Landau. On the theory of superfluidity, mar 1949.
[81] Heather Jean Lewandowski. Coherences and correlations in an ultracold Bose gas. PhD thesis, University of Colorado at Boulder, 1997.
[82] X. Li, M. Ke, B. Yan, and Y. Wang. Reduction of interference fringes in absorption imaging
of cold atom cloud using eigenface method. Chinese Optics Letters, 5(3):128–130, mar 2007.
[83] Naoto Masuhara, John M. Doyle, Jon C. Sandberg, Daniel Kleppner, Thomas J. Greytak,
Harald F. Hess, and Greg P. Kochanski. Evaporative Cooling of Spin-Polarized Atomic
Hydrogen. Physical Review Letters, 61(8):935–938, aug 1988.
[84] T. Matsubara. Quantum-Statistical Theory of Liquid Helium. Progress of Theoretical Physics,
6(5):714–730, oct 1951.
[85] N. D. Mermin. Crystalline order in two dimensions. Physical Review, 176(1):250–254, 1968.
[86] N. D. Mermin and H. Wagner. Absence of ferromagnetism or antiferromagnetism in oneor two-dimensional isotropic Heisenberg models. Physical Review Letters, 17(22):1133–1136,
nov 1966.
[87] Harold J. Metcalf and Peter van der Straten. Laser Cooling and Trapping. Graduate Texts
in Contemporary Physics. Springer New York, New York, NY, oct 1999.
[88] P. Meystre and Stig Stenholm. Introduction to feature on the mechanical effects of light.
JOSA B, Vol. 2, Issue 11, pp. 1706-1706, 2(11):1706–1706, nov 1985.
[89] C R Monroe, E A Cornell, C A Sackett, C J Myatt, and C E Wieman. Measurement of Cs-Cs
elastic scattering at T =30 µK. Physical Review Letters, 70(4):414–417, jan 1993.
[90] Linxiao Niu, Xinxin Guo, Yuan Zhan, Xuzong Chen, W. M. Liu, and Xiaoji Zhou. Optimized
fringe removal algorithm for absorption images. Applied Physics Letters, 113(14):144103, oct
2018.
[91] M. A. Olshany, S. U. Shulga, and V. G. Minogin. Laser cooling of atoms below the singlephoton classical limit. JOSA B, Vol. 6, Issue 11, pp. 2108-2111, 6(11):2108–2111, nov 1989.

168[92] H K Onnes. The Liquefaction of Helium: Communication No. 108 from the Physical Laboratory at Leiden.
1908.
[93] H K; Onnes. Research notebooks 56 and 57. Kamerlingh Onnes Archive, Boerhaave Museum,
Leiden, the Netherlands, 1911.
[94] N. G. Parker, S. L. Cornish, C. S. Adams, and A. M. Martin. Bright solitary waves andtrapped solutions in Bose–Einstein condensates with attractive interactions. Journal of
Physics B: Atomic, Molecular and Optical Physics, 40(15):3127, jul 2007.
[95] N. G. Parker, A. M. Martin, C. S. Adams, and S. L. Cornish. Bright solitary waves of trapped
atomic Bose-Einstein condensates. Physica D: Nonlinear Phenomena, 238(15):1456–1461, jul
2009.
[96] Oliver Penrose and Lars Onsager. Bose-einstein condensation and liquid helium. PhysicalReview, 104(3):576–584, 1956.
[97] Víctor M. Pérez-García, Humberto Michinel, and Henar Herrero. Bose-Einstein solitons inhighly asymmetric traps. Physical Review A - Atomic, Molecular, and Optical Physics,
57(5):3837–3842, 1998.
[98] C. J. Pethick and H. Smith. Bose–Einstein Condensation in Dilute Gases. Cambridge University Press, Cambridge, 2008.
[99] D. S. Petrov, D. M. Gangardt, and G. V. Shlyapnikov. Low-dimensional trapped gases.
Journal De Physique. IV : JP, 116:5–44, 2004.
[100] William D. Phillips and Harold Metcalf. Laser deceleration of an atomic beam. PhysicalReview Letters, 48(9):596–599, mar 1982.
[101] Lev Pitaevskii and Sandro Stringari. Bose-Einstein Condensation and Superfluidity.
Bose-Einstein Condensation and Superfluidity, mar 2016.
[102] Lev Pitaevskii and Sandro Stringari. Bose-Einstein Condensation and Superfluidity. Oxford
University Press, Oxford, 2016.
[103] Tadas Pyragius. Developing and building an absorption imaging system for Ultracold Atoms.
arXiv, sep 2012.
[104] A. Rançon, Chen Lung Hung, Cheng Chin, and K. Levin. Quench dynamics in Bose-Einsteincondensates in the presence of a bath: Theory and experiment. Physical Review A - Atomic,
Molecular, and Optical Physics, 88(3):031601, sep 2013.
[105] J. L. Roberts, N. R. Claussen, S. L. Cornish, E. A. Donley, E. A. Cornell, and C. E. Wieman.
Controlled Collapse of a Bose-Einstein Condensate. Physical Review Letters, 86(19):4211,
may 2001.
[106] P. A. Ruprecht, M. J. Holland, K. Burnett, and Mark Edwards. Time-dependent solutionof the nonlinear Schrödinger equation for Bose-condensed trapped neutral atoms. PhysicalReview A, 51(6):4704–4711, 1995.
[107] O. Sackur. Die Anwendung der kinetischen Theorie der Gase auf chemische Probleme.
Annalen der Physik, 341(15):958–980, jan 1911.

169[108] Hiroki Saito and Masahito Ueda. Bright solitons in a two-dimensional bose-einstein condensate. Journal of the Physical Society of Japan, 72(SUPPL.C):140–143, 2003.
[109] J. J. Sakurai and Jim Napolitano. Modern Quantum Mechanics. Cambridge University Press,
sep 2017.
[110] R. E. Sapiro, R. Zhang, and G. Raithel. 1-D Mott insulator transition of a Bose-Einsteincondensate. may 2008.
[111] Isaac F. Silvera and J. T.M. Walraven. Stabilization of Atomic Hydrogen at Low Temperature.
Physical Review Letters, 44(3):164–168, jan 1980.
[112] Bo Song, Chengdong He, Zejian Ren, Entong Zhao, Jeongwon Lee, and Gyu Boong Jo.
Effective statistical fringe removal algorithm for high-sensitivity imaging of ultracold atoms.
Physical Review Applied, 14(3):034006, sep 2020.
[113] Dan M. Stamper-Kurn. Peeking and poking at a new quantum fluid: Studies of gaseous
Bose-Einstein condensates in magnetic and optical traps, 2000.
[114] Kevin E. Strecker, Guthrie B. Partridge, Andrew G. Truscott, and Randall G. Hulet. Formation and propagation of matter-wave soliton trains. Nature, 417(6885):150–153, may 2002.
[115] Willian C. Stwalley and L. H. Nosanow. Possible ”new” quantum systems. Physical Review
Letters, 36(15):910–913, apr 1976.
[116] Willian C. Stwalley and L. H. Nosanow. Possible ”New” Quantum Systems. Physical Review
Letters, 36(15):910, apr 1976.
[117] C. E. Tanner, C. I. Westbrook, P. D. Lett, R. N. Watts, S. L. Rolston, and W. D. Phillips.
Optical molasses. JOSA B, Vol. 6, Issue 11, pp. 2084-2107, 6(11):2084–2107, nov 1989.
[118] H. Tetrode. Berichtigung zu meiner Arbeit: „Die chemische Konstante der Gase und daselementare Wirkungsquantum”. Annalen der Physik, 344(11):255–256, jan 1912.
[119] H. Tetrode. Die chemische Konstante der Gase und das elementare Wirkungsquantum.
Annalen der Physik, 343(7):434–442, jan 1912.
[120] Stavros Theodorakis and Andreas Hadjigeorgiou. The wavefunction of the collapsing BoseEinstein condensate. Journal of Physics B: Atomic, Molecular and Optical Physics, 50(23),2017.
[121] L. Tisza. On the theory of superfluidity. Physical Review, 75(5):885–886, mar 1949.
[122] Tobias Gerard Tiecke. Feshbach resonances in ultracold mixtures of the fermionic quantum
gases Li and 40 K, 2009.
[123] W F Vinen. The Physics of Superfluid Helium. CAS - CERN Accelerator School on
Superconductivity and Cryogenics for Accelerators and Detectors, pages 363–373, 2002.
[124] V. Weisskopf. Zur Theorie der Kopplungsbreite und der Stoßdämpfung. Zeitschrift für Physik
1932 75:5, 75(5):287–301, may 1932.

170[125] V. WEISSKOPF. The Width of Spectral Lines in Gases. Atomic Spectra, pages 328–363,
1967.
[126] V. Weisskopf and E. Wigner. Berechnung der natürlichen Linienbreite auf Grund der Diracschen Lichttheorie. Zeitschrift für Physik 1930 63:1, 63(1):54–73, jan 1930.
[127] V. WEISSKOPF and E. WIGNER. Calculation of the Natural Line Width on the Basis of
Dirac’s Theory of Light. Atomic Spectra, pages 304–327, 1967.
[128] D. J. Wineland and Wayne M. Itano. Laser cooling of atoms. Physical Review A, 20(4):1521,
oct 1979.
[129] D.J. Wineland and H Dehmelt. Proposed 1014∆ν < ν Laser Fluorescence Spectroscopy on
Tl + Mono-Ion Oscillator. Bulletin of the American Physical Society, (20):637, 1975.
[130] David Wood. The computation of polylogarithms. 0:1–19, 1992.
[131] Xin Xie. Precise Calibrations of Few-Body Physics in Potassium-39: Experiment and Theory.
PhD thesis, University of Colorado at Boulder, 2020.
[132] Xin Xie, Michael J Van De Graaff, Roman Chapurin, Matthew D Frye, Jeremy M Hutson,
José P D’incao, Paul S Julienne, Jun Ye, and Eric A Cornell. Observation of Efimov Universality across a Nonuniversal Feshbach Resonance in 39 K. PHYSICAL REVIEW LETTERS,
125:243401, 2020.
[133] Feng Xiong, Yun Long, and Colin V. Parker. Enhanced principle component method forfringe removal in cold atom images. Journal of the Optical Society of America B, 37(7):2041,
jul 2020.
[134] I Yavin, M Weel, A Andreyuk, and A Kumarakrishnan. A calculation of the time-of-flight
distribution of trapped atoms. American Journal of Physics, 70(2):149–152, 2002.
[135] Clarence Zener. Non-adiabatic crossing of energy levels. Proceedings of the Royal
Society of London. Series A, Containing Papers of a Mathematical and Physical Character,
137(833):696–702, sep 1932.

Appendix A
Density of States

A.1

Density of states for a free particle
the volume on an d-ball and the surface area of the corresponding d-sphere in d dimensions

is given by
Vd (R) =π d/2RdΓ(d/2 + 1)

Sd (R) =2π d/2 d−1RΓ(d/2)(A.1)

p
A free particle has energy p = 2m
hence
2

dp
=
d

r

m
.
2

In general the density of states for a system isg() ≡dd
d

where d is the total number of states with energy less than . The phase space associated with aparticle with d degrees of motional freedom confined to a box volume V is
Z
V=

d

Z ∞

d r
V

dd p.

0

One could explicitly solve the Schrodinger equation in d dimensions to proceed, but here we will
use the fact that there is on average one quantum state per volume (2π~)d in phase space.
As a brief aside, I learned recently that this value did not require quantum mechanics in its
modern (i.e. post 1925) form but was in fact determined by Otto Sackur and Hugo Tetrode in1912 [107, 119, 118, 58]. Their method relied on deriving the entropy of a monoatomic ideal gas

172
under the assumption that elementary volume of phase space was zhd for a system with d degrees
of freedom. Classical physics does not provide an a priori justification for this ansatz, but in theircomparison to experimental data then available for mercury vapor, they found agreement with thequantization of phase space. The descretization of pohase space into elementary cells is a top-down
method of calculating the actual number of microstates available to a system; later methods using
standarad quantum mechanics provides a bottum-up approach. Tetrode actually fit the constantz and determined it to be consistent with unity, prompting him to propose it may be precisely
unity. I find it interesting how much this approach reflected Planck and Einstein’s own approaches
in phenomenologically determining the necessity, utility, and numerical value of Planck’s constanth, particularly considering Sackur and Tetrode applied h in the problem of massive particles as
opposed to photons. In a paper celebrating the 100th anniversary of the Sackur-Tetrode equation
for the absolute entropy, Grimus repeats Tetrode’s fit to mercury data and finds it accurately
determines h to within 1%.
The total number of states N with energy E <  and thus with momentum p < p =√2mis thus
V()
V
N =
=
d
(2π~)
(2π~)d

Z p
0

V
d p=
(2π~)d
d

Z 

d

r

d 
0

m
2

(A.2)

We compute the integral using Eq. ?? for the volume in p-space.
N =

√
V
V (2πm)d/2
V
(
2m
)
=
d
(2π~)d(2π~)d Γ(d/2 + 1)

(A.3)

and thus the density of states (per unit volume) is
D() =√
1d(2πm)d/2
V
(
2m
)
=
d
(2π~)d2(2π~)d Γ(d/2 + 1)

and usingΓ(m + 1/2) =(2m)! √
π
4m m!

(A.4)

173
we can see that

m3/2 √
D3U () = √
2 π 2 ~3m
D2U () =
2π~2
√
m 1
√
D1U () = √2 π~ 

(A.5)

where we have used the subscript U to denote the system is uniform.

A.2

Density of States for harmonic confinement
For the harmonic oscillator in d dimensions, the eigenstates have energies
{n} = ~X

ωi n iwhere {n} denotes the set of occupation numbers {n1 , ..., nd } for each dimension and we are ignoring
the zero-point energy. We make the assumption that all the ωi are equal to proceed under thePsimplifying relation nT =ni is the total occupation of all modes. We determine the numberof states N with energy E <  by counting the available combinations of the ni . If d = 1 thenN1D = nT + 1 since n1 may vary from 0 to nT . Then for d = 2 and fixed n ≤ nT we have n + 1
options for n1 and n + 1 − n1 options for n2 , making the multplicity
N2D =

n=n
1 =n
XT nX
n=0 n1 =0

n + 1 − n1 =

n1X
=nT

n2

n1 =0

Since n1 and n2 have exactly the same set of options (as we could have fixed n2 instead) we canPthT
replace n2 → n1 in nn11 =n=0 n2 and thus N2D = Tn where Tn is the n triangular number. Adding
a third dimension and considering how for each unique choice of n3 the multiplicity of states is
TnT −n3 , and recognizing our sums for each dimension are identical, the multiplicity of the 3D QHOwith total occupancy less than or equal to nT is the sum of the first nT triangular numbers, also
known as the tetrahedral number T en . This can be generalized to high dimensions.
We are interested in d = 1, 2, 3 only. For these cases the density of states in a harmonic

174
oscillator is

1
~ω

(A.6)
D2H () = 2 2~ ω2
D3H () = 3 3
2~ ω
where the subscript H indicates harmonic confinement. Although these results were derived forD1H () =isotropic potentials, they generalize to anisotropic traps with the substitution ω 2 → ωx ωy ωz .

Appendix B
Lattice potentials

B.1

Interference of two laser fields
Optical lattices are the result of interference of optical fields. While optical interference isa general phenomenon, in the context of optical lattices as used in AMO physics typically theinterfering fields are coherent. In our laboratory each lattice is produced by two laser fields if thesame wavelength. We may therefore write the two fields as
E1 = E01 cos(k1 · r − ωt)(B.1)

E2 = E02 cos(k2 · r − ωt)where ω is the frequency of the laser and |k1 | = |k2 | = ωc are the lasers’ wavevectors of equalmagnitude which differ in direction. The polarization and magnitude of the fields are determinedby E0,i . When the polarization are identical, we will get an intensity of the form [65]:I = I1 + I2 + 2

(B.2)

p
I1 I2 cos δ

Where δ = (k1 − k2 ) · r and we have switched from electric vector fields to scalar intensities. In
the case of equal intensities I1 = I2 = I0 we get
I = 2I0 (1 + cos δ) = 4I0 cos2

δ
2

(B.3)

and we can see that we have maxima whenever δ = 2πn for integer n. Defining the lattice spacing as
2π
2πd, the maxima occur when 2πn = 2π
d r. So if k1 = λ (cos θx̂+sin θŷ) and k2 = λ (− cos θx̂+sin θŷ)176
in terms of the laser wavelength λ then we see that
δ = (k1 − k2 ) · r =4π
x cos(θ).
λ

Hence we can see that the lattice spacing is given byd=λ.
2 cos θ(B.4)

Now let us define the lattice wavevector kL as the half the magnitude of the reciprocal lattice
wavevector:
kL =

δ
π
2π sin(θ)= ŷ =
.
2
d
λ

(B.5)

The single photon recoil (such as from an off-resonant excitation of the atom by the lattice beams)
is determined by Eγrec = ~2 k02 /2m where k0 = 2π/λ. The subscript L denotes the lattice while thesubscript 0 denotes the light field. The lattice recoil energy is given by(B.6)

ELrec = ~2 kL2 /2m

The lattice recoil energy ELrec and the photon recoil Eγrec are the same in the case of counterpropagating beams, for which θ = π/2. Future reference to the lattice recoil ELrec will omit the recsubscript for concision.

B.2

Partial derivation of Raman-Nath approximation
The generating function for Bessel function of the first kind is (9.4.41 of [1])
 

∞
X
z
1
exp
t−
=
tn Jn (z),
2
t
n=−∞

(B.7)

a form amenable to the exponential representation of sin(θ), thus yielding:ia sin(θ)

e

=e

a iθ
(e −e−iθ )
2

=

∞
X

einθ Jn (a).

(B.8)

n=−∞

With the substitution θ → θ + π/2 we get the Jacobi-Anger expansione

ia cos(θ)

=

∞
X
n=−∞

in Jn (a)einθ

(B.9)

Appendix C
Thermal distribution of Bosons

A Maxwell-Boltzmann ensemble in a harmonic trap has the form
n(x, y, z) =Nth ωx ωy ωz m3/2 2k−mT (x2 ωx2 +y2 ωy2 +z 2 ωz2 )e b(2πkb T )3/2

=n0,th,M B,3D e



2
2
2
− x2+ y2+ z 2
2σx

2σy

(C.1)

2σz

N
q thp
=p
e
2σx2 2σy2 2σz2 π 3/2



2
2
2
− x2+ y2+ z 2
2σx

2σy

2σz

where in the second line n0,th,M B,3D is the peak (n0 ) (ie at the center), thermal (nth ), MaxwellBoltzmann (nM B ), three-dimensional (n3D ) density:n0,th,M B,3D =Nth ωx ωy ωz m3/2.
(2πkb T )3/2

(C.2)

The third line of C.1 is deliberately written in a manner emphasizing that the 2 often seen in
normalized Gaussian prefactors is due to the σ terms, not the π. So considering that the RHS ofP
zn
C.1 is given as the argument of the Bose function gν (z) = ∞n=1 nν , we can thus think of the Bosefunction as a series of Gaussians with there size and amplitude rescaled in particular ways. How
precisely? Let’s simplify to one dimension and consider the base Gaussian:G1 (x) = √

1
2σ 2

−x2

−x2

e 2σ2 = Ae 2σ2 .

(C.3)

This is the argument of gν (z), so the first term is of course unmodified, but the nth term gives aprefactor of 1/nν and transforms the exponential like:
e

−x2
2σ 2

⇒e

−nx2
2σ 2

=e

−x2
√
2(σ/ n )2178
Thus the nth term of the series of Gaussians is
A 2(σ/−x√2n )2Gn (x) = ν e
.
n

(C.4)

√
This is a new Gaussian with σnthterm = σ/ n , but since we are discussing density to get√number we integrate over all space, which gives a factor of 1/ n (per dimension!) from the rescaled√
R −x2
Gaussian since e 2σ2 = σ 2π and a factor of 1/nν from the denominator in ??. More explicitly:Z ∞Gn (r)dd r =

−∞

1

Z ∞

nν+d/2−∞G1 (r)dd r.

(C.5)

Thus, for three dimension, ν = 3/2 so the total normalization isn0,th,BE,3D = n0,th,M B,3D

∞
X
1
n=1

C.0.1

n3

.

(C.6)

Time of flight analysis

Following [134] we can derive the density profile of a thermal cloud obeying Maxwell-Boltzmann
statistics for simplicity. The probability column density of the cloud is given by


x2
z2
p2x
p2z
n(r0 , p0 )dx0dz0dpx0 dpz0 = C exp − 2 − 2 − 2 − 2 dr02 dp20 .
2σx 2σz
2px0 2pz0(C.7)

Let us suppose the cloud is in thermal equilibrium so px0 = py0 = pz0 = p0 and that the constant C fixes the normalization. Now lets consider where the atoms will be after a time t andchange n(r0 , p0 ) = n(x0 z0 , p0 , p0 ) into n(x0 , z0, x(t), z(t)). A particle initially at r0 = (x0 , z0 ) with
momentum p0 = (px0 , pz0 ) will move after a time t such that:
px0
t
m
pz0
z = rz0 +
t
m

(C.8)

x − rx0mt
x − rx0
px0 =
m
t

(C.9)

x = rx0 +where we ignore gravity. Inverting yields:
px0 =

179
We can use this to define the change of coordinates via
ZZ

ZZ
f (a, b)dadb ⇒f (c, d)(C.10)

To simplify, consider the initial probability density of the cloudmω 2 x20
p2
+
2
2m
h
i
x2
The position-space form of the cloud is thus Gaussian n(x) ∝ exp − 2σ02 with1
n(x0 , p) = C exp
kb T

1
σ0 =
ω



r

kb T
.
m

(C.11)(C.12)

When the cloud is released, the particles move ballistically. We measure the position of the particles
0)
after a time t at a position x, and thus we know it’s momentum must be p = m(x−xSo integratingtpt
over all initial positions requires a delta function to enforce momentum conservation δ(x−(x0 + m)).

The delta function in p space is thusptδ(x − (x0 + m)) = mt δ(p − mt (x − x0 ))and then the position after the time of flight is done by integrating over the initial probability
density so constrained:
ZZ
ntof (x, t) = C



−1
exp
kb T



mω 2 x20
p2
+
2
2m



t
m
m δ(p − t (x − x0 ))dx0 dp.

(C.13)

Applying the delta function and simplifying yields:

Z




t
−mx2
−m
1
2x
2
2
ntof (x, t) = C exp
exp
ω + 2 x0 − 2 x0
m
2kb T
2kb T
t
t
(C.14)√

2πkb T m
mx2
ω2
ntof (x, t) = C √exp −.
2kb T 1 + ω 2 t2
1 + ω 2 t2p
Hence we see that the Gaussian cloud with initial rms size σ0 = kb T /m /ω expands and remains√
Gaussian but with σ(t) = σ0 1 + ω 2 t2This is convenient because although an ideal, thermal Bose gas is not Gaussian, but instead
has a density distribution given by the first term of Eq. 3.28, which is the Bose function g3/2 ,180
the argument of the Bose function gν (z) =zn
n=1 nν , is itself a Gaussian.

P∞

So every term of thatfunction has the same form as:ZZ
ntof (x, t) = C
C.0.1.1

exp

−x20
−p2
+
2σ 2
2mkb Ttmm δ(p − t (x − x0 ))dx0 dp(C.15)

Thomas Fermi Expansion

The parabolic shape of a condensate in the Thomas Fermi regime also expands by a simple
rescaling of its widths, conveniently retaining its original functional form [73, 26, 40]

ProQuest Number: 29321798
INFORMATION TO ALL USERS
The quality and completeness of this reproduction is dependent on the quality
and completeness of the copy made available to ProQuest.

Distributed by ProQuest LLC ( 2022 ).
Copyright of the Dissertation is held by the Author unless otherwise noted.

This work may be used in accordance with the terms of the Creative Commons licenseor other rights statement, as indicated in the copyright statement or in the metadata
associated with this work. Unless otherwise specified in the copyright statement
or the metadata, all rights are reserved by the copyright holder.

This work is protected against unauthorized copying under Title 17,
United States Code and other applicable copyright laws.

Microform Edition where available © ProQuest LLC. No reproduction or digitizationof the Microform Edition is authorized without permission of ProQuest LLC.

ProQuest LLC
789 East Eisenhower Parkway
P.O. Box 1346
Ann Arbor, MI 48106 - 1346 USA