ABSTRACT

Title of Dissertation: TURBULENCE AND SUPERFLUIDITY
IN THE ATOMIC BOSE-EINSTEIN CONDENSATE

Mingshu Zhao
Doctor of Philosophy, 2024

Dissertation Directed by: Professor Ian Spielman
Joint Quantum Institute,
National Institute of Standards and Technology,
and Department of Physics, University of Maryland College Park

In this dissertation I investigate turbulence in atomic Bose-Einstein condensates (BECs),

focusing on the challenge of quantifying velocity field measurements in quantum fluids. Turbulence,

a universal phenomenon observed across various scales and mediums – from classical systems

like Earth’s oceans and atmosphere to quantum fluids including neutron stars, superfluid helium,

and atomic BECs – exhibits complex fluid motion patterns spanning a wide range of length

scales. While classical turbulence has been extensively studied, quantum systems present many

open questions, particularly regarding the existence of an inertial scale and the applicability of

Kolmogorov scaling laws.

I introduce a novel velocimetry technique, analogous to particle image velocimetry (PIV),

using spinor impurities as tracer particles. This method enables the direct measurement of the

velocity field and thereby the velocity structure functions (VSFs) in turbulent atomic BECs.

The technique overcomes limitations of existing experimental approaches that rely on time of



flight (TOF) measurements, offering a clearer connection to VSFs and enabling a more direct

comparison of turbulence in atomic gases with other fluids.

The cold-atom PIV technique enables directly measuring the velocity field, leading to a

detailed analysis of both VSFs and the velocity increment probability density functions (VI-

PDF). Key findings include the observation of superfluid turbulence conforming to Kolmogorov

theory from VSFs, and intermittency from high order of VSFs and the non-Gaussian fat tail in

the VI-PDF.



TURBULENCE AND SUPERFLUIDITY IN THE ATOMIC
BOSE-EINSTEIN CONDENSATE

by

Mingshu Zhao

Dissertation submitted to the Faculty of the Graduate School of the
University of Maryland, College Park in partial fulfillment

of the requirements for the degree of
Doctor of Philosophy

2024

Advisory Committee:
Professor Daniel Lathrop, Chair
Professor Ian Spielman, Advisor
Professor Nathan Schine
Professor Thomas Antonsen
Professor Johan Larsson



© Copyright by
Mingshu Zhao

2024



Preface

In this thesis, I aim to present a comprehensive and accessible exploration of turbulence in

atomic gases, particularly through the lens of cold atom Particle Image Velocimetry (PIV).

The inception of this technique came from an unlikely source: a longstanding issue with

the bipolar current servo used for magnetic field control in our lab. This defect, when combined

with high-power radiofrequency (RF) signals, led to an uncontrolled DC offset, impacting the

lab’s RF controls. My advisor, Ian, once suggested using two-tone copropagating Raman beams

as an RF alternative, which, while effective, was not widely adopted.

The breakthrough came when I considered the use of a spatial light modulator (SLM) to

control the spatial pattern of the Raman beam. This would allow for spatially-dependent RF

coupling and, consequently, spatially-dependent spinor transfers. This concept paralleled the

principles of Particle Imaging Velocimetry (PIV) in classical fluids, with spatially-dependent

spinor transfers acting as tracer particles.

I am hopeful that this innovative technique will reinvigorate interest in hydrodynamic

studies within the experimental physics community. By enabling direct measurements of velocity

fields in cold atom experiments, we open new doors to understanding and exploring the rich

dynamics of quantum fluids.

ii



Acknowledgments

I want to start by thanking my advisor, Ian Spielman, for his support and encouragement,

especially when I was exploring the idea of measuring velocity fields, even though it was a tough

challenge.

A big thanks goes to Alessandro Restelli for helping me build the servo and all the other

electronic parts we needed. My lab mate, Junheng Tao, deserves a special mention for all the

hard work in setting up our lab and making our equipment work better.

I’m also grateful to Dimitris Trypogeorgos, Ana Valdés Curiel, and Qiyu Liang for showing

me how to use the RbLi apparatus. Ana was really helpful with fixing our coils and getting our

optics and lasers right. Qiyu helped a lot with our microscope system and inspired us with her

drive. Thanks to Francisco Salces Cárcoba and Chris Billington for designing our new setup and

teaching us about baking the system.

I want to acknowledge my friends from other labs in the PSC basement - Peter Elgee,

Ananya Sitaram, Yanda Geng, Shouvik Mukherjee, Hector Sosa Martı́nez, Monica Gutierrez

Galan, Swarnav Banik, Madison Anderson, Sarthak Subhankar, Tsz-Chun Tsui, Kevin Weber,

Patrick Banner, Deniz Kurdak, Yaxin Li and James Maslek - for lending me equipment and

giving advice whenever I was stuck.

A shoutout to the folks from NIST - Graham Reid, Alina Pineiro, Amilson Fritsch, Mingwu

Lu and Emmanuel Mercdo, Emine Altuntas, Shangjie Guo, Yuchen Yue, Micheal Doris, Dario

iii



D’Amato, João Braz and Davis Garwood - for their helpful advice in our group meetings.

Mental health is really important, so I’m incredibly thankful for my friends Zhiyu Yin,

Yihang Wang, Shuyang Wang, Xiangyu Li, Yalun Yu, Wance Wang, Haonan Xiong, Tianyu Li,

Ziyue Zou, Gong Cheng, Yixu Wang, Jingnan Cai, Haoying Dai, Haining Pan and Kaixin Huang

and others for being there for me.

And lastly, none of this would have been possible without the support from my family.

iv



Table of Contents

Preface ii

Acknowledgements iii

Table of Contents v

List of Tables viii

List of Figures ix

List of Abbreviations xiv

Chapter 1: Introduction 1
1.1 Thesis overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

Chapter 2: Atom Light Interaction and Bose-Einstein Condensate 6
2.1 Two-level systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.2 Coherent transfer between two levels . . . . . . . . . . . . . . . . . . . . . . . . 8

2.2.1 Rotating Wave Approximation . . . . . . . . . . . . . . . . . . . . . . . 9
2.2.2 Rabi Oscillation and Adiabatic Rapid Passage . . . . . . . . . . . . . . . 9

2.3 Optical Stark shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.3.1 Tune-out (magic) Wavelength . . . . . . . . . . . . . . . . . . . . . . . 16
2.3.2 AC vector light shift as RF . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.4 Laser Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.4.1 Optical Molasses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.4.2 Magneto-optical Trap . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.4.3 Polarization gradient cooling . . . . . . . . . . . . . . . . . . . . . . . . 25
2.4.4 Evaporative Cooling . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.5 Atomic Bose-Einstein Condensates . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5.1 Gross-Pitevskii equation . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5.2 Bogoliubov-de-Gennes excitations . . . . . . . . . . . . . . . . . . . . . 32

Chapter 3: Review of Turbulence and Velocimetry 34
3.1 Incompressible Turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.1.1 Basic features of turbulence . . . . . . . . . . . . . . . . . . . . . . . . 36
3.1.2 The Reynolds-Averaged Navier–Stokes Equations . . . . . . . . . . . . . 39
3.1.3 Vortex Stretching . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

v



3.1.4 Equilibrium range theories . . . . . . . . . . . . . . . . . . . . . . . . . 41
3.2 Compressible turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

3.2.1 Governing dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
3.2.2 Coarse graining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.2.3 Scale Decomposition . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53

3.3 Quantum turbulence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.3.1 Superfluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
3.3.2 Theoretical models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.3.3 Coarse Graining and Richardson Energy cascade . . . . . . . . . . . . . 58
3.3.4 Dissipation mechanism . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.3.5 Quantum Turbulence Experiments . . . . . . . . . . . . . . . . . . . . . 61

3.4 Velocimetry . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.1 Velocimtery in classical fluid . . . . . . . . . . . . . . . . . . . . . . . . 64
3.4.2 Velocimetry in cold gases . . . . . . . . . . . . . . . . . . . . . . . . . . 67

Chapter 4: The RbRb apparatus 71
4.1 Experimental Layout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
4.2 Vacuum system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4.2.1 Vacuum bakeout . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
4.2.2 Atom source control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

4.3 Experimental optical setups . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3.1 MOT optics geometry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.3.2 Laser Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.3.3 Laser locking . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

4.4 Magnetic field Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.4.1 Coils in RbRb . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.4.2 Coil winding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.4.3 Current servos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.4.4 Magnetic transport . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91
4.4.5 Bias and gradient field control . . . . . . . . . . . . . . . . . . . . . . . 93

4.5 Digital micromirror device control . . . . . . . . . . . . . . . . . . . . . . . . . 96

Chapter 5: Velocimetry 98
5.1 PIV in BEC . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

5.1.1 Details of the Raman PIV technique . . . . . . . . . . . . . . . . . . . . 100
5.2 Velocity field measurement of benchmark flow in BEC . . . . . . . . . . . . . . 102

5.2.1 Dipole mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103
5.2.2 Scissors mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105
5.2.3 Rotating Trap . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

Chapter 6: Turbulence experiments in atomic BEC 112
6.1 Turbulence Generation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
6.2 Velocity Increments Dataset Collection . . . . . . . . . . . . . . . . . . . . . . . 113
6.3 Structure functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

6.3.1 Third order structure function (S3) . . . . . . . . . . . . . . . . . . . . . 115

vi



6.3.2 Second order of structure function (S2) . . . . . . . . . . . . . . . . . . . 116
6.3.3 Sn(l) for n ≤ 7 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120

6.4 Velocity increments PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 122
6.4.1 PDF from Deconvolution . . . . . . . . . . . . . . . . . . . . . . . . . . 123
6.4.2 Non-Gaussian Statistics . . . . . . . . . . . . . . . . . . . . . . . . . . . 127
6.4.3 Structure Function after Deconvolution . . . . . . . . . . . . . . . . . . 129

6.5 Numerical Simulation of the stirring turbulence . . . . . . . . . . . . . . . . . . 129
6.5.1 Numerical dissipative GPE scheme . . . . . . . . . . . . . . . . . . . . . 131
6.5.2 Kinetic energy spectrum and structure functions . . . . . . . . . . . . . . 132

Chapter 7: Measurement of Superfluid Density 136
7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 136
7.2 Leggett’s Formula in anisotropic superfluids . . . . . . . . . . . . . . . . . . . . 138
7.3 Superfluid Hydrodynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 141

7.3.1 Hamiltonian fluid mechanics . . . . . . . . . . . . . . . . . . . . . . . . 141
7.3.2 Coarse Graining . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 143
7.3.3 Coarse-grained superfluid hydrodynamics . . . . . . . . . . . . . . . . . 146

7.4 Superfluid sum-rule and sound velocity . . . . . . . . . . . . . . . . . . . . . . . 157
7.4.1 Josephson sum-rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 157
7.4.2 Speed of sound from the superfluid hydrodynamics . . . . . . . . . . . . 160

7.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
7.5.1 Bragg Spectroscopy . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 161
7.5.2 Scissors mode . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 164
7.5.3 Rotation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 167

Chapter 8: Conclusion and Outlook 169
8.1 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 169
8.2 Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 170

Appendix A: List of Publications 171
A.1 List of Peer-Reviewed Publications . . . . . . . . . . . . . . . . . . . . . . . . . 171
A.2 Manuscripts in Preparation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 171

Bibliography 173

vii



List of Tables

6.1 Velocity increments dataset. The dataset covers six distinct tracer separations
from the squared array pattern, i.e., (l1, l2, l3, l4, l5, l6) = (10.6, 11.4, 12.6, 15.0, 16.1, 17.8)µm.
Velocity increments data is labeled as δ[i](j)vnj

, where i = 1..44 is the index
of the experimental run, and j = 1..6 is the index of the tracer separation lj .
nj = 1.. ≈ 200 for j = 1, 2, 3 and nj = 1.. ≈ 100 for j = 4, 5, 6 since in each
experimental run, the measurement is repeated for ≈ 50 times for each tracer
pattern which consists 4 data points for the square side (j = 1, 2, 3) and 2 for the
square diagonal (j = 4, 5, 6). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 114

viii



List of Figures

1.1 (a) An example of velocity field. (b) Tracer particles (bright pink) are injected
into the system at t = 0. (c) At t = ∆t, tracer particles (bright pink) move from
their initial position (dark pink). . . . . . . . . . . . . . . . . . . . . . . . . . . 3

2.1 87Rb D line. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.2 Frictional force of Doppler cooling when s(0) = 2 and δ = Γ. The blue (purple)

dashed line is when k and v have the same (opposite) direction, and the red solid
line is the sum of the two. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23

2.3 MOT scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
2.4 Sisyphus cooling scheme. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27

3.1 A representation of the chaotic and random nature of the velocity field over time. 37
3.2 Illustration of the light dragging effect in a moving medium with vertical motion. 68

4.1 The RbRb apparatus. Reproduced from [1]. For clarity the science cell optics is
not shown here, and only half of the transport coils are shown. . . . . . . . . . . 71

4.2 The RbRb vacuum chamber. (a), (b) are reproduced from [1]. (c) shows the
physical vacuum chamber. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73

4.3 (a) shows the bakeout technique we use, in which the bellow is connected to a
turbo pump, and the vacuum system is wrapped by the UHV aluminum foils. (b)
shows a typical bakeout log. The current reading is from two ion pumps and the
spikes are temperature sudden increasing of the baking during the process. . . . . 75

4.4 MOT optics. Reproduced from [1]. (a) Three pairs of beams that form the MOT
region within the glass cell. The pair that goes in and out of the plane is not
labeled in (a). In (b), MOT5 is labeled, and MOT6 is hidden behind the coils.
In (b) an additional optical pumping beam labeled by OptPump has orthogonal
polarization to the MOT5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

4.5 Laser frequency diagram. The unit is in MHz. The first row shows the frequency
rough difference of the cooling, master, and repump laser. The second row gives
detailed beam frequencies relative to the lock. . . . . . . . . . . . . . . . . . . . 78

4.6 Master Laser SatAbs configuration. . . . . . . . . . . . . . . . . . . . . . . . . . 79
4.7 Cooling laser configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81
4.8 Repump laser configuration. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.9 Master laser locking diagram . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83

ix



4.10 Beatnote lock diagram. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
4.11 Coils in RbRb. Reproduced from [1]. . . . . . . . . . . . . . . . . . . . . . . . . 86
4.12 Coils winding. (a) shows the coil winding pipeline. (b) shows the coils with

epoxy in a vacuum bubble. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
4.13 (a) Coil surface after lathing under microscope. (b) A cloverleaf coil with epoxy

covered. (c) A round shape coil with epoxy removed. . . . . . . . . . . . . . . . 88
4.14 (a) The metal bars are the ”bus” of the +15V , ground and −15V for the current

control in the lab. Below is the transistor bank for the transport coil current. (b)
Shows the 11 high current cable connectors for the transport coil. (c) Shows the
diagram of the unipolar and bipolar current control in the lab. . . . . . . . . . . 88

4.15 Bipolar servo schematics. (a) High-level modules including current sensing,
feedback and current generation. (b) Bipolar current control schematics. The
control signal is input on the left and divided into positive and negative channels
(top and bottom respectively) that control separate banks of NMOS and PMOS
transistors before being delivered to the load (far right) and sensed. . . . . . . . . 90

4.16 Time traces of current control during a round-trip magnetic transport. . . . . . . . 93
4.17 Bias coil and current configuration. . . . . . . . . . . . . . . . . . . . . . . . . . 94
4.18 Gradient coil and current configuration. . . . . . . . . . . . . . . . . . . . . . . 95
4.19 Calibration of B field using microwave resonance. Blue points are measured from

ARP resonance, and the orange is fitted from the square root of a parabola. . . . . 96

5.1 Concept. (a) An example of velocity field. (b) Tracer particles (bright pink)
are injected into the system at t = 0. (c) At t = ∆t, tracer particles (bright
pink) move from their initial position (dark pink). (d) Spatially-resolved Raman
technique to create localized tracer particles. . . . . . . . . . . . . . . . . . . . 101

5.2 Velocity measurement of the dipole mode using PIV. (a) shows the ground state
BEC in a harmonic trap. (b) shows the PIV tracer patterns. Colorbars in (a-b)
show the optical density (OD). (c) shows the position change of a tracer under
dipole mode excitation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

5.3 Scissors mode measurement using PIV. (a) shows the velocity field of a rotating
harmonic trap. In the PIV, the red dots move to the pink dots with the angle
between two ”arms” θ changing dynamically. (b) shows the scissors mode oscillation
of the θ using PIV, and the inlet shows the tracers’ pattern used in the experiment. 106

5.4 Rotating harmonic trap flow field measurement using PIV. (a) shows measured
the velocity field of a rotating harmonic trap. The inlet shows the tracer particles
motion of the black arrow. The arrow end is the tracer’s initial position and
the arrow head is the final position after 1.5ms. The color scale shows the
experimental optical density from PTAI. (b) shows the velocity field from the
GPE simulation under the same case as (a). The color scale is rescaled to agree
with the measurement in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110

6.1 Turbulence initialization. Atomic density measured at at six times during the
excitation process. Two counter-rotating stirring rods are moved in the condensates.
The color bar is the optical density from PTAI measurements. . . . . . . . . . . . 113

x



6.2 An example of Tracer patterns. Tracers are positioned on a squared array, each
side measuring 12.6µm. The left shows the position where the tracers are originally
injected and the right represents the tracers’ position with a time interval ∆t =
0.3ms compared to the left. The red squares outline each region of interest
(ROI), within which the red dots indicate the center of mass used for velocity
field calculation. The color bar is the optical density from the PTAI measurements. 114

6.3 Measured S3(l). Data point results from the average of 44 experimental runs,
each of which derived S3(l) from about 50 nominally identical experimental
repetitions. The uncertainties are the two-sigma standard error of the mean across
the set of experimental runs, and lines are fitted to the data plotted along with their
2− σ uncertainty band. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 115

6.4 Measured S2(l). Data point results from the average of 44 experimental runs,
each of which derived S2(l) from about 50 nominally identical experimental
repetitions. The uncertainties are the two-sigma standard error of the mean across
the set of experimental runs, and lines are fitted to the data by l2/3 plotted along
with their 2− σ uncertainty band. . . . . . . . . . . . . . . . . . . . . . . . . . 119

6.5 Measured SS
n (l). (a) Log-Log plot of SS

n (l), n = 1..7 fitted to the dashed lines
anl

εn . The error bar shows the two-sigma standard error of the mean. (b) intermittency
correction εn−n

3
versus n. The error bar is the two-sigma uncertainty of the fitting

in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 120
6.6 Measured SL

n (l). (a) Log-Log plot of SS
n (l), n = 1..7 fitted to the dashed lines

anl
εn . The error bar shows the two-sigma standard error of the mean. (b) intermittency

correction εn−n
3

versus n. The error bar is the two-sigma uncertainty of the fitting
in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

6.7 Measured ST
n (l). (a) Log-Log plot of SS

n (l), n = 1..7 fitted to the dashed lines
anl

εn . The error bar shows the two-sigma standard error of the mean. (b) intermittency
correction εn−n

3
versus n. The error bar is the two-sigma uncertainty of the fitting

in (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121
6.8 Histograms of longitudinal velocity increments at l = 10.6µm for unstirred

(blue) and stirred (red) Bose-Einstein condensates (BECs). Panel (a) shows normalized
PDFs in linear scale, while panel (b) displays PDFs rescaled to peak at 1 in log
scale. Error bars indicate

√
n statistical counting errors. . . . . . . . . . . . . . . 124

6.9 Histograms of longitudinal velocity increments at l = 10.6µm for stirred BECs
(red) and their deconvolved versions (light blue), with the red curve illustrating
the convolution of instrumental noise with the deconvolved PDF. Panel (a) displays
PDFs normalized in linear scale, while panel (b) shows PDFs rescaled to peak
at 1 in log scale, with velocity increments rescaled by the standard deviation.
The black curve represents the normal distribution. Error bars denote statistical
counting errors, given by

√
n. . . . . . . . . . . . . . . . . . . . . . . . . . . . 127

6.10 Kurtosis and Jarque-Bera (JB) Statistic of the deconvolved PDF as Functions
of Positional Separations. Red and blue denote the transverse and longitudinal
direction, respectively. The error bars indicate the standard error [2]. . . . . . . . 128

xi



6.11 (a) Log-Log plot of deconvolved SS
n (l), n = 1..7 fitted to the dashed lines anlεn .

The error bars are estimated from the
√
n statistical counting error. (b) intermittency

correction εn − n
3

versus n. The error bar is the two-sigma uncertainty of the
fitting in (a). The black dashed line represents the fit to the K62 theory, expressed
as −µn(n− 3), where the pink shaded area denotes the standard error of the fitting.130

6.12 Kinetic energy spectrum corresponding to the case of turbulence freely decaying
for 40ms after stirring. The blue dots represent the kinetic energy spectrum. The
green line is fitted to the inertial range and the red line is fitted to the crossover
range. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133

6.13 (a)Third order scalar structure function (blue) from direct calculation, and the
green line shows the linear fit in the inertial range. (b)Second order scalar structure
function (blue) from direct calculation, and the green curve shows the r2/3 fit in
the inertial range. The uncertainties are the two-sigma standard error of the mean
across the set of numerical runs. . . . . . . . . . . . . . . . . . . . . . . . . . . 134

6.14 (a) Log-Log plot of SS
n (l), n = 1..7 fitted to the dashed lines anlεn . The error

bar shows the two-sigma standard error of the mean. (b) intermittency correction
εn − n

3
versus n. The error bar is the two-sigma uncertainty of the fitting in (a). . 134

7.1 Concept. (a) A BEC is confined in a harmonic trap superimposed with a 1D
optical lattice (along ex, green), spatially modulating the condensate density
(red). The dashed and dotted lines call out a region of nominally constant mean
density and the left and right columns indicate the (b) state of the condensate and
(c) SF in the presence of a current. These were computed for a 5Er deep lattice
and plot: i. density (red), ii. current (green), iii. phase (orange), and iv. local
velocity (blue). The red dashed line plots the mean density ρ̄. . . . . . . . . . . 140

7.2 Moment of inertia in rotating systems computed using 2D GPE simulations. The
left column (a, c) indicates simulations in which the lattice is static while in
the right column (b, d) the lattice co-rotates with the confining potential. (a,
b) Angular momentum density for trap frequencies 2π×(56,36) and U0 = 10Er.
The colormap ranges from negative to positive, by normalizing to the largest
absolute angular momentum density. (c, d) Total momentum of inertia in traps
with frequencies 2π×(56,36) (top, green) and 2π×(36, 56) Hz (bottom, blue). In
(c) and (d), the cross markers are GPE simulated results of superfluid contribution
to the moment of inertia Isf/Ic. This is identified by calculating the gradient of
phase coarse-grained across a unit cell. The triangle makers are GPE simulated
results of the total moment of inertia I/Ic including the normal and superfluid
contributions. Dashed curves plot Isf/Ic and the solid curve plots I/Ic both
analytically derived from the superfluid hydrodynamics formalism. . . . . . . . 155

7.3 Modification of the phonon spectrum by the a = 266 nm optical lattice via BdG
calculation. U0 = 3Er. Dashed black and solid red curves mark excitations
created along ex and ey respectively. . . . . . . . . . . . . . . . . . . . . . . . . 160

xii



7.4 Bragg spectroscopy. Black and red symbols mark excitations created along ex
and ey respectively. (a) Transferred population fraction p as a function of frequency
difference δω with wavevetor δk/2π = 0.26 µm−1 and lattice depth U0 = 5.7Er.
The solid curve is a Lorentzian fit, giving the resonance frequency marked by the
vertical dashed line. (b) Phonon dispersion obtained from Bragg spectra. The
bold symbols resulted from (a) and the linear fit (with zero intercept) gives the
speed of sound. (c) Anisotropic speed of sound. The bold symbols are derived
from (b) and the solid curves are from BdG simulations (no free parameters [3]).
(d) SF density obtained from speed of sound measurements (blue markers, error
bars mark single-sigma statistical uncertainties). We compare with two models:
the red dashed curve plots a homogeneous gas BdG calculation, and the solid
black curve plots the result of Eq. (7.1). The simulations used our calibrated
experimental parameters. In (a)-(c) each point has uncertainty as shown in the
last point. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 162

7.5 Moment of inertia from scissors mode. (a-inset) Measured dipole mode frequencies
(circles) along with fits (curves) where the bare trap frequency is the only free
parameter for each curve. (a) Normalized scissors mode frequency. Blue and
green correspond toU0 = 0 trap frequencies (34, 51) Hz and (54, 36) Hz respectively.
(b) Moment of inertia in units of Ic. In (a) and (b) each point has uncertainty as
shown on the first point. Symbols are the data computed as described in the text,
and the solid curves are GPE predictions. . . . . . . . . . . . . . . . . . . . . . 165

xiii



List of Abbreviations

BEC Bose-Einstein condensate
VSF Velocity Structure Function
PIV Particle Image Velocimetry
VI Velocity Increments
PDF Probability Density Function
ARP Adiabatic Rapid Passage
AC Alternating Current
DC Direct Current
RF Radio Frequency
RWA Rotating Wave Approximation
OD Optical Depth
ToF Time of Flight
GPE Gross–Pitaevskii Equation
BdG Bogoliubov-de-Gennes
MOT Magneto-Optical Trap
PGC Polarization Gradient Cooling
AOM Acousto-optic Modulator
TTL Transistor-Transistor Logic
Re Reynolds Number
RANS Reynolds-averaged Navier-Stokes Equations
LHS Left hand side
RHS Right hand side
LIA Local Induction Approximation
LDV Laser Doppler Velocimetry
UHV Ultrahigh Vacuum
TEC Thermoelectric Cooler
UV Ultraviolet
CMOT Compressed Magneto-Optical Trap
SatAbs Saturated Absorption
FM Frequency modulation
DDS Direct Digital Synthesis
PLL Phase Lock Loop
UMD University of Maryland
MOSFET Metal-Oxide-Semiconductor Field Effect Transistor
NMOS N-Channel MOSFET
PMOS P-Channel MOSFET
DMD Digital Micromirror Device
LED Light-emitting Diode
EM Electromagnetic

xiv



PBS Polarizing Beam Splitter
PTAI Partial Transfer Absorption Imaging
NA Numerical Aperture
K41 Kolmogorov 1941
K62 Kolmogorov 1962
COM Center of Mass
ROI Region of Interest
QP Quadratic Programming
JB Jarque-Bera
SF Superfluid
UC Unit Cell

xv



Chapter 1: Introduction

Bose-Einstein condensates (BECs) represent one of the most intriguing states of matter

in quantum physics. First predicted by Satyendra Nath Bose and Albert Einstein in the early

twentieth century, BECs are formed when particles known as bosons are cooled to temperatures

near absolute zero. Under these extreme conditions, a significant fraction of the bosons occupy

the lowest quantum state, leading to the emergence of quantum phenomena on a macroscopic

scale. This unique state of matter provides an unparalleled platform for exploring a range of

quantum phenomena.

One of the most fascinating aspects of BECs is their behavior as quantum fluids. Unlike

classical fluids, where the flow is governed by the Newtonian mechanics, quantum fluids exhibit

collective behavior governed by quantum mechanics. The particles in a BEC act coherently,

displaying properties such as superfluidity, a flow with zero viscosity, and quantized vortices.

These characteristics make BECs an ideal system for studying fluid dynamics in a regime where

classical intuition that circulation of a vortex can vary continuously gives way to quantum behavior.

Despite the extensive research on BECs, the study of turbulence within these quantum fluids

remains relatively unexplored territory.

Turbulence is a fundamental phenomenon encountered in a wide range of fluids and at all

scales: from classical systems such as the earth’s oceans and atmosphere [4, 5]; confined and solar

1



plasmas [6, 7]; and the self-gravitating media of the large-scale universe [8] to quantum fluids

such as neutron stars [9], superfluid 4He [10] and atomic Bose-Einstein condensates (BECs) [11,

12]. All of these are characterized by complex patterns of fluid motion that span a wide range of

length scales. While the understanding of classical turbulence has matured in the past century [13],

that of quantum systems has many open questions [14]. For example, in BECs does there

exist a range of length scales—often termed the inertial scale—in which kinetic energy cascades

from large to small scale in accordance with a Kolmogorov scaling law? Although this scaling

was predicted only for incompressible fluids, it has been observed in virtually all turbulent

fluids [13]. Kolmogorov scaling is generally quantified in terms of velocity structure functions

(VSFs) that require knowledge of the fluid velocity field, which is difficult to measure in quantum

gas experiments. In this dissertation I present a velocimetry technique, analogous to particle

image velocimetry (PIV) [15, 16] employing spinor impurities as tracer particles, and thereby

obtain VSFs in turbulent atomic BECs in agreement with the Kolmogorov scaling.

Existing experimental evidence for turbulence in atomic BECs relies on time-of-flight

(TOF) measurements that have contributions from interaction-driven expansion [11] and the

momentum distribution [17]. This has no clear connection to VSFs, where the order-p VSF is

defined as Sp(l) = ⟨|δvl(x)|p⟩ as a function of displacement l which describes the typical change

in velocity

δvl(x) = [v(x+ l)− v(x)] · el (1.1)

along the direction of the displacement el 1. Without access to the VSF, the turbulence in

1The average ⟨· · · ⟩ is the ensemble (and ergodic for classical fluids) average over all positions x and displacement
directions el. Since the longitudinal and transverse VSFs are expected to be equal in isotropic systems such as ours,

2



atomic gases lacks a direct point of comparison with other fluids.

Unlike classical fluid flow, superfluid flow is strictly irrotational with a velocity field governed

by the phase of the superfluid order parameter ϕ via v = ℏ∇ϕ/m. Despite this, it is generally

believed that superfluid turbulence obeys the same scaling Sp(l) ∝ l(p/3) as classical fluids,

described by the initial K41 Kolmogorov theory [18, 19, 20]; in the case of 4He this has been ex-

perimentally verified [21, 22] for p ≤ 3. The more complete K62 theory [23] adds an intermittency

exponent that becomes important for large p and also predicts that the ensemble probability

density function of velocity increments (VI-PDF) is non-Gaussian, with “fat-tails.” Power-

law scaling behavior and the energy cascade have been observed in the momentum distribution

of homogeneously trapped BECs undergoing relaxation [12]; while this was interpreted in the

context of Kolmogorov-type scaling for order p = 2, the observed exponent departed from the

prediction of K41 theory and was instead interpreted using a wave turbulence model.

x position

y 
po

si
ti

on

Figure 1.1: (a) An example of velocity field. (b) Tracer particles (bright pink) are injected into
the system at t = 0. (c) At t = ∆t, tracer particles (bright pink) move from their initial position
(dark pink).

The cold-atom PIV technique, developed in our group, schematically illustrated in Fig. 1.1(a)-

we focus our exposition on the longitudinal VSF.

3



(c), enables directly measuring the velocity field and thereby both Sp(l) and the VI-PDF. In

this technique we prepare an initial velocity distribution [representative depiction in Fig. 1.1(a)],

then create localized “tracer particles” consisting of atoms in a different hyperfine state using

a spatially resolved technique [Fig. 1.1(b)], and, after a ∆t delay, measure the displacement of

the tracers [Fig. 1.1(c)]. This then leads directly to the local fluid velocity. The detail of this

technique is introduced in the Chap. 5.

1.1 Thesis overview

This thesis covers three broad topic areas: the first is background, including topics from

atomic physics that form the basis for our experimental techniques and turbulence theory that is

related to the science studied in the context. The second area is technical details, including some

details of our experimental apparatus, as well as how we use it to implement the cold atom PIV

technique. The final category is scientific results, including our measurement of VSFs from VI

which agrees well for low order of VSFs and intermittency is observed from high order of VSFs

and scale-dependent VI-PDF in a turbulent condensate.

Chapter 2 introduces important tools from atomic physics, including the light matter interactions

in which the tensor light shift is discussed in detail since it is directly related to the way we

generate the tracer particles in the system. I also discuss the laser cooling and evaporative cooling

techniques we apply to produce BEC in the lab.

Chapter 3 introduces the basic theory of classical and quantum turbulence, including the

Kolmogorov theory in classical incompressible turbulence, and classical compressible turbulence

using coarse graining, and finally the theoretical models and key experimental findings of quantum

4



turbulence. Velocimetry techniques are also discussed in this chapter. For readers who are only

interested in the atomic physics, this chapter can be skipped.

Chapter 4 describes the apparatus for producing 87Rb BEC, including the vacuum system,

magnetic field control, and optical control. Readers who are not interested in experimental details

can skip this chapter.

Chapter 5 presents the implementation of cold-atom PIV in velocity field measurement.

Benchmark examples include dipole mode and scissors mode in the harmonic trap and the quadrupular

irrotational flow field in a rotating harmonic trap.

Chapter 6 introduces the application of the cold-atom PIV in a turbulent condensate. The

velocity structure function is measured and agreed well with the Kolmogorov scaling. The

intermittency is also observed from the fat tail of the velocity increments probability distribution

function. Numerical simulation from a dissipation Gross Pitaevskii equation is also discussed.

Chapter 7 introduces a strange behavior in superfluid hydrodynamics, in which by imprinting

anisotropic spatial-dependent potential, the superfluid density becomes a tensor related to the

speed of sound. The theory and experimental details are discussed.

5



Chapter 2: Atom Light Interaction and Bose-Einstein Condensate

In this chapter, I will discuss the basic light atom interaction formalism used in the ultracold

gas community, from two-level atoms to the coherent transfers, in which the light shift is discussed

in detail. Then we discuss its application of laser cooling and evaporative cooling. The realization

in the lab is introduced in Chap. 4. Please note that the ultracold gas is considered as a single

particle in the first four sections of this chapter, which means that the scattering between atoms

is neglected. In the last section, we will add the interaction when we discuss the Bose-Einstein

condensate.

2.1 Two-level systems

A two-level atom is the simplest configuration for quantum mechanics, where the Hamiltonian

can be expressed by

Hatom = ℏω|e⟩⟨e|. (2.1)

Here, the energy of the ground state |g⟩ is set to 0 for simplicity, and the energy difference

between the excited state |e⟩ and the ground state |g⟩ is ℏω. Now we introduce the light to the

system. In a semi-classical picture, the light can be expressed by its electrical field.

E(r, t) = ϵ(r, t) exp (iνt− ik · r) + C.C, (2.2)

6



where ϵ is a slowly varing amplitude compared to the spatial frequency k, and C.C is the complex

conjugate. In atomic physics, under dipole approximation, we write the atom-light interaction as

a perturbative term −d · E(r = r0, t), where d = qer is the dipole moment of the atom with

charge qe, and r0 is the position of the atom. A rigorous derivation can be obtained from the

Hamiltonian of an electron with charge qe and mass m moving in an electromagnetic (EM) field.

H =
1

2m
[pe − qeA(r, t)]2 + qeϕ(r, t) + V (r, t)− µ ·B(r, t), (2.3)

where A and ϕ are the vector and scalar potential of the EM field, B is the magnetic field, µ is

the magnetic dipole moment, and V is the core Coulomb potential. Using the Coulomb gauge,

that is, ∇·A = 0; ϕ = 0, the vector potential satisfies the wave equation and can be decomposed

as A(r, t) = A(t) exp (ik · r). Expand the vector potential around the position of the atom r0,

A = A(t)eik·r0 [1 + ik · δr]. (2.4)

Under the dipole approximation k · δr ≪ 1, the vector potential becomes spatially independent,

that is, A ≈ A(t)eik·r0 . Since the position dependence of A is gone, thereby A commutes to pe.

So we basically treat the EM field as a classical field and the Hamiltonian becomes the following.

H =
1

2m
[pe

2 − 2qepe ·A+ q2eA
2] + V − µ ·B. (2.5)

The vector potential coupled to the momentum can be transferred away by a unitary transformation

|ψ′⟩ = R|ψ⟩ = exp [−iqe
r ·A(t)

ℏ
]|ψ⟩, (2.6)

7



and the Hamiltonian becomes

H ′ = RHR† + iℏ
∂R

∂t
R† =

pe
2

2m
+ iℏ(

−iqer
ℏ

· ∂A
∂t

) + V (r)− µ ·B, (2.7)

where RpeR
† = pe + qeA is used. The second term on the right-hand side (RHS) is exactly the

interaction term Hint = −d · E(r0, t). So, we can decompose the Hamiltonian into two terms.

Hatom =
pe

2

2m
+ V (r), Hint = −d · E(r0, t)− µ ·B(r0, t) (2.8)

The atom-light interaction consists of the electrical dipole interaction d · E(r0, t) and the

magnetic dipole interaction −µ · B(r0, t). Electric dipole transitions only have a non-vanishing

matrix element between quantum states with different parity. Magnetic dipole transitions in

contrast couple states with the same parity. The response of the magnetic dipole transition is

much weaker than that of electric dipole transitions, so normally we consider only the electric

dipole interaction unless it is forbidden by selection rules.

2.2 Coherent transfer between two levels

Rabi oscillations [24] and adiabatic rapid passage (ARP) [25] are two frequently used

methods for transferring atoms between different states. In this section, the rotating wave approximation

(RWA) [26] is applied to obtain an effective Hamiltonian without a fast oscillation in the optical

frequency, from which we discuss Rabi flopping and ARP.

8



2.2.1 Rotating Wave Approximation

For a two-level atom, the electric dipole moment d can be decomposed under the eigenbasis

of Hatom

d = dge|g⟩⟨e|+ deg|e⟩⟨g|, (2.9)

where for spherical symmetric atoms dgg = dee = 0. Under the dipole approximation E(t) =

ϵ(t)e−iνt + C.C. So, the Hamiltonian becomes

H = ℏω|e⟩⟨e| − (dge|g⟩⟨e|+ deg|e⟩⟨g|) · (ϵ(t)e−iνt + ϵ∗(t)eiνt). (2.10)

To obtain an effective Hamiltonian without the fast oscillating frequency ν, we go to the rotating

frame by applying U = |g⟩⟨g|+ e−iνt|e⟩⟨e|. The effective Hamiltonian then becomes

Heff = U †HU − iℏ(
d

dt
U)U † = ℏ(ω − ν)|e⟩⟨e|+ deg · ϵ(t)|e⟩⟨g|+ dge · ϵ∗(t)|g⟩⟨e|, (2.11)

where the fast oscillation with frequency 2ν is ignored under RWA.

2.2.2 Rabi Oscillation and Adiabatic Rapid Passage

When the excitation is on resonance (ν = ω), and the atom is initially on the ground

state, the atom population will oscillate between the ground and excites states by cos2 Ωt
2

, where

Ω = 2deg·ϵ
ℏ is the Rabi frequency. Introducing the detuning δ = ω− ν, the effective Hamiltonian

9



is

Heff = ℏ

 0 −Ω∗
2

−Ω
2

δ

 . (2.12)

with eigenenergies ϵ± = ℏδ
2
∓ ℏ

2

√
|Ω|2 + δ2 and eigenstates which often called dressed states

|ϕ±⟩ =
1

N±

(
Ω∗|g⟩+ (−δ ±

√
|Ω|2 + δ2)|e⟩

)
. (2.13)

Note that |g⟩ and |e⟩ are already rotated; for convenience, I use the same label. The normalized

factors are N2
± = |Ω|2 +

(
−δ ±

√
|Ω|2 + δ2

)
.

When the excitation is far off-resonant (|δ| ≫ ∥Ω|), the dressed state reduces to the bare

state without the coupling, i.e., for the blue-detuned case

|ϕ+⟩ ≈
Ω∗

2δ
|g⟩+ |e⟩ ≈ |e⟩, |ϕ−⟩ ≈ |g⟩ − Ω

2δ
|e⟩ ≈ |g⟩. (2.14)

ϵ+ ≈ ℏ(δ +
|Ω|2
4δ

), ϵ− ≈ −ℏ
|Ω|2
4δ

. (2.15)

We see that the eigenenergy of the ground and excited states shifts by ℏ |Ω|2
4δ

in opposite directions;

this effect is called the AC stark shift. If Ω is inhomogeneous, it can be used to create a trapping

potential for optical tweezers and dipole trap.

Consider the case that we turn on the coupling when the atoms are initially in the ground

state; then the detuning is swept adiabatically from far blue detuned to far red detuned, where the

Landau-Zener transition is negligible if the swept is slow compared to the Rabi frequency, and

we transfer the state from |g⟩ along the channel |ϕ−⟩ to the case δ ≈ +∞ where |ϕ−⟩ reduces to

|e⟩. Finally, if we turn off the coupling, the atoms will end up in the bare state |e⟩. The above

10



process is called an adiabatic rapid passage.

2.3 Optical Stark shift

In the last section, we have seen the effect of the AC stark shift when the excitation is far

from resonant. In principle this light shift can be derived from the second-order time-dependent

perturbation theory without RWA, i.e.,

∆Eα = −
∑
β ̸=α

2ωβα|⟨α|d · ϵ̂|β⟩|2|ϵ|2
ℏ(ω2

βα − ω2)
, (2.16)

where ωβα = (Eβ − Eα)/ℏ. For the large detuning and two-level case, this reduces to the AC

Stark shift in the previous section, and the ωβα on the numerator shows the opposite shift for the

ground and excited states. We can obtain the polarizability α(ω) of the system by −∆E/|ϵ|2, and

the result

α(ω) =
∑
β ̸=α

2ωβα|⟨α|d · ϵ̂|β⟩|2
ℏ(ω2

βα − ω2)
(2.17)

is Kramers-Heisenberg formula for the polarizability.

However, since the electric field is a vector, the polarizability should be a tensor which

satisfies dµ(ω) = αµν(ω)ϵν , so the light shift in general should be written in the form ∆E =

−αµν(ω)ϵµϵν with

αµν(ω) =
∑
β ̸=α

2ωβα⟨α|dµ|β⟩⟨β|dν |α⟩
ℏ(ω2

βα − ω2)
. (2.18)

Note that the above expression assumes linear polarization, if the light is circularly polarized the

ωβα in the numerator should be replaced by ω.

Since the state in the atomic BEC is typically described under the hyperfine basis |F,mF ⟩,

11



we will rewrite the polarizability tensor for the state |F,mF ⟩

αµν(ω) =
∑
F ′,m′

F

2ωF ′F ⟨F,mF |dµ|F ′,m′
F ⟩⟨F ′,m′

F |dν |F,mF ⟩
ℏ(ω2

F ′F − ω2)
≡

∑
F ′

2ωF ′FTµν
ℏ(ω2

F ′F − ω2)
, (2.19)

where the angle dependence is merely on the dipole product tensor

Tµν =
∑
m′

F

⟨F,mF |dµ|F ′,m′
F ⟩⟨F ′,m′

F |dν |F,mF ⟩. (2.20)

It is straightforward to decompose the rank-2 tensor into scalar T (0), vector T (1), and tensor parts

T (2) using the irreducible tensor formalism, that is,

Tµν =
1

3
T (0)δµν +

1

4
T (1)
σ ϵσµν + T (2)

µν . (2.21)

T (0) is the trace of Tµν , that is,

T (0) = Tµµ (2.22)

T (1) is the anti-symmetric part of the tensor, i.e.,

T (1)
σ = (Tµν − Tνµ)ϵσµν . (2.23)

T (2) is the rest part.

We now try to express the scalar and vector part by the inner product and cross product of

12



two vectors A and B, i.e.,

T
(0)
0 = − 1√

3
A ·B, T (1)

q =
i√
2
A×B, (2.24)

where the prefactors are given from the ⟨1, q; 1,−q|0, 0⟩ = −(−1)q/
√
3 and |⟨1, q′; 1, q−q′|1, q⟩| =

1/
√
2 unless q = q′ = 0. So T (0) becomes

T (0) =
∑
m′

F

⟨F,mF |dµ|F ′,m′
F ⟩⟨F ′,m′

F |dµ|F,mF ⟩

= −
√
3⟨F,mF |

∑
m′

F

d|F ′,m′
F ⟩⟨F ′,m′

F |d

(0)

|F,mF ⟩

= −
√
3⟨F ||

∑
m′

F

d|F ′,m′
F ⟩⟨F ′,m′

F |d

(0)

||F ⟩ ⟨F,mF |F,mF ; 0, 0⟩

= −
√
3(−1)2F

√
2F ′ + 1


1 1 0

F F F ′

 ⟨F ||d||F ′⟩⟨F ′||d||F ⟩

= −
√
3(−1)F+F ′√

2F + 1


1 1 0

F F F ′

 |⟨F ||d||F ′⟩|2

= |⟨F ||d||F ′⟩|2.

(2.25)

We used the Wigner-Eckert theorem in the third row, and further express the reduced matrix

13



elements by Wigner-6j symbol using

⟨F ||

∑
m′′

F

d|F ′′,m′′
F ⟩⟨F ′′,m′′

F |d

(k)

||F ′⟩

=(−1)k+F+F ′√
(2F ′′ + 1)(2k + 1)


1 1 k

F ′ F F ′′

 ⟨F ||d||F ′′⟩⟨F ′′||d||F ⟩.

(2.26)

The proof can be found in chapter 7 of [27], where the 1 in the Wigner-6j indicates the rank of the

dipole operator, and the Wigner-6j value is (−1)−F−F ′−1/
√
3(2F + 1). In the second last row of

Eq. (2.25) we also used the conjugate of the reduced matrix elements

⟨F ′||T (k)||F ⟩ = (−1)F
′−F

√
2F + 1

2F ′ + 1
⟨F ||T (k)||F ′⟩∗. (2.27)

Similarly, we can express the vector part

T (1)
q = 2

∑
m′

F

⟨F,mF |d|F ′,m′
F ⟩ × ⟨F ′,m′

F |d|F,mF ⟩

= −i2
√
2⟨F,mF |

∑
m′

F

d|F ′,m′
F ⟩ × ⟨F ′,m′

F |d

(1)

q

|F,mF ⟩

= −i2
√
2⟨F ||

∑
m′

F

d|F ′,m′
F ⟩ × ⟨F ′,m′

F |d

(1)

||F ⟩ ⟨F,mF |F,mF ; 1, q⟩

= −i2
√
2(−1)2F+1

√
3(2F ′ + 1)


1 1 1

F F F ′

 ⟨F ||d||F ′⟩⟨F ′||d||F ⟩ mF δq0
F (F + 1)

= i2
√
2(−1)F+F ′√

3(2F + 1)


1 1 1

F F F ′

 |⟨F ||d||F ′⟩|2 mF δq0
F (F + 1)

.

(2.28)

14



Note that the vector part is linear to mF , and the only non-vanishing part is T (1)
0 , thereby the

electric field vector should be ∼ (ϵ× ϵ∗)(1)0 to contract the polarizabiity.

Finally, following the similar procedure, the tensor part can be expressed as

T (2)
q = (−1)F+F ′

√
5(2F + 1)

F (F + 1)(2F − 1)(2F + 3)


1 1 2

F F F ′

 |⟨F ||d||F ′⟩|2[m2
F−F (F+1)]δq0

(2.29)

Notice that all non-vanishing irreducible components T (k)
q should satisfy q = 0 to preserve

mF , and thereby contract the T (k)
q to a rank-0 energy shift, the corresponding electric field tensor

(EE)
(k′)
q′ should satisfy k′ = k and q′ = 0. So, the total energy shift can be expressed by

∆E(F,mF , ω)

=−
∑
F ′

2

ℏ(ω2
F ′F − ω2)

[
ωF ′F

3
T (0)|ϵ|2 + ω

4
T

(1)
0 (ϵ∗ × ϵ)z +

ωF ′F√
6
T

(2)
0 (3|ϵz|2 − |ϵ|2)

]
.

(2.30)

We note that the first term in the square bracket is independent of the polarization of the electric

field, so it is called the scalar light shift, and we will show below that it reduces to the familiar

AC Stark shift for large detunings. The second term is called the vector light shift, since the rank

of the tensor is 1, where ω in the second term indicates that the vector light shift comes from the

circular polarized light. The last term is related to the rank-2 tensor, so we call it tensor light

shift.

Writing out the mF dependence, we have

δE(F,mF , ω) =− α(0)(F, ω)|ϵ|2 − α(1)(F, ω)(iϵ∗ × ϵ)z
mF

F

− α(2)(F, ω)
3|ϵz|2 − |ϵ|2

2

3m2
F − F (F + 1)

F (2F − 1)
,

(2.31)

15



where the scalar, vector and tensor polarizabilities are

α(0)(F, ω) =
∑
F ′

2ωF ′F |⟨F ||d||F ′⟩|2
3ℏ(ω2

F ′F − ω2)
,

α(1)(F, ω) =
∑
F ′

(−1)F
′+F+1

√
6F (2F + 1)

F + 1


1 1 1

F F F ′


ω|⟨F ||d||F ′⟩|2
ℏ(ω2

F ′F − ω2)
,

α(2)(F, ω) =
∑
F ′

(−1)F
′+F

√
40F (2F + 1)(2F − 1)

3(F + 1)(2F + 3)


1 1 2

F F F ′


ωF ′F |⟨F ||d||F ′⟩|2
ℏ(ω2

F ′F − ω2)
.

(2.32)

It is worth noting that the vector light shift is linear to mF , thereby it is similar to the

weak field Zeeman shift, and we can view the AC electric field as an effective magnetic field, i.e.

Beff ∼ (iϵ∗ × ϵ)z. Clearly, only a circularly polarized light contributes to the vector light shift;

for the case ϵ̂ = (êx + iêy)/
√
2, the effective B field is along (iϵ̂∗ × ϵ̂)z = −êz.

2.3.1 Tune-out (magic) Wavelength

In this section, we will focus on the case of large detuning, where the detuning is much

greater than the hyperfine splitting (∼GHz) so we can ignore the hyperfine structure. In this

case the light shift can be calculated from the two-level AC Stark shift −ℏ|Ω|2/4δ, but if more

energy levels are involved, it is possible to cancel out the scalar light shift, and the wavelength

of the incident light is referred to as the tune-out (magic) wavelength. To see this, we first

decompose the hyperfine structure reduced matrix elements into fine structure reduced matrix

16



elements. Using

⟨F ||d||F ′⟩ ≡ ⟨J, I;F ||d||J ′, I ′;F ′⟩

= ⟨J ||d||J ′⟩(−1)F
′+J+I+1

√
(2F ′ + 1)(2J + 1)


J J ′ 1

F ′ F I

 ,

(2.33)

and Biedenharn—Elliott sum rule [27], the polarizabilities become

α(0)(F, ω) ≈
∑
J ′

2ωJ ′J |⟨J ||d||J ′⟩|2
3ℏ(ω2

J ′J − ω2)
,

α(1)(F, ω) ≈
∑
J ′

(−1)−2J−J ′−F−I+1

√
6F (2F + 1)

F + 1
(2J + 1)

ω|⟨J ||d||J ′⟩|2
ℏ(ω2

J ′J − ω2)


1 1 1

J J J ′



J J 1

F F I

 ,

α(2)(F, ω) ≈
∑
J ′

(−1)−2J−J ′−F−I

√
40F (2F + 1)(2F − 1)

3(F + 1)(2F + 3)
(2J + 1)

ωJ ′J |⟨J ||d||J ′⟩|2
ℏ(ω2

J ′J − ω2)


1 1 2

J J J ′



J J 2

F F I

 ,

(2.34)

where the ωF ′F is replaced by ωJ ′J due to large detuning.

Now we apply this to the atom we are interested, i.e., 87Rb. The energy level of the ground

state and D line is shown in Fig. 2.1. The dashed level corresponds to the tune-out wavelength

which makes the α(0)(F, ω) vanishing.

From large detuning we have 2ωJ ′J/(ω
2
J ′J − ω2) ≈ 1/δ, thereby

α(0)(F, ω) ≈ |⟨J = 1
2
||d||J ′ = 1

2
⟩|2

3ℏδ1
+

|⟨J = 1
2
||d||J ′ = 3

2
⟩|2

3ℏδ2
, (2.35)

17



52S1/2

52P3/2

52P1/2

795
nm

780
nm

δ2

δ1

Figure 2.1: 87Rb D line.

where ⟨J = 1
2
||d||J ′ = 1

2
⟩ = 2.992ea0, ⟨J = 1

2
||d||J ′ = 3

2
⟩ = 4.227ea0 [28]. Here, a0 is the

Bohr radius. Note that the ratio of the reduced matrix element is approximately
√
2. Therefore,

at the tune-out wavelength λmagic, the δ2 + 2δ1 = 0, that is,(1/780 − 1/λmagic) + 2(1/795 −

1/λmagic) = 0. So, the tune-out wavelength is λmagic ≈ 790 nm.

Normally we don’t need to consider the effect of the vector light shift and tensor light shift

because the scalar light shift dominates, however, around the tune-out wavelength the higher rank

light shift is not negligible.

It is interesting to see that the tensor polarizability α(2)(F, ω) is also vanishing around

λmagic for the ground state of the 87Rb, owing to the vanishing Wigner-6j, i.e.,


1 1 2

1
2

1
2

1
2

 =


1 1 2

1
2

1
2

3
2

 = 0. So at the tune-out wavelength, if the incident light is not purely linear

polarized, the only nonvanishing light shift comes from the vector light shift, and its effect is

analogous to an effective magnetic field.

18



2.3.2 AC vector light shift as RF

In the previous section, we have seen that the vector light shift can be analogous to a DC

effective magnetic field. Can we make it an AC effective magnetic field so that it can drive the

magnetic dipole transition? The time dependence can be introduced from the two-tone circularly

polarized light E = ϵ[(êx + iêy)/
√
2](e−i(ω+∆ω)t + e−i(ω−∆ω)t) + C.C., where ω can be selected

to tune out the scalar light shift and ∆ω is the modulation frequency that can be easily added

in the experiments by an acousto-optic modulator (AOM). Consequently we have an effective B

field Beff ∼ −2[cos (2∆ωt) + 1]êz, thereby apart from the DC effective magnetic field along

the light propagation, we also obtain an effective AC magnetic field with angular frequency 2∆ω

along the same direction. We can view it as an effective RF field, which is capable of driving the

magnetic dipole transition between different mF states in the 52S1/2 manifold of 87Rb, where the

electric dipole transition is forbidden due to J = J ′ = 0.

2.4 Laser Cooling

Previous sections neglect spontaneous emission and dissipation in the system. In this

section, we add the dissipation into the system and discuss the laser cooling techniques from

the optical Bloch equation.

19



2.4.1 Optical Molasses

With the spontaneous emission, the pure state decoheres to a mixed state typically described

by a density matrix ρ satisfying

∂tρee = i
Ω

2
(ρeg − ρge)− Γρee,

∂tρgg = −iΩ
2
(ρeg − ρge) + Γρee,

∂tρge = −(
Γ

2
− iδ)ρge − i

Ω

2
(ρee − ρgg),

∂tρeg = −(
Γ

2
+ iδ)ρeg + i

Ω

2
(ρee − ρgg),

(2.36)

where Ω and δ are the Rabi frequency and the detuning, and Γ is the excited state decay rate [29].

The force on an atom can be derived using Heisenberg equation of motion, i.e.,

F =
d

dt
⟨p⟩ = i

ℏ
⟨[H,p]⟩ = −⟨∇H⟩ = −Tr(ρ∇H). (2.37)

Since we are interested in the force induced by the laser-atom interaction, and under dipole

approximation the force becomes

F = −ℏ
2
(ρ∗eg∇Ω + ρeg∇Ω∗). (2.38)

We know that one part of the force is related to the scalar light shift, which only depends on the

magnitude of Rabi frequency; thus we decompose it to be a phase irrelevant and a phase relevant

20



term using Ω = |Ω|eiϕ, i.e.,

∇Ω = Ω

(
∇|Ω|
|Ω| + i∇ϕ

)
⇒ ∇(lnΩ) = ∇(ln |Ω|) + i∇ϕ. (2.39)

We then use the steady-state solution of Eq. (2.36)

ρeg =
−iΩ

2(Γ
2
+ iδ)(1 + s)

, s =
|Ω|2

2[(Γ
2
)2 + δ2]

(2.40)

to replace the ρeg in Eq. (2.38), where s is the saturation parameter. We then obtain the following.

F =
ℏs

1 + s
(δ∇ ln |Ω|+ Γ

2
∇ϕ). (2.41)

The first term, phase irrelevant, is the dipole force used to trap atoms; the second term, phase

relevant, is the radiation-pressure force.

The previous argument assumes the atom is at rest; now we add velocity into the expression.

Suppose that the atom is moving with velocity v, the Rabi frequency Ω now becomes time

dependent, that is,

∂Ω

∂t
= v · ∇Ω. (2.42)

The steady-state solution of the optical Bloch equations has to change due to this time dependence.

Treating the change as a perturbation and keeping the terms to linear order of v, we have

∂ρeg
∂t

= ρeg

(
1

Ω

∂Ω

∂t
− 1

1 + s

∂s

∂t

)
,

∂(ρee − ρgg)

∂t
= −(ρee − ρgg)

1

1 + s

∂s

∂t
. (2.43)

21



The final steady-state solution becomes

ρeg =
−iΩΓγ∗v/2

Γv|γv|2 + |Ω|2Re[γv]
, (2.44)

where Re stands for the real part and

Γv = Γ− s

1 + s
v ·

(
∇s

s

)
, γv =

Γ

2
+ iδ +

1− s

1 + s
v · ∇s

2s
+ iv ·∇ϕ. (2.45)

We can obtain the total force on the moving atom by putting this in Eq. (2.38). Now we study

the case that the incident light is a plane wave with wavenumber k, so the ∇s = ∇|Ω| = 0 and

∇ϕ = k. We then obtain

ρeg =
−iΩ

2[Γ/2 + iδ(v)](1 + s(v))
, s(v) =

|Ω|2
2[(Γ

2
)2 + δ(v)2]

(2.46)

where δ(v) = δ + k · v. We notice that the only difference between Eq. (2.46) and Eq. (2.40) is

a Doppler shift added to the detuning. And the radiation-pressure force becomes

Frad =
ℏkΓ
2

s(v)

1 + s(v)
. (2.47)

For a red-detuned light (δ > 0), the force is frictional, as can be seen from Fig. 2.2. And if two red

detuned beams with wave vector k and −k are involved, the deceleration effect is symmetric and

can be used to slow the motion of the atomic gas along the direction of k. Furthermore, if three

orthogonal pairs of these Doppler cooling beams are added to the system, the atom’s deceleration

will be along every direction, and the temperature of the gas can be largely decreased. This

22



technique is typically referred to as optical molasses [30]. The limit of this cooling is T =

−6 −4 −2 0 2 4 6

−0.4

−0.2

0

0.2

0.4

v[Γ/k]

F
r
a
d
[ℏ
k
Γ
]

1D optical molasses

Figure 2.2: Frictional force of Doppler cooling when s(0) = 2 and δ = Γ. The blue (purple)
dashed line is when k and v have the same (opposite) direction, and the red solid line is the sum
of the two.

ℏΓ/2kB, which is called the Doppler limit.

2.4.2 Magneto-optical Trap

The optical molasses can cool the temperature of the gas, but is unable to trap the atoms.

In this section, I will introduce a technique that trapping and cooling can be simultaneously

achieved by using circular polarized red-detuned light with anti-Helmholtz coil pairs, and the

trap is referred to as the magneto-optical trap (MOT) [31].

For an anti-Helmoholtz pair with the symmetry axis along z, around the center of the trap,

the magnetic field produced can be expressed by B(x, y, z) = Bxxx + Byyy + Bzzz, and due to

23



symmetry Bxx = Byy and Bzz = −2Bxx derived from ∇ ·B = 0.

Now let us consider the cooling along the z direction. The magnetic gradient along z gives

a spatially dependent Zeeman shift as shown in Fig. 2.3, where for simplicity the ground and

excited states are |F = 0⟩ and |F = 1⟩, so we only need to consider the Zeeman shift on the

excited states. Meanwhile, a pair of red-detuned circular-polarized counter-propogating beams

propogates along the z direction. The polarizations of the two beams are opposite compared to

the quantization axis êz.1

B B

νν

ν0
σ+ σ−

F = 1

F = 0

mF = 0

mF = 0

mF = 1

mF = −1

z

Figure 2.3: MOT scheme.

We now calculate the radiation force for an atom located at z with velocity v along êz. The

detuning for the σ− and σ+ light can be expressed by

δ(−)(z, v) = −ν + (ν0 −∆B(z))− kv, δ(+)(z, v) = −ν + (ν0 +∆B(z)) + kv, (2.48)

where ℏ∆B(z) = gµBBzzz is the Zeeman splitting at z. We can obtain the radiation-pressure

force by plugging this into Eq. (2.47), and the only difference compared to the optical molasses is
1But the polarization is the same in the optical reference frame, i.e., if we ride on the beam the polarization would

be the same.

24



the kv should be replaced by kv+∆B(z), therefore the radiation-pressure force can be expressed

by

Frad(v, z) = −α(v +∆B(z)/k), (2.49)

where α is the damping coefficient. In addition to the cooling force −αv, the velocity-independent

term −α∆B(z)/k is the trapping force that pushes the atoms back to z = 0. So if we have

this setup in 3 orthogonal directions, we can obtain the 3D optical molasses and a 3D trap

simultaneously, and it is called 3D MOT.

2.4.3 Polarization gradient cooling

The MOT and optical molasses cannot beat the Doppler limit because of random heating

from the spontaneous emission. In this section, I will introduce a sub-Doppler cooling scheme

which exploits the internal structure of atoms and the fact that light can exert forces on atoms

beyond the simple Doppler effect. The force comes from the spatial dependence of the polarization;

thus, it is also called polarization gradient cooling (PGC) [32]. Generally speaking, the PGC

requires two red-detuned counterpropogating beams with orthogonal polarization. In this section,

I will discuss two cases: (i) the orthogonal linear polarized configuration and (ii) the orthogonal

circular polarized configuration.

2.4.3.1 Linear polarization configuration

In the section. 2.3, we have seen that the AC Stark shift can be decomposed to the scalar,

vector, and tensor parts. Here, we will use the vector light shift that lifts the degeneracy of

ground and excited states of the atom analogous to the Zeeman effect under a spatial-dependent

25



polarization setup.

Consider the two counter-propogating beams are polarized along êx and êy, i.e.,

E(z, t) = E0e
ikz−iωtêx + E0e

−ikz−iωtêy + C.C.

= E0e
−iωt[eikz êx + e−ikz êy] + C.C.

(2.50)

Note that the |E| is spatial independent, but the polarization is spatial independent, thereby

the scalar light shift cannot break the degeneracy and the vector light shift will lift the degeneracy

if the effective magnetic field is non-vanishing.

Using Eq. (2.31) the effective magnetic field becomes

Beff (z) ∼ (iϵ∗ × ϵ)z êz = −2E2
0 sin(2kz)êz, (2.51)

where ϵ = E0[e
ikz êx + e−ikz êy].

For simplicity, we assume that the ground state is |F = 1⟩, and the spatial-dependent shift

can be seen in Fig. 2.4.

Here, we will explain the cooling process qualitatively. Due to the optical pumping effect,

the steady-state solution at z = (1
8
+ n

2
)(2π/k) would be |1,−1⟩ where the light is σ−. When

the |1,−1⟩ atom moves by 1
4
(2π/k), it reaches the valley of |1, 1⟩, where the light is σ+, then

again due to optical pumping, the atom will be transferred to |1, 1⟩ in an anti-stokes process

which reduces the kinetic energy of the system. This argument is similar to the Greek mythology

of Sisyphus who rolling a stone upward to the peak of a mountain but it rolls off to the valley

”spontaneously”, and this process repeats over and over, so we typically refer to it as Sisyphus

cooling. A more rigorous derivation using the optical Bloch equation is given in [33].

26



−1
2
−3

8
−1

4
−1

8
0 1

8
1
4

3
8

1
2

σ−σ− Lin LinLin σ+ σ+

z[2π/k]

E
ne

rg
y

1D Sisyphus Cooling

|1,−1⟩
|1, 1⟩
|1, 0⟩

Excited state

Figure 2.4: Sisyphus cooling scheme.

2.4.3.2 Circular polarzation configuration

Consider the light is incident the same as the MOT scheme but without any external

magnetic field, then the electric field can be written as

E(z, t) = E0e
ikz−iωt(êx + iêy) + E0e

−ikz−iωt(êx − iêy) + C.C.

= 2E0e
−iωt[cos(kz)êx − sin(kz)êy] + C.C.

(2.52)

The |E| is spatial-independent, thereby the scalar light shift does not lift the degeneracy.

Note that the vector light shift is also 0 since ϵ = ϵ∗, so the AC Stark shift does not break the

degeneracy, thereby the cooling doesn’t come from the Sisyphus effect.

The polarization of the electric field at each point is linear. Consider an atom moving

27



along êz with velocity v, then if we go to the frame rotating with −kv about êz, the polarization

is fixed. If the atom’s initial position is z = 0, then in the rotating frame the electric field

is constantly along êx. In the rotating frame, an additional gauge transformation term would

modify the Hamiltonian by ∆H = kvFz, which breaks the degeneracy by a purely motion-

induced effect. Again consider the |F = 1⟩ atoms, for the steady state the |1,−1⟩ atoms tend to

move along êz, while the |1, 1⟩ tends to move along −êz. If the excited state has cyclic transitions

with the |F = 1⟩ manifold, then the atom tends to absorb more photons from the cyclic transition

because of the greater Clebsch-Gordon coefficient. In this case, the |1,−1⟩ (|1, 1⟩) tends to absorb

more photons from σ−(σ+) light which propagates along −êz (êz). Thus, an atom moving right

(left) tends to absorb photons moving left (right), which results in a cooling effect. A rigourous

derivation is given in [33].

Note that the above argument relies on a zero external magnetic field. A nonzero magnetic

field results in a nonzero shift vB in the radiation-pressure force Frad = −α(v − vB), where α is

the damping coefficient. Under a zero external magnetic field, the temperature limit of the PGC

is the recoil limit Tr, i.e., kBTr = ℏk2
m

, where m is the mass of the atom and k is the wavenumber

of the light.

2.4.4 Evaporative Cooling

The laser cooling technique has a temperature limit close to the recoil limit, but still not

cold enough to reach the Bose-Einstein condensation. To lower the temperature, one has to

sacrifice the number of atoms using evaporation which is similar to the evaporation of hot water,

where the steam with high kinetic energy goes away and the remaining water molecules with

28



less kinetic energy stay in the cup. This process does not conserve the number of particles, so to

reach a high phase-space density one cannot evaporate too aggressively. Atomic physicists apply

this evaporative cooling technique to the cold atoms [34]. In this section, I will introduce two

commonly used evaporative cooling techniques: (i) RF evaporation in the magnetic trap and (ii)

dipole evaporation in optical dipole trap.

2.4.4.1 RF Evaporation

Consider atoms loaded in the quadrupole magnetic trap, where the atoms with high energy

are further away from the trap center. Using 87Rb as an example, in the |F = 1⟩ manifold only

the state |1,−1⟩ is magnetic trappable. So, evaporation can be realized by pumping the atoms

with high kinetic energy into magnetic untrappable states, that is, |1, 0⟩ and |1, 1⟩. These states

can be coupled by the RF field, and we can understand it using the dressed state picture. The

coupled RF transition to the untrapped state is prominent when the detuning is within the Rabi

frequency, that is, |δ| ∼ |Ω|.

Again we only study the 1D case, the energy level of a magnetic trap along (0, 0, z) is

H = gmzµBBzzz, so the energy splitting is δω = gµBBzzz/ℏ. Then given the RF frequency

ωRF , some of the atoms in the spatial range ℏ
gµBBzz

(ωRF ± Ω) will become untrappable. Then

atoms with less energy that are located closer to the trap center will collide with the remaining

trappable atoms and rethermalize to a lower temperature. This RF effect is sometimes called the

”RF knife”. If we continuously (adiabatically) reduce the ωRF , we placed the RF knife close to

the center of the trap, and the temperature will be much lower at the cost of fewer atoms.

29



2.4.4.2 Dipole Evaporation

Apart from the magnetic trap, the atoms can be loaded into optical dipole trap which

typically uses the dipole force (scalar light shift) from a red-detuned Gaussian beam. Around

the center of the trap, it can be approximated to a harmonic trap. By reducing the light intensity,

the trap depth decreases, so that the atoms away from the trap center with high kinetic energy

go out of the trap and become untrappable. If this trap depth reduction is done adiabatically, the

atoms can collide and re-thermalize to a lower temperature.

2.5 Atomic Bose-Einstein Condensates

In the context of dilute Bose gases, the interactions between particles, though weak due to

the low density of the gas, are crucial for understanding the system’s behavior. These interactions

are commonly approximated using a pseudo-potential that is characterized by the s-wave scattering

length. This approximation is particularly useful because it simplifies the complex nature of

interatomic forces into a manageable form, especially at low temperatures where s-wave scattering

predominates. In this section, I will introduce the Gross-Pitevskii equation (GPE) [35, 36] to

describe the ground state of BEC, and the study the elementary excitation using Bogoliubov-de-

Gennes (BdG) formalism [37, 38].

2.5.1 Gross-Pitevskii equation

In a completely condensed system, each particle is in an identical single-particle state,

denotedψ(r), fulfilling the normalized condition
∫
d3r|ψ(r)|2 = 1. The many-particle condensate’s

wavefunction in a mean-field approach is represented as the symmetrized product of these single-

30



particle wavefunctions due to its bosonic nature, expressed as Ψ(r1, r2, ..., rN) =
∏N

i=1 ψ(ri).

This wavefunction is influenced by three key components: kinetic energy, potential energy,

and interaction energy. The mean-field interaction energy g is dependent on the s-wave scattering

length a in a weakly interaction Bose gas, i.e. g = 4πℏ2a
m

. Consequently, the system’s Hamiltonian

is formulated as:

H =
N∑
i=1

(
p2i
2m

+ V (ri)

)
+ g

∑
i<j

δ(ri − rj) (2.53)

Solving for the BEC’s wavefunction from this Hamiltonian is a challenging task. Thus, a

variational approach is employed, where the condensate wavefunction is approximated as ψ(r) =

√
Nϕ(r), and the particle density is n(r) = |ψ(r)|2. For large numbers of atoms, terms on the

order of 1/N are negligible, leading to the energy functional for the N-particle wavefunction:

E(ϕ, ϕ∗) = N

∫
d3r

(
ℏ2

2m
|∇ϕ(r)|2 + V (r)|ϕ(r)|2 + 1

2
Ng|ϕ(r)|4

)
(2.54)

The solution of the ground state is obtained by minimizing this energy functional under

variations of ϕ with the constraint
∫
d3r|ϕ(r)|2 = 1, maintaining a constant total number of

atoms. This minimization involves a Lagrange multiplier µ:

δ

[
E(ϕ, ϕ∗)− µN

∫
d3r|ϕ(r)|2

]
= 0, (2.55)

which leads to the time-independent Gross-Pitaevskii equation (GPE).

iℏ∂tϕ(r) =
[
− ℏ2

2m
∇2 + V (r) + gN |ϕ(r)|2

]
ϕ(r) = µϕ(r). (2.56)

31



This equation, a type of non-linear Schrödinger equation, incorporates the external potential V (r)

and a non-linear term gN |ϕ(r)|2 that describes the mean-field potential of other atoms. The

eigenvalue µ is the chemical potential, which differs from the mean energy per particle found in a

linear Schrödinger equation. Using the condensate wavefunction ψ(r) with normalized condition∫
d3r|ψ(r)| = N , the GPE becomes

iℏ∂tψ(r) =
[
− ℏ2

2m
∇2 + V (r) + g|ψ(r)|2

]
ψ(r) = µψ(r). (2.57)

2.5.2 Bogoliubov-de-Gennes excitations

The Bogoliubov-de-Gennes (BdG) equations, while initially formulated for superconductivity [39],

are also applicable in the analysis of Bose-Einstein condensates (BECs). In the context of BECs,

the BdG framework is utilized to examine the spectrum of elementary excitations. It provides

a mathematical structure to describe the behavior of small perturbations from the ground state

within the condensate. This application is particularly relevant in discovering the dispersion

relations of these excitations and understanding the interactions within the condensate. The BdG

equations offer a method for quantifying and analyzing the collective modes and the dynamics of

particles in the condensate, contributing to the theoretical understanding of various phenomena

observed in BECs, such as superfluidity and the formation of vortices.

From Eq. (2.57) we can obtain the ground stateψG and the corresponding chemical potential

µG =

∫
d3r ψ∗

G

[−ℏ2∇2

2m
+ V (r) + g|ψG|2

]
ψG. (2.58)

For a BEC, the wavefunction ψ with ”weak” excitations should also satisfy the GPE, thereby we

32



can decompose it into ψ = ψG + δψ. Under the BdG formalism, the excitation is expressed by u

and v with frequency ω, i.e.,

ψ(r, t) = exp(−iµGt) [ψG(r) + u(r) exp(−iωt)− v∗(r) exp(iω∗t)] . (2.59)

Inserting it into Eq. (2.57) and separating the terms with exp(−iωt) and exp(iω∗t), we obtain the

following BdG equations.

ωu =

[
− ℏ2

2m
∇2 + V (r) + 2g|ψG|2 − µG

]
u− gψ2

Gv, (2.60)

−ωv =

[
− ℏ2

2m
∇2 + V (r) + 2g|ψG|2 − µG

]
v − gψ∗

G
2u. (2.61)

This set of equations is linear and can be solved by diagonalizing.

33



Chapter 3: Review of Turbulence and Velocimetry

Turbulence, characterized as chaotic and unpredictable flow patterns in fluids, has intrigued

scientists for centuries as a result of its complexity and elusive nature. Horace Lamb once said:

”I am an old man now, and when I die and go to heaven there are two matters on which I

hope for enlightenment. One is quantum electrodynamics, and the other is the turbulent motion

of fluids. And about the former I am rather optimistic.” Werner Heisenberg said something

similar afterwards: ”When I meet God, I am going to ask him two questions: Why relativity?

And why turbulence? I really believe he will have an answer for the first.” These reflective

comments emphasize the mysterious nature of turbulence, illustrating its status as one of the

most challenging and unresolved puzzles in the scientific world.

The nonlinearity inherent in turbulent flows is a fundamental aspect that makes turbulence

a formidable challenge in fluid dynamics. Similar to all chaotic systems, turbulence is extremely

sensitive to the initial condition, making prediction a challenging task. Furthermore, when

attempting to model turbulence using mathematical equations, one encounters more unknowns

than available equations. This imbalance necessitates the introduction of approximations or

models to ”close” the system of equations, a significant challenge that brings uncertainties and

limitations to the predictive capabilities of turbulence models. The rich physics also lies in the

nonlinearity. The most profound findings in the turbulence are related to energy cascade, which is

34



vividly described by Lewis Fry Richardson’s poetry, ”Big whorls have little whorls, Which feed

on their velocity; And little whorls have lesser whorls, And so on to viscosity.” Energy injected

into the system on a large scale by external force such as stirring and shaking, and the kinetic

energy transferred from the large scale to small scales in a self-similar manner, as the ”big whorls

break into small whorls”, and finally it is dissipated on the small scale due to viscosity.

Although turbulence is characterized by randomness and unpredictability, there is a fundamental

aspect that exhibits statistical regularity: the scaling law for the kinetic energy spectrum. In

the subinertial range, where the energy cascade predominantly occurs, this scaling law becomes

prominent, standing as one of the most universal and consistent findings in the study of turbulence.

In 1941, Andrey Kolmogorov, a prominent mathematician in turbulence research, introduced

groundbreaking theories that helped unveil the statistical nature of turbulent flows [18, 19, 20].

Kolmogorov’s theories elucidated the multiscale nature of turbulence, shedding light on the

energy cascade process and the manifestation of the scaling laws in the kinetic energy spectrum.

His work laid a robust foundation for understanding the complexities inherent in turbulent flows,

steering subsequent research toward a more profound comprehension of this chaotic yet structured

phenomenon.

This chapter gives a review for the turbulence and velocimetry, from incompressible classical

turbulence to compressible turbulence as well as quantum turbulence. The Kolomogorov theory

for the incompressible turbulence is discussed in detail. Then the compressible turbulence is first

introduced under a coarse-grained picture. Finally, the quantum turbulence in both superfluid

Helium and atomic BEC is discussed. This chapter also gives a review for the velocimtery

techiniques in both classical and quantum fluids.

35



3.1 Incompressible Turbulence

In this section, the fluid considered is incompressible fluid, so the density of the fluid is

constant and the velocity field satisfies ∇ · v = 0 due to continuity equation.

3.1.1 Basic features of turbulence

Before discussing the details of the features, it is necessary to introduce the concept of

”eddy”, which is defined as the region of significant self-correlation. If we stir a fluid to create

turbulence, the largest eddy would be the system size, where it shows the macroscopic flow.

In Richardson’s poems, this large eddy then breaks into smaller eddies which have smaller

correlation length, and finally to a scale that dissipation dominates.

3.1.1.1 Irregular, random and chaotic

The most annoying feature of turbulence is that its velocity field v(x, t) varies significantly

and irregularly. Suppose that you place a velocity sensor in the fluid and record the velocity.

You would obtain a trace as in the Fig. 3.1. And if you repeat the measurement, the traces are

unrepeatable even though you set all the experimental parameters the same, which means that at

least some of the parameters are not actually under control.

36



0 2 4 6 8 10
−10

−5

0

5

10

Time (s)

V
el

oc
ity

(m
m

/s
)

Velocity Field vs Time

Figure 3.1: A representation of the chaotic and random nature of the velocity field over time.

3.1.1.2 High Reynolds number

The Reynolds number (Re) [40] is a dimensionless quantity that is used to predict flow

patterns in fluid mechanics. It is defined as

Re =
UL

ν
, (3.1)

where U is the characteristic velocity, L is the characteristic length, and ν is the kinematic

viscosity.

The Reynolds number is actually the ratio of convection to diffusion, or the inertial force to

viscous force, which can be seen from the Navier-Stokes equation [41] for incompressible fluid

∂v

∂t
+ (v · ∇)v = −1

ρ
∇p+ ν∇2v + fext, (3.2)

37



where ρ is the density, p is the pressure, and fext is the external force. The term (v · ∇)v is the

inertial force, which also describes the convection of the momentum. And ν∇2v is the viscous

force, also known as the diffusion force. The ratio of these two terms is aprroximately to be UL
ν

in the scale you are interested in. So high Re (Re ≫ 1) typically stands for the case when the

scale is set to the large scale, where the dissipation from the viscous force can be ignored.

3.1.1.3 Dissipation

Dissipation occurs always in the fluid. On a small scale, dissipation dominates (Re ≪ 1).

For the steady state, the energy injection rate must be equal to the dissipation rate

ϵ ≈ u20
τ0

=
u30
L0

, (3.3)

where u0, L0 are the charactersitic velocity and length for the largest eddy, τ0 = L0

u0
is its

characteristic time. Note that only in the steady state we can write the dissipation rate in a form

independent of ν.

The dissipation range η must be related to the dissipation rate ϵ and the viscosity ν. Dimensional

analysis gives η = (ν
3

ϵ
)
1
4 .

3.1.1.4 Efficient mixing

For a regular nonturbulent fluid, suppose that at t = 0 we have two closely located tracer

particles, its relative distance |∆x| at time twould satisfy ⟨|∆x|2⟩ ∝ Dt, whereD is the diffusion

coefficient. This
√
t diffusion rate can be derived from Fick’s second law of diffusion [42].

In turbulence, since dissipation (viscosity) is closely related to diffusion, |∆x|2 must be a

38



function of the dissipation rate ϵ and t. From the dimensional analysis we have ⟨|∆x|2⟩ = ϵt3.

The chaotic motion in the turbulent fluid accelerates diffusion. This t3/2 suggests superdiffusion

in turbulence, which is called Richardson’s law [43].

3.1.1.5 Vortices

A turbulent fluid is always combined with vortices. Vortices in incompressible turbulence

are fundamental elements, marking regions where the fluid circulates around a core area. Understanding

the dynamics of these vortices is crucial in the study of incompressible turbulence, as they are

integral to the energy transfer processes and the overall chaotic motion observed in such fluid

flows.

3.1.2 The Reynolds-Averaged Navier–Stokes Equations

Reynolds decomposition [44] involves breaking down the instantaneous flow variables into

mean and fluctuating components. Consider a flow variable, such as velocity u, which varies with

time and space due to turbulence. It can be decomposed as follows: u(x, t) = ū(x) + u′(x, t),

where ū(x) is the velocity field averaged over time and u′(x, t) is the fluctuating part. This

decomposition simplifies the equations of motion and aids in understanding and analyzing the

turbulent characteristics of the flow. Applying the Reynolds decomposition to the Navier-Stokes

equations [Eq. (3.2)] leads to the Reynolds-averaged Navier-Stokes equations (RANS).

∂(v̄ + v′)

∂t
+ ((v̄ + v′) · ∇) (v̄ + v′) = −1

ρ
∇(p̄+ p′) + ν∇2(v̄ + v′) + fext, (3.4)

39



Time average to the entire equation, we obtain

v̄ · ∇v̄ + v · ∇v′ = −1

ρ
∇p̄+ ν∇2v̄ + fext. (3.5)

For incompressible fluid, we have ∇ · v̄ = 0 and ∇ ·v′ = 0, so the equation can be written as the

RANS form, that is,

∇ · v̄2 +∇ · v′2 = −1

ρ
∇p̄+ ν∇2v̄ + fext, (3.6)

where on the left hand side v2 is a short hand notation for the rank 2 tensor with elements v2
ij =

vivj . And the v′2 is called the Reynolds stress tensor. Note that the above equation contains

more unknowns than the number of equations. The equation of Reynolds stress tensor actually

introduces a higher-order tensor of v′, and the number of unknowns is always greater than the

number of equations. This difficulty is always called the ”Turbulence closure problem”.

3.1.3 Vortex Stretching

Vorticity is defined by ω = ∇× v. From Eq. (3.2), we have

∂v

∂t
+

1

2
∇(v · v)− v × ω = −1

ρ
∇p+ ν∇2v + fext, (3.7)

where we have used 1
2
∇(v ·v) = (v ·∇)v+v×(∇×v). Now take the curl of the entire equation

∂ω

∂t
−∇× (v × ω) = ν∇2ω +∇× fext. (3.8)

40



The second term on the left hand side (LHS) can be expanded

∇× (v × ω) = (ω · ∇)v − (v · ∇)ω + v∇ · ω − ω∇ · v, (3.9)

where the last two terms on the right hand side (RHS) are 0. We then obtain the vorticity transport

equation

∂ω

∂t
+ (v · ∇)ω = (ω · ∇)v + ν∇2ω +∇× fext. (3.10)

The first term on the RHS represents ”vortex streching”, i.e., velocity gradient leads to the change

in the rate of vorticity. The second term is the diffusion term, and we observe that the kinematic

viscosity is actually the diffusion coefficient for the diffusion of vorticity.

The vortex streching is once believed to be the cause of energy cascade from an angular

momentum conservation argument; however, it might not be the real cause [45], so I will only

briefly state its implication to the cascade here. Consider the fluid element with vorticity ω,

the fluid element tends to be more elongated under strain, and the diameter of the fluid element

decreases. Assuming that the local angular momentum is conserved, the vorticity has to increase.

Thus, the kinetic energy is transported from large-scale to small-scale.

3.1.4 Equilibrium range theories

Kolmogorov’s theory assumes homogeneous and isotropic flow, which might be quite

different from the real cases, but it is a good starting point and gives us scaling laws agrees

well even with some complicated inhomogenous and anisotropic flow.

Much similar to Richardson’s poet, Kolmogorov explains the idea more rigorously. It is

41



built on two similarity hypothesis.

3.1.4.1 Hypothesis of Kolmogorov

Kolmogorov presented a theory in which large anisotropic eddies are the primary energy

carriers in turbulent flows, transferring energy across various scales. These eddies gradually lose

their distinct shapes and characteristics, transitioning into a state where they are homogenous and

isotropic. In this state, the energy of the eddies depends solely on the energy they receive from

the larger eddies and the energy dissipation rate due to the smaller eddies.

Kolmogorov’s first similarity hypothesis: At large Reynolds numbers, the local average properties

of the small-scale components of any turbulent motion are determined entirely by kinematic

viscosity and average rate of dissipation per unit mass.

Following his hypothesis, we can derive the characteristic length η, time τη and velocity uη for

the ”small scale”, since it only depends on ϵ and ν, i.e.,

η =

(
ν3

ϵ

) 1
4

, τη =
(ν
ϵ

) 1
2
, uη = (ϵν)

1
4 . (3.11)

It is easy to verify that the Reynolds number of the scale η is 1, which suggests a balance of the

convection to diffusion.

Kolmogorov’s second similarity hypothesis: There is an upper subrange (the inertial subrange)

in this bandwidth of small eddies in which the local average properties are determined only

by the rate of dissipation per unit mass.

42



On this scale, both energy injection and dissipation are negligible, and the kinetic energy is

transferred from the large scale to the small scale. In steady state, the energy injected rate should

be equal to the dissipation rate, and on a large scale, the dissipation is negligible, so it depends

only on the energy injected rate ϵ and not on the detail of the dissipation mechanism given by ν.

In the inertial subrange, the energy spectrum depends only on k and ϵ. The kinetic energy

has the unit of L3τ−2. ϵ has the unit of L2τ−3. From dimensional analysis we obtain E(k) =

Cϵ2/3k−5/3, where C is a constant and lies in the range of 1.5 to 2.5 for most fluids.

3.1.4.2 Velocity Structure function

The energy spectrum is not a quantity that can be measured with ease, so realistically it

is worth using a more experimental friendly metric, that is, the velocity field structure function

Sp(l), which is defined in Eq. (3.12) based on the velocity increments δvl,

δvl(x⃗) = ⟨|v⃗(x⃗+ l⃗)− v⃗((x⃗))|⟩ · l̂, Sp(l) = ⟨|δvl(x⃗)|p⟩; (3.12)

where l⃗ is the positional separation vector of the two velocity field sample points, and l̂ is its unit

vector, and p is the order of the structure function. Note that the structure function can also be

defined without the absolute value, and Kolmogorov’s four-fifth law is built upon the S3 of the

longitudinal velocity increments without taking the absolute value. On the experimental side, it

is more convenient to use the absolute value, so the gap between the magnitude of the odd and

even order of structure function is less prominent.

Historically, experimentalists only use a single velocity sensor to read the streamline velocity

43



on a single spot, record its readings at different times, and then map the temporal difference to

the spatial difference making use of the Galilean invariance. The second-order structure function

was found to follow a scaling of 2/3 over l, i.e. S2(l) ∝ l2/3. This empirical law is often called

the ”two-thirds law”, which seems to hold for all turbulent flow, at least approximately.

From Kolmogrov’s second similarity hypothesis, the velocity structure function should

solely depends on the dissipation rate ϵ if the relative distance l is in the inertial subrange, then

using dimensional analysis, Sp has the unit of Lpτ−p, thus

Sp(l) ∝ (ϵl)p/3 , (3.13)

which agrees with the two-thirds law when p = 2.

In 1941, Kolmogorov gave an exact result for the third order longitudinal structure function

S3(l) = −4

5
ϵl, (3.14)

which is called the ”four-fifths law”. I will introduce Kolmogorov’s derivation in the following

section.

44



3.1.4.3 Karman-Howarth-Monin relation

Kolmogorov derives the four-fifths law based on the Karman-Howarth-Monin relation [46],

i.e.,

∂t⟨v(x) ·v(x+l)⟩ = 1

2
∇l ·⟨δv2δv⟩+⟨v(x) ·(fext(x+ l) + fext(x− l))⟩+2ν∇2

l ⟨v(x) ·v(x+l)⟩.

(3.15)

Proof: We will use the shorthand notation prime to denote the variables at x + l. From

Eq. (3.2), we have

∂tvi = −∂j(vjvi)− ∂iP + fext,i + ν∂j∂jvi, (3.16)

and an equation for primed variable. Then we have

∂t⟨viv′i⟩ =− ⟨v′i∂j(vjvi)⟩ − ⟨vi∂′j(v′jv′i)⟩ − ⟨v′i∂iP ⟩ − ⟨vi∂′iP ′⟩

+ ⟨v′ifext,i⟩+ ⟨vif ′
ext,i⟩+ ν⟨v′i∂j∂jvi⟩+ ν⟨vi∂′j∂′jv′i⟩,

(3.17)

where the last two terms on the first line is zero due to integral by parts and incompressibility. The

viscous terms can be expressed by the derivative of the separation distance l, that is, ν⟨v′i∂j∂jvi⟩ =

ν⟨vi(x)∂j∂jvi(x − l)⟩ = ν∇2
l ⟨viv′i⟩, where homogeneity is used. The forcing term can be

expressed as ⟨vi(x)(fext,i(x+ l) + fext,i(x− l))⟩.

What remains to prove is that the first two terms in the RHS of the above equation are equal

to the first term on the RHS of Eq. (3.15).

⟨v′i∂j(vjvi)⟩+ ⟨vi∂′j(v′jv′i)⟩ = −∂lj⟨v′ivjvi⟩+ ∂lj⟨viv′jv′i⟩ = ∂lj⟨v′iviδvj⟩. (3.18)

45



Then using v′ivi = (vivi + v′iv
′
i − δviδvi)/2,

∂lj⟨v′iviδvj⟩ =
1

2
⟨v2i ∂ljδvj⟩+

1

2
∂lj⟨v′iv′iδvj⟩ −

1

2
∂lj⟨δviδviδvj⟩, (3.19)

where ⟨v2i ∂ljδvj⟩ = ⟨v2i ∂jv′j⟩ = 0 because of incompressibility, and ∂lj⟨v′iv′iδvj⟩ = −⟨vivi∂ljvj(x−

l)⟩ = ⟨vivi∂jvj(x − l)⟩ = 0 because of homogenity and incompressibility. Thereby we obtain

Eq. (3.15).

3.1.4.4 Derivation of four-fifths law

Considering the fully developed turbulence (ν → 0) and in the stationary limit (t → ∞),

the time dependence of Eq. (3.15) is removed. In the inertial subrange, the correlation length

of the external force is much greater than the separation l, so we can write the forcing term to

2⟨v(x)fext(x)⟩ = 2ϵ, where ϵ is the energy injection rate, which is also equal to the dissipation

rate in steady state. In the high Renyolds number limit (ν → 0), the viscous term is negligible in

the inertial subrange. Thus we have

∇l · ⟨δv2δv⟩ = −4ϵ. (3.20)

Now we consider the third-order longitudinal structure function S3(l), where the longitudinal

notation is removed for simplicity.

S3(l) =
liljlk
l3

⟨δviδvjδvk⟩ = 2
liljlk
l3

(Sij,k + Ski,j + Sjk,i), (3.21)

46



where Sij,k = ⟨vivjv′k⟩. Using homogeneity and isotropy, we can expand the tensor Sij,k to

Sij,k(l) = A(l)δij
lk
l
+B(l)

(
δki
lj
l
+ δjk

li
l

)
+ C(l)

liljlk
l3

. (3.22)

And S3(l) becomes a function solely depends on A(l), B(l) and C(l), i.e.,

S3(l) = 2
liljlk
l3

[
(A+ 2B)

(
δij
lk
l
+ δjk

li
l
+ δki

lj
l

)
+ 3C

liljlk
l3

]
= 6(A+ 2B) + 6C. (3.23)

Now we consider connecting the structure function to the Eq. (3.20), using incompressibility

∇l · ⟨δv2δv⟩ = 4∂jSij,i = 4∂j

[
(A+ 4B + C)

lj
l

]
. (3.24)

Using the relation ∂kf(l) = (lk/l)f
′(l) and ∂k(li/l) = δik/l − lilk/l

3, we have

∇l · ⟨δv2δv⟩ = 4(A′ + 4B′ + C ′) +
8

l
(A+ 4B + C), (3.25)

which links the A, B and C with the dissipation rate ϵ but introducing new differential variables.

Using the incompressibility of ∂kSij,k = 0 to obtain more relations,

∂kSij,k =

[
A′(l) + 2

A(l) +B(l)

l

]
δij +

[
2B′(l) + C ′(l) + 2

C(l)−B(l)

l

]
lilj
l2

= 0, (3.26)

where the two terms in the square bracket should be both vanishing. Multiplying the first bracket

47



by 3 and adding the second bracket, we obtain the following.

d

dl
[3A+ 2B + C] +

2

l
[3A+ 2B + C] = 0. (3.27)

The general solution is 3A + 2B + C = const/l2, since the solution should be finite at l → 0,

we obtain

3A+ 2B + C = 0, B = −A− l

2
A′, C = −A+ lA′. (3.28)

Eq. (3.23) and Eq. (3.25) become

S3(l) = −12A, −ϵ = −lA′′ − 7A′ − 8

l
A. (3.29)

Under the change of variable y = A/l, x = ln l, the differential equation becomes

y′′ + 6y′ + 15y = ϵ. (3.30)

The general solution is as follows.

y = α exp [(−3 + 6j)x] + β exp [(−3− 6j)x] +
ϵ

15
. (3.31)

Under the limit of l → 0 (x→ −∞), y should be finite, so y = ϵ/15 and

S3(l) = −4

5
ϵl. (3.32)

48



3.1.4.5 Velocity increments PDF and the intermittency

The velocity increments probability distribution function (PDF) is an essential metric for

intermittency, which describes the anomolous scaling of the velocity increments compared to the

Kolmogorov theory.

Kolmogorov theory is based on the self-similarity assumption, which leads to scale-invariant

scalings for the longitudinal structure function. So, if we study the PDF of the longitudinal

velocity increments δv(l), it should also be scale-invariant, which means the large-scale velocity

increment PDF should have the same shape as the small-scale velocity increment PDF in the

inertial range. However, it is found that the fat-tailed feature becomes more and more clear when

the spatial separation l is decreased. When the separation is large, the PDF is Guassian or near

Guassian, but when the separation is small, the rare event with large velocity increments becomes

more and more frequent. This scale dependent phenomenon is the symbol of intermittency, which

suggests that the Kolmogorov theory needs revision to the self-similarity hypothesis.

In 1962, Kolmogorov and Obukhov gave a refined the self-similarity hypothesis [23]. The

original averaged dissipation rate ϵ is replaced by a scale-dependent dissipation rate ϵl, which can

be viewed as the filtered version of the dissipation field. The previous K41 calculation is still

assumed to be valid for ϵl, but the whole system’s expectation value should depend on the PDF

of ϵl, which is assumed to follow a log-normal distribution. Then a refined structure function is

given by

Sn(l) ∝ (E[ϵl]l)
n/3(l/L)−µn(n−3)/18, (3.33)

where E stands for the expectation value, L is the system size and µ is experimentally measured to

49



be close to 0.23. So this anomalous scaling is apparent for high-order velocity structure function,

and for low-order the correction is negligible since the rare event’s effect to the structure function

is not prominent.

3.2 Compressible turbulence

The previous sections mainly discuss incompressible turbulence, where the density of the

fluid is spatially independent, so we have the incompressible condition ∇ · v = 0. However,

realistically, the fluid is compressible, especially when the fluid’s velocity is higher than the

speed of sound. In this section we will discuss the compressible fluid, where the density now

becomes a spatial-dependent variable.

3.2.1 Governing dynamics

From conservation law, we can write the transport equation of density, momentum, and

kinetic energy [47], i.e.,

∂tρ+ ∂j(ρvj) = 0, (3.34)

∂t(ρvi) + ∂j(ρvivj) = −∂iP + ∂jσij + ρfext,i, (3.35)

∂t(ρ
|v|2
2

) + ∂j

[
(ρ
|v|2
2

+ P )vj − 2µ(viSij −
1

3
vjSkk)

]
= P∂jvj − 2µ

(
|Sij|2 −

1

3
|Skk|2

)
+ ρvifext,i,

(3.36)

where Eq. (3.35) is the compressible Navier-Stokes equation. µ is the dynamic viscosity, Sij =

(∂jvi + ∂ivj)/2 is the symmetric strain tensor, and σij = 2µ(Sij − 1
3
Skkδij) is the viscous stress.

Note that the Stokes hypothesis is used so that the bulk viscosity is ignored.

50



3.2.2 Coarse graining

Following the coarse graining and filtered appraoch used in the large eddy simulation, we

can decompose any field a(x) into the large-scale component āl(x) and sub-grid component

a′
l(x), where the coarse-grained field contains modes at scales greater than l, i.e.,

āl(x) =

∫
d3rGl(r)a(x+ r), (3.37)

where G(r) is the Friedrichs mollifier satisfies normalized condition and centered at r = 0. And

Gl(r) = l−3G(r/l). Note that this filtering operation is linear and commutes with spatial and

time derivatives. The residue field is the sub-grid component

a′
l(x) = a(x)− āl(x). (3.38)

Applying the filtering operation to the continuity and Navier-Stokes equation, we have

∂tρ̄l + ∂j(ρvj)l = 0. (3.39)

∂t(ρvi)l + ∂j(ρvivj)l = −∂iP l + ∂jσij l + (ρfext,i)l (3.40)

To reduce number of terms, we use the density-weighted Favre average [48], i.e.,

ãl =
(ρa)l
ρl

. (3.41)

We will remove the subscript l from now on when there is no confusion. Then the continuity and

51



Navier-Stokes equation become

∂tρ̄+ ∂j(ρṽj) = 0. (3.42)

∂t(ρṽ) + ∂j(ρṽiṽj) = −∂j(ρτ̃(vi, vj))− ∂iP + ∂jσij + ρf̃ext,i, (3.43)

where τ̃(vi, vj) = ṽivj − ṽiṽj is the turbulent stress accounting for the effect from the scale

smaller than l, and it is the only term couples the large scale momentum to the small scale. In

addition, we can define the 2nd-order generalized central moments of any fields f(x) and g(x)

by τ(f, g) = (fg)l−f lgl. Using τ and τ̃ to describe the small-scale effect, the transport equation

for large-scale kinetic energy can be written as

∂t(ρ
|ṽ|2
2

) + ∂jJj = −Πl − Λl + P∂jvj −Dl + ϵinjl , (3.44)

where

Πl(x) = −ρ∂jviτ̃(vi, vj), Λl(x) =
1

ρ
∂jPτ(ρ, vj) (3.45)

are the subgrid scale kinetic energy flux to scales less than l,

Dl(x) = ∂jũi[2µSjj −
2

3
µSkkδij], ϵinjl (x) = ṽiρf̃ext,i (3.46)

are the dissipation and energy injection by the external force on the large scale,

J(x) = ρ
|ṽ|2
2
ṽj + Pvj + ṽiρτ̃(vi, vj)− ṽiσij. (3.47)

is the large-scale transport of kinetic energy, and −P∂jvj is the large-scale pressure dilatation.

52



3.2.3 Scale Decomposition

The viscous range lµ is defined as the scale at which viscous effects become significant in

kinetic energy balance. In [47] it has been shown that Dl is negligible on scales l ≫ lµ, and

injection of kinetic energy can be localized on large scales L≫ l with proper stirring. Therefore,

there exists an intermediate scale L ≫ l ≫ lµ where dissipation and external injection are

negligible, which is the inertial range of compressible turbulence.

The subgrid scale flux terms in Eq. (3.45) transfer large-scale kinetic energy to scales

smaller than l. H. Aluie argues in [49] that this subgrid scale flux cascade is localized based on

empirical facts that the structure function is localized and the scalings are weak. This result has

been verified by a numeric study [50]. If we use a density-weighted velocity field, the statistics

of the compressible turbulence is very similar to the incompressible turbulence.

3.3 Quantum turbulence

Quantum turbulence refers to the turbulence of quantum fluid or superfluid, where the fluid

is purely inviscid or mixed with some portion of viscid classical fluid. The quantum fluid in the

literature emphasizes quantized vortices more than the invisicid fluid. But normally we do not

distinguish them very carefully, since the quantized vortices come from the irrotational nature of

superfluid.

In this section, I will first introduce the history of superfluid and theoretical models of its

fluid dynamics, then discuss the energy cascade and dissipation mechanism.

53



3.3.1 Superfluid

The superfluid, named by Pyotr Kapitsa, was first observed by him in 1937 in liquid

4He [51]. When the liquid 4He is cooled below the λ point (2.17 K), he observed that the liquid

could flow with no resistance through the capillaries and over the edges of containers. He also

observed that a thin layer of helium below the λ point (helium II) could climb up and out of

the container. Around the same time, John Allen and Don Misener observed an unusually high

thermal conductivity [52]. A