ABSTRACT

Title of Dissertation: FLOQUET HEATING AND RELAXATION
OF INTERACTING BOSE EINSTEIN CONDENSATES

James Maslek
Doctor of Philosophy, 2022

Dissertation Directed by: Professor Trey Porto
Department of Physics

Floquet’s theorem says that any unitary, periodically driven system can be described by an

effective time-independent Hamiltonian, where the effective Hamiltonian can have completely

different properties than the static, undriven system. Floquet engineering makes use of this idea

to simulate new Hamiltonians that would otherwise not be possible in the undriven case. For

interacting systems, this approach can be used to realize interesting correlated many-body states,

but drive-induced heating must be understood and mitigated. Cold atoms in optical lattices

provide a controllable, well-isolated system in which these ideas can and have been realized. I

describe research into two areas of Floquet engineering for interacting Bose-Einstein condensates

in periodically driven optical lattices.

The first half of this thesis focuses on the study of heating mechanisms for condensates

in periodically driven lattices. In the weakly interacting limit, one might expect that heating

could be described with a Fermi Golden Rule approach. Parametric driving of fluctuations in the



condensate, however, can lead to runaway heating that cannot be described perturbatively. We ex-

perimentally study heating in shaken 2D square lattices and demonstrate heating consistent with

the theoretical predictions of parametric instabilities.

The second half of this thesis describes experiments that realize Floquet-induced effec-

tive staggered magnetic fields, and the relaxation dynamics of interacting particles subject to

these fields. Interestingly, we observe pre-thermal relaxation dynamics, where an initially heated

cloud suddenly subject to the effective Hamiltonian condenses into a state governed by the drive-

induced effective Hamilton on a timescale faster than heating.



FLOQUET HEATING AND RELAXATION
IN INTERACTING BOSE EINSTEIN CONDENSATES

by

James Maslek

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

2022

Advisory Committee:
Professor Steven Rolston, Chair
Professor Trey Porto, Co-Chair/Advisor
Professor Gretchen Campbell
Professor Victor Galitski.
Professor Thomas Murphy (Dean’s Representative)



© Copyright by
James Maslek

2022



Acknowledgments

I would like to start off by thanking Trey for the opportunity to join the Rubidium 1 when I

was still an undergrad at Rochester. I appreciate your guidance and support throughout the years

and I will always be amazed that you turned your house into a sundial. I will take to my grave

that python is superior to Mathematica. Additionally, thank you to Steve, Gretchen and Alicia

for help and mentorship through the CPRK group.

I would like to acknowledge my lab mate Carlos, it was a pleasure to work with over the

last many years. From building 80/20 in an empty room and figuring out if a laser that fell in a car

was broke to jerry-rigging the vacuum shutter minutes before opening it up, I enjoyed learning

and growing through grad school with you. Additionally, throughout my time, I have had the

pleasure to work with a great bunch of people on the Rb1 experiment: Thomas, Eric, Elizabeth,

Daniel, AJ, Pierre.

I would like to thank my family for their continued support over the last 20 years of school.

Finally, to my wife, Marisa. It was fun getting to have two weddings over the pandemic,

but you were always there for me. And you let me get a doggo, Cooper, so thanks for that.

Go Bills.

ii



Table of Contents

Acknowledgements ii

Table of Contents iii

List of Figures v

Chapter 1: Introduction 1

Chapter 2: Background 4
2.1 Floquet Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Time Evolution Operator . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.3 High Frequency Expansion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.4 Mean Field Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.5 Weakly Interacting Condensates : GP Limit . . . . . . . . . . . . . . . . . . . . 9
2.6 Bogoliubov Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.7 Optical Lattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14

2.7.1 Band Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
2.7.2 2D Checkerboard Band Structure . . . . . . . . . . . . . . . . . . . . . . 20

2.8 Bose-Hubbard Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
2.8.1 Checkerboard Tight Binding . . . . . . . . . . . . . . . . . . . . . . . . 28

2.9 Artificial Magnetic Fields on a Lattice . . . . . . . . . . . . . . . . . . . . . . . 31
2.9.1 Harper-Hofstadter Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . 32

2.10 Periodically Driven Optical Lattices . . . . . . . . . . . . . . . . . . . . . . . . 33
2.10.1 Tight Binding Floquet Treatment . . . . . . . . . . . . . . . . . . . . . . 33
2.10.2 Trotter Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
2.10.3 Extended Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39

Chapter 3: Experiment 44
3.1 Vacuum System . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.1.1 Atom Source and Oven . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
3.1.2 UHV Side . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

3.2 Magnetic Fields . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
3.3 Laser Systems . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.3.1 Cooling Light . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.4 Microwave and RF . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
3.5 Imaging . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

iii



3.6 Computer Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
3.7 Optical Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.7.1 Optical Lattice Calibration . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.7.2 Tilt Calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

3.8 Tripod Piezo Mounts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
3.8.1 Lattice Displacement . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.8.2 Design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
3.8.3 Piezo Characterization . . . . . . . . . . . . . . . . . . . . . . . . . . . 60

Chapter 4: Parametric Heating in Periodically Driven Optical Lattice 66
4.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

4.1.1 Driving in 2D . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.2 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
4.3 Experimental Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78

4.3.1 Lattice Depth Dependence . . . . . . . . . . . . . . . . . . . . . . . . . 78
4.3.2 Drive Amplitude Dependence . . . . . . . . . . . . . . . . . . . . . . . 79
4.3.3 Frequency Dependence . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.3.4 Comparison With Fermi Golden Rule . . . . . . . . . . . . . . . . . . . 86

Chapter 5: Relaxation into Floquet Prethermal Condensates 88
5.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

5.1.1 Next Nearest Neighbor Tunneling . . . . . . . . . . . . . . . . . . . . . 95
5.1.2 Extended Basis Picture . . . . . . . . . . . . . . . . . . . . . . . . . . . 97
5.1.3 Mapping Tight Binding to Extended Basis . . . . . . . . . . . . . . . . . 101

5.2 Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.2.1 Adiabatic Loading . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104
5.2.2 Heating from Band Sweeps . . . . . . . . . . . . . . . . . . . . . . . . . 109
5.2.3 Pre-thermal Relaxation . . . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.2.4 Relaxation Between Band Minima . . . . . . . . . . . . . . . . . . . . . 116

Chapter 6: Outlook 119
6.1 Topology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119

6.1.1 Theory . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 119
6.1.2 Experimental Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

Bibliography 127

iv



List of Figures

2.1 Bogoliubov mode dispersion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.2 Wannier functions for 1D optical lattice . . . . . . . . . . . . . . . . . . . . . . 19
2.3 Dynamically controllable 2D optical lattice . . . . . . . . . . . . . . . . . . . . 21
2.4 Matrix plot of Hamiltonian in plane wave basis . . . . . . . . . . . . . . . . . . 23
2.5 2D Band Structure for checkerboard lattice . . . . . . . . . . . . . . . . . . . . . 24
2.6 Band structure along symmetry line for different out of plane lattice depths . . . . 25
2.7 Superfluid to Mott insulator transition . . . . . . . . . . . . . . . . . . . . . . . 28
2.8 Tight Binding approximation to plane wave Hamiltonian. . . . . . . . . . . . . . 30
2.9 Floquet modification of static band structure . . . . . . . . . . . . . . . . . . . . 35
2.10 Band structure for amplitude modulated optical lattice with the Trotter method . . 39
2.11 Band structure for shaken lattice in the extended basis. . . . . . . . . . . . . . . 43

3.1 87Rb atomic levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
3.2 Experimental µW spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48
3.3 Optical lattice configuration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
3.4 Optical lattice diffraction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52
3.5 Piezo layout for driven optical lattice . . . . . . . . . . . . . . . . . . . . . . . . 54
3.6 Flat mirror displacement geometry . . . . . . . . . . . . . . . . . . . . . . . . . 55
3.7 Lattice shaking trajectories . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
3.8 Tripod piezo mirror design . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59
3.9 Electronics board for piezo control . . . . . . . . . . . . . . . . . . . . . . . . . 60
3.10 Mechanical resonance experiment and calibration . . . . . . . . . . . . . . . . . 61
3.11 Piezo displacement calibration . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
3.12 Frequency dependent piezo displacement . . . . . . . . . . . . . . . . . . . . . . 63
3.13 DC Piezo steering angle . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64
3.14 Piezo calibration via lattice depth measurement . . . . . . . . . . . . . . . . . . 65

4.1 Lattice driving procedure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77
4.2 Exponential decay of condensate fraction . . . . . . . . . . . . . . . . . . . . . 78
4.3 Dependence of heating rate on lattice depth . . . . . . . . . . . . . . . . . . . . 79
4.4 Heating rate dependence on drive strength . . . . . . . . . . . . . . . . . . . . . 81
4.5 Single experimental run for Kc

0 . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.6 Critical drive strength vs frequency . . . . . . . . . . . . . . . . . . . . . . . . . 84
4.7 Frequency dependence of heating rates . . . . . . . . . . . . . . . . . . . . . . . 85

v



4.8 ω−1 vs ω−2 dependence of heating rates . . . . . . . . . . . . . . . . . . . . . . 86

5.1 Effective staggered magnetic field setup . . . . . . . . . . . . . . . . . . . . . . 91
5.2 Structure factors arising from Tight Binding description . . . . . . . . . . . . . . 94
5.3 Frequency dependence of extended basis calculation . . . . . . . . . . . . . . . . 99
5.4 K0 shifts of the energy levels . . . . . . . . . . . . . . . . . . . . . . . . . . . . 100
5.5 Comparing extended basis calculations to Tight Binding . . . . . . . . . . . . . . 103
5.6 Adiabatic loading into effective flux configurations . . . . . . . . . . . . . . . . 105
5.7 X± condensate dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107
5.8 K0 dependence of the 2π flux model . . . . . . . . . . . . . . . . . . . . . . . . 108
5.9 Thermalization dependence on ramp time . . . . . . . . . . . . . . . . . . . . . 109
5.10 Landau-Zener sweeps starting below resonance . . . . . . . . . . . . . . . . . . 110
5.11 Landau-Zener sweeps starting above resonance . . . . . . . . . . . . . . . . . . 111
5.12 Coupling to higher bands at q = Γ . . . . . . . . . . . . . . . . . . . . . . . . . 112
5.13 Prethermal relaxation for the staggered 2π flux . . . . . . . . . . . . . . . . . . . 113
5.14 Prethermal relaxation for the staggered π flux . . . . . . . . . . . . . . . . . . . 115
5.15 Single shot image for data analysis . . . . . . . . . . . . . . . . . . . . . . . . . 116
5.16 Elliptical Band Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117
5.17 Relaxation between X+ and X− . . . . . . . . . . . . . . . . . . . . . . . . . . 118

6.1 Dirac points in the extended basis picture . . . . . . . . . . . . . . . . . . . . . 121
6.2 Raw images of q dependent oscillations . . . . . . . . . . . . . . . . . . . . . . 124
6.3 q dependent oscillations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 125
6.4 Amplitude and phase of Bloch vector reconstruction. . . . . . . . . . . . . . . . 126

vi



Chapter 1: Introduction

Floquet engineering is a powerful tool that makes use of a periodic modulation to simulate

an effective time independent system, where the resulting time independent system can exhibit

behaviors that are not present in the static system. Cold atoms in optical lattices provide a con-

trollable, well isolated system in which these ideas can and have been realized. Circular driving

of a honeycomb lattice has been performed to realize the Haldene model [1], which exhibits a

topological band structure. Periodic driving in optical superlattices has been used to simulate the

Harper-Hofstadter model with neutral atoms [2, 3, 4]. The simulation of the Harper-Hofstadter

Hamiltonian is a specific realization of artificial gauge fields[5], which have been studied in Flo-

quet engineered systems [6, 7, 8, 9]. Additionally, Floquet engineering has been proposed as a

method to generate fractional quantum hall states [10, 11], and experimentally simulate Z2 lattice

gauge theories [12].

Despite the interesting new possibilities that emerge in time periodic systems, the addition

of energy, particularly in interacting systems, inevitably results in heating, which in itself is an

outstanding theoretical and experimental area of investigation. Although the continual addition

of energy eventually heats the system to infinite temperature [13, 14], there can exist a range

of parameters for which there is different behavior on intermediate timescales. These “prether-

mal” states emerge as an approximate equilibrium condition where the system has thermalized

1



according to the effective Hamiltonian associated with the periodic modulation, which is dif-

ferent from the full time dependent Hamiltonian [15, 16, 17]. For lattice systems with locally

bounded Hilbert spaces, such “Floquet pre-thermalization” is generally expected in the limit of

high frequency driving [15, 18, 19, 20]. While the Hilbert spaces of real systems are not locally

bounded, if the states of interest are separated from all other states by a sufficient energy gap,

Floquet pre-thermalization can occur [21].

In this thesis, I describe research into two areas of Floquet engineering for Bose-Einstein

condensates in periodically driven optical lattices. The first is a mechanism of heating on conden-

sates in an optical lattice. Previous theoretical work presented two main descriptions of heating

dynamics for bosons in driven optical lattices. One is a Fermi golden rule approach where the

periodic driving of the lattice couples atoms to many other states, which is expected to be ac-

curate for weakly heating systems. The second, which is the focus of our work, is applicable

for interacting mean-field systems with bosonic excitations. This strong heating arises due to

the existence of parametric instabilities in the Floquet system. Driving resonantly excites modes,

which can then grow exponentially under the right conditions and result in a depletion of the

condensate.

The second aspect of this thesis is the realization of effective staggered fields in a Floquet

system and the study of interaction driven dynamics in the effective system. Circularly driving

the optical lattice, under the right conditions, leads to an effective system with complex tunneling

parameters that can be described by an effective magnetic field. The effective magnetic field is

similar to the Harper-Hofstadter model, which describes particles in a lattice subject to a gauge

field. The Harper-Hofstadter model was motivated by the band structure of an electron on a lattice

in a magnetic field, which results in the fractal Hofstadter butterfly band structure [22]. One of

2



the signatures of the gauge field is that the bands display minima at non-zero quasimomentum. In

our work, we study the relaxation of atoms into these non-zero momentum states and the extent

to which the full Hamiltonian for atoms in the circularly driven lattice can be approximated by

the effective tight binding model.

The outline of this thesis is as follows. Chapter 2 starts with the theoretical background

relevant to the main experimental results in this work. Chapter 3 goes into the experimental ap-

paratus used for the work, with the primary focus on the design of our tripod-piezo mirrors that

allowed us to perform Floquet experiments. Chapters 4 and 5 go into the main experiments per-

formed, the study of parametric instabilities and the study of prethermal BEC relaxation, respec-

tively. Chapter 6 provides an outlook into future work that can be explored with our experimental

techniques.

3



Chapter 2: Background

This section provides some of the background physics that will be discussed throughout

the rest of this thesis. The primary focus is the emergence of heating and relaxation in Floquet-

engineered systems.

The starting point for all my research was a Bose Einstein Condensate (BEC) loaded into

an optical lattice. The physics of BECs and their production has been extensively studied and

documented elsewhere (see [23, 24]). I will limit the discussion of background topics to the

aspects needed for the works in this thesis

2.1 Floquet Theory

Periodically driven systems exhibit properties similar to those of systems with a periodic

spatial potential. Instead of repeating in space, the Hamiltonian repeats in time, with a period T ,

H(t+ T ) = H(t). The starting point is the Schrödinger equation:

ih̄
∂ψ(x, t)

∂t
= H(t)ψ(x, t) (2.1)

Because the Hamiltonian is time periodic, the Floquet theorem says that we can expand ψ(x, t) in

terms of Floquet modes, φn(x, t), which have the same periodicity as the Hamiltonian, φn(x, t) =

4



φn(x, t+ T ) [25]

ψ(x, t) =
∑
n

φn(x, t)e−iεnt/h̄ (2.2)

The quantities εn are known as quasi-energies, and they resemble the quasi-momentum that ap-

pear in the Bloch expansion of periodic spatial potentials. Substituting this expansion into (2.1)

gives an eigenvalue equation for the Floquet modes

(
H(t)− ih̄ ∂

∂t

)
φn(x, t) = εnφn(x, t). (2.3)

The Floquet modes are time periodic, with frequency ω = 2π/T , so they can be expanded in

harmonics of ω, namely φn(x, t) =
∑

m χnm(x)eimωt, where χnm are the spatial modes of the

mth fourier component. At this point, the total wavefunction ψ is given by

ψ(x, t) =
∑
n,m

e−iεnt/h̄eimωtχnm(x) (2.4)

The Floquet mode, φnj = φne
ijωt can be substituted into Eq. 2.3, which yields an identical

solution, but with a shifted quasienergy, namely εnj = εn + jh̄ω [26]. This implies that the

Floquet modes φn and eijωtφn are the same physical state, but they differ in quasienergy by j

photons (jh̄ω). The quasi-energy spectrum resembles a Brillouin zone, where the quasienergies

can take values in a reduced energy zone, [−h̄ω/2, h̄ω/2].

5



2.2 Time Evolution Operator

Going back to the Schrödinger equation, if we have an initial state ψ(x, t0), then we can

time-evolve this state using the operator U(t, t0), such that at time, t, (writing ψ(x, t) as |ψ(t)〉)

|ψ(t)〉 = U(t, t0)|ψ(t0)〉 (2.5)

U is the operator which is defined to be a solution to the Schrödinger equation, namely

ih̄∂tU(t, t0) = H(t)U(t, t0). (2.6)

This has the general solution

U(t, t0) = T exp

(
−i/h̄

∫ t

t0

H(t′)dt′
)
, (2.7)

where the T indicates a time ordered product. Since our Hamiltonian is periodic in time, the time

evolution operator is time periodic as well, namely evolving from a time t0 for an integer number

of periods, nT , results in U(nT, 0) = U(T, 0)U(nT − T, 0)... = (U(T ))n (where t0 is taken to

be 0 without loss of generality).

According to Floquet theory, we can now express this time evolution as

U(t2, t1) ≡ e−iK(t2)e−i(t2−t1)HeffeiK(t1), (2.8)

whereK(t) is a “Kick” operator, which is periodic with zero average value andHeff is an effective

6



Hamiltonian which, for properly chosen K(t), has no time dependence [6]. The kick operator

can also be interpreted as an effective initial condition for the drive. The effective Hamiltonian,

Heff can display properties different from the time dependent H , such as topological structure

or interesting energy dispersions. This is the premise of Floquet engineering, in that interesting

effective Hamiltonians can be simulated through a periodically driven system.

2.3 High Frequency Expansion

In most systems, Eq. 2.8 will not be useful for finding analytical results for Heff. In many

cases, a high frequency expansion in 1/ω = T/2π, can be used to accurately describe the sys-

tem. Ignoring the kick operators, the time evolution is governed by U(T ) = exp(−iTHeff) =

exp(−2πiHeff/ω). In realistic experimental setups, a smooth turning on of the drive can largely

mediate the effect of the kick operator. The time evolution operator can be expanded in terms of

the inverse frequency,

Heff =
∑
n=0

H [n]

(h̄ω)n
. (2.9)

Since the time dependent Hamiltonian is periodic, it can be expanded in fourier terms, H(t) =∑
` h`e

i`ωt. The first two terms in this high frequency expansion of the effective Hamiltonian are

given as [27, 28]

H [0] =
1

T

∫ T

0

H(t)dt = h0 (2.10)

H [1] =
∑
`

1

`
[h`, h−`], (2.11)

7



which is referred to either as the High Frequency expansion (HFE) or the Magnus Expansion [29]

and can be extended to higher orders in (1/ω).

2.4 Mean Field Theory

Having done a quick overview of Floquet theory, I now discuss the relevant physics of

cold atoms. In the second quantization description of our system of particles, the field operators,

Ψ̂(x)(Ψ̂†(x)) are defined to be the annihilation (creation) of a boson at position x. The Hamilto-

nian for N bosons in a confining potential, Vext, subject to two body interactions, V (xi − xj), is

given as

H =

∫
d3x Ψ̂†(x)

(
−h̄2

2m
∇2 + Vext(x)

)
Ψ̂(x)

+
1

2

∫
d3x d3x′ Ψ̂†(x)Ψ̂†(x′)V (x− x′)Ψ̂(x)Ψ̂(x′)

(2.12)

Properties of this Hamiltonian can be computed from this expression, but as the number of

particles increases, the complexity of the Hamiltonian increases, and calculations become com-

putationally expensive. In order to formulate a framework that can simplify these calculations,

we introduce a mean-field approach, where we separate out the contribution from the condensate

to the field operator, and treat interactions with the condensate with an average mean field. The

field operator can be expanded into a mode basis as Ψ̂ =
∑

n ψnan, which annihilates a single

particle wave function, ψn and the index n runs over a complete set of orthonormal basis states.

The mode annihilation and creation operators satisfy the usual bosonic commutation relations:

[an, a
†
m] = δn,m.

The condensate can be described as a macroscopic occupation of a single particle state and

8



adding a particle will not significantly change the physical configuration, meaning that N0 ≈

N0 + 1. In this situation, the annihilation and creation operators can be expressed as just a

complex number, which we denote as a = a† =
√
N0, such that N = a†a. The BEC mean field

single-particle state, ψ0 =
√
N0/V , confined to a volume V , has zero momentum and its density

is given by |ψ0|2, which is proportional to the total number of condensed particles, N0. The total

state in mean this mean field approach is then represented as the BEC plus a small fluctuation,

Ψ = ψ0 + δψ.

2.5 Weakly Interacting Condensates : GP Limit

Given the classical field Ψ(r, t), we first assume that there are no fluctuations and that the

single particle state is time dependent, Ψ(r, t) = φ(r, t) and look at the time evolution from the

non-linear Schrödinger equation:

ih̄∂tφ(r, t) = − h̄
2∇2

2m
φ(r, t) + Vext(r) φ(r, t)

+

(∫
d3x′φ†(r, t)V (x′ − x)φ(r, t)

)
φ(r, t).

(2.13)

For weakly interacting condensates in a dilute, low temperature system, the interaction is dom-

inated by s-wave scattering, so we replace the interaction with a delta function with strength

proportional to the scattering length, a:

V (x′ − x) =
4πh̄2a

m
δ(x′ − x). (2.14)

9



Substitution of this into (2.13) results in the Gross-Pitaevskii (GP) equation:

ih̄∂tφ(r, t) =

(
− h̄

2∇2

2m
+ Vext(r) +

4πh̄2a

m
|φ(r, t)|2

)
φ(r, t). (2.15)

This nonlinear equation governs the evolution for the condensate under the assumption that there

is a macroscopic number of particles, the scattering length is smaller than the average distance

between bosons, and the interactions are weak enough that the ground state is not highly corre-

lated.

The ground state of the mean-field Hamiltonian can be written φ(r, t) = Φ(r)e−iµt/h̄ where

Φ is a real function with normalization given by
∫
dxΦ(x) = N and µ is the chemical potential.

Defining the interaction strength to be g = 4πh̄2a/m, the GP equation for the ground state

simplifies to the following nonlinear Schrödinger equation:

µΦ(x) =

(
−h̄2∇2

2m
+ Vext(x) + gΦ2(x)

)
Φ(x). (2.16)

In the limit where the kinetic energy can be neglected with respect to interactions, known

at the Thomas-Fermi regime, the atomic density, n(x) = Φ(x)2 can be written as

n(x) =
1

g
(µ− Vext(x)) . (2.17)

For a harmonic external potential (such as an optical dipole trap), this is an inverted parabola

for points where the chemical potential is greater than the external potential, and 0 otherwise

(since Φ(x) is real). In realistic experimental setups, there is a thermal component as well, so the

total atomic distribution can be approximated as a sum of a thermal gaussian distribution and this

10



Thomas Fermi inverted parabola.

2.6 Bogoliubov Theory

I now consider excitations of the BEC, which is assumed to be homogenous and lattice

free, in the presence of weak interactions. To characterize interactions, we choose the plane-

wave basis for expanding the field operator, Ψ̂†(x) = âpe
ipx/
√
V , where V is the volume of the

uniform gas, which results in the following momentum space Hamiltonian:

H =
∑
p

p2

2m
a†pap +

g

2V

∑
p1,p2,q

a†p1+qa
†
p2−qap1ap2 . (2.18)

In the limit of no interactions, the free particle dispersion, p2/2m, is recovered. To include

the interaction terms, we can apply perturbation theory, starting with atoms condensed in the

BEC (p = 0) state. Allowing for excitations into a plane wave of momentum k, the total number

of atoms will be given as N = N0 +
∑
a†kak ≈ N0. If we expand the Hamiltonian in powers

of a(†)
k and ignore any higher order excitations, everything remains condensed in p = 0, and the

contribution is g
2V
a†0a

†
0a0a0 =

gN2
0

2V
.

First order processes in the Hamiltonian violate conservation of momentum, as they contain

terms such as a†ka
†
0a0a0 = N0a

†
ka0. In order to conserve momentum, terms that include second

order processes have to be included, which results in:

g

2V

∑
p1,p2,q

a†p1+qa
†
p2−qap1ap2 →

gN0

2V

∑
k

(
aka−k + a†ka

†
−k + 4a†kak

)
, (2.19)

where these come from the following combinations of (p1, p2, q) : (0,0,k), (0,k,0), (k,0,0),

11



(0,k,k), (k,0,−k), (k,−k,k). All of these terms result in momentum conservation and are pro-

portional to N0. In this Bogoliubov approximation, we can make the assumptions that N0 ≈ N

and N2
0 ≈ N2 − 2N

∑
k a
†
kak and rewrite the Hamiltonian as:

H =
gN2

2V
+
∑
k

(
k2

2m
+
gN

V

)
a†kak +

gN

2V
(a†ka

†
−k + aka−k), (2.20)

which is not diagonal in the ak basis. The linear transformation that diagonalizes Eq. 2.20 is

given by [30]:

ak = ukbk + vkb
†
−k

a†k = ukb
†
k + vkb−k

(2.21)

To maintain the same commutation relations, [bk, b
†
k′ ] = δkk′ , the relationship |uk|2 − |vk|2 = 1

must hold, which is done by defining a parameter θk, in which uk = cosh(θk), vk = sinh(θk).

Diagonalizing Eq. 2.20 requires that the b†kb
†
−k + bkb−k terms vanish, and that the following

relationship holds. (
k2

2m
+
gN

V

)
ukvk +

gN

2V
(u2

k + v2
k) = 0 (2.22)

Substitution of Eq. 2.21 into Eq.2.22 defines our parameter θk:

coth(2θk) = −
k2

2m
+ gN

V
gN
V

(2.23)

The diagonal contribution of 2.20 is given by

∑
k

[(
k2

2m
+
gN

V

)
(u2

k + v2
k) +

gN

V
ukvk

]
b†kbk =

∑
k

E(k)b†kbk, (2.24)

12



where

E(k) =

√(
k2

2m
+
gN

V

)2

−
(
gN

V

)2

(2.25)

The collective mode dispersion, Eq. 2.25, is shown in Fig. 2.1. In the limit of no interac-

tions, the quadratic free particle dispersion is recovered and in the presence of interactions, for

small momentum, the dispersion is linear, with a speed of sound given by vs =
√
gN/V m [31].

As the interaction parameter is increased, the speed of sound, or the slope of the dispersion, is

increased as shown in the figure. The appearance of this linear dispersion can be interpreted as

a Goldstone mode arising from the broken U(1) symmetry of the BEC [32]. In the limit of large

momentum (low wavelength), the parabolic free particle dispersion is recovered, with an offset

equal to the interaction parameter, U .

4 2 0 2 4
k/2m

E(
k)

U=0
U=5
U=10
U=50

50 0 50
k/2m

E(
k) 49 50

a) b)

U

Figure 2.1: Bogoliubov mode dispersion for different interaction parameters, U = gN/V . (a)
Low momentum, where the dispersion becomes linear with increasing interactions due to the
Goldstone mode appearing from the broken U(1) symmetry. (b) High momentum, where the
behavior is quadratic, just offset by the interaction parameter U, which is shown in the inset.

13



2.7 Optical Lattices

Optical lattices are periodic structures that arise from the interference between electric

fields of propagating light. For an arbitrary number of beams, the intensity is given as the absolute

square of the electric fields, i.e. I(r, t) = |
∑N

i
~Ei(r, t)|2. In the simplest case of a 1D lattice,

where a beam is reflected on itself, the total Electric field is given as E = E0(eikx + e−ikx) =

2E0 cos(kx). The intensity is just the magnitude squared of this, I = 4|E0|2 cos2(kx) = I0(1 +

cos(2kx)), which is just a standing wave with a spatial period of half the wavelength.

Optical lattices are ideal for quantum simulation tasks because different geometries can be

formed by engineering different electric fields. Multiple lattices can be superimposed, resulting

in so-called superlattices, such as checkerboard, kagome or honeycomb lattices. Examples of

quantum simulation can include studies of magnetization effects, or band structure studies, which

may exhibit rich topological features.

For neutral atoms subject to a standing optical wave, their dynamics are governed by the

optical dipole force. If we consider our atom to be a two level system, with energy states denoted

|e〉 and |g〉 for the excited and ground state, respectively, then the energy shift, to second order in

the dipole interaction, Hint = −d · E, is given by [33]:

δε
(2)
i =

∑
i 6=j

|〈j|dE|i〉
εi − εj

. (2.26)

Defining the Rabi frequency as h̄Ω ≡ d·E which couples the two states, and the energy difference

14



as εi − εj ≡ h̄∆, this can be written as:

V = δε =
Ω2

∆
. (2.27)

There are a few important things to note. Since Ω is proportional to the electric field, the potential

felt by the atoms will be proportional to the intensity, which means that for a periodic wave

(lattice), the atoms will see a potential of the same period. Another thing to consider is the

Rabi frequency. The orientation of the dipole matrix operator was ignored in this simplified

derivation, but in the full picture, this contains information regarding the strengths of certain

transitions. They are effected by Clebsch-Gordon coefficients, so the same intensity of light can

couple different states differently, or not at all. Transitions are also dependent on the orientation

of the Electric fields, as angular momentum must be conserved, meaning certain transitions may

require circular polarization, whereas others require a linear polarization. The final thing to

consider is that atoms in the lab are not two level systems (as I have been told, ”There are no two

level atoms, and Rubidium is not one of them”). In order to account for the atomic structure, the

sum in Eq. 2.26 should be carried out for all dipole allowed transitions.

2.7.1 Band Structure

Now that we have a form for the potential felt by neutral atoms, we wish to solve the

Schrödinger equation for atoms in an optical lattice:

− h̄2p2

2m
ψ + V ψ = Eψ (2.28)

15



Taking the 1D case as an example, we substitute in the lattice potential, V = V0 sin2(2kx), where

V0 will be a function of the atomic dipole matrix elements and the detuning, and k ≡ 2π/λ:

− h̄2p2

2m
ψ + V0 sin2(kx)ψ = Eψ (2.29)

Given the periodicity of the lattice potential, we can employ the Bloch theorem, which

states the solutions of this equation will be the product of a plane wave, exp(iqx), and a func-

tion unq (x), which will have the same periodic dependence as the lattice. The n subscript now

represents the band index and q represents the quasi-momentum which is defined in the range

[−π/λ, π/λ]. We now rewrite the Schrödinger equation in terms of this ansatz (working with

h̄ = 1):

Hφnq = En
q φ

n
q(

− p2

2m
− V0

4

(
e2ikx + e−2ikx − 2

))
eiqxunq (x) = En

q e
iqxunq (x) (2.30)

Due to the periodicity, it is helpful to Fourier transform un in terms of the lattice momentum

harmonics

unq =
∑
`

ale
2i`kx (2.31)

The kinetic term in the Hamiltonian will just shift p to p + q, where p will be the momentum of

the `-th component and the lattice potential shifts the Fourier component by±2k. This will result

in the following eigenvalue equation

16



∑
`

(2k`+ q)2

2m
a`e

2ik` −
∑
`

V0

4
a`
(
e2ik(`+1) + e2ik(`−1)

)
= En

q

∑
`

a`e
2ik` (2.32)

We can match plane wave solutions, which results in the following eigenvalue problem for the

coefficients

a`
(2k`+ q)2

2m
− V0

4
(a`−1 + a`+1) = En

q al (2.33)

Defining recoil units to be qr = h̄kr, ER = h̄2k2r
2m

, with kr = 2π/λ we can recast this as:

(2`+ q/qr)
2a` −

V0

4ER
(a`−1 + a`+1) =

En
q

ER
al (2.34)

In order to numerically solve this, we truncate the sum to a finite number of components, `max,

and compute the plane wave states for ` = (−`max, `max) which is an eigenvalue problem, where

the Hamiltonian is just a tri-driagonal matrix with the following structure.

17





. . . . . . 0 0

. . . (2 ∗ (`− 1) + q/qr)
2 −V0/4 0

0 −V0/4 (2 ∗ `+ q/qr)
2 . . .

0
. . . . . . (2 ∗ (`+ 1) + q/qr)

2





a−`max

...

a`−1

a`

a`+1

...

a`max



= E



a−`max

...

a`−1

a`

a`+1

...

a`max



(2.35)

Wannier Functions

In the previous section, the energy band structure was solved for in the presence of an op-

tical lattice. The next step is to determine what the eigenvectors, or the eigenstates are for an

atom in an optical lattice. The matrix equation, 2.35, solves for the coefficients introduced in the

bloch function (Eq. 2.31), and the bloch wave associated with this solution is delocalized over all

sites, due to the a†iai+1 terms, which are off-diagonal in the site location. In order to determine

18



a localized wave function for an atom in a given site, we introduce the Wannier functions, which

make up an orthonormal basis for an atom in the nth energy band. The Wannier functions repre-

sent the most localized function one can construct from the Bloch solutions confined to a single

band, and can be expressed as:

Wn(x− xi) =
∑
q

e−iqxiφnq =
∑
`

∑
q

a` e
−iqxie2ikx` (2.36)

Figure 2.2 shows the lowest Wannier orbital for increasing lattice depths. As the lattice depth

increases, the Wannier function becomes more localized to a specific site, and there is a decreased

lobe on the next nearest site, indicating a reduced probability of tunneling to the next site. We

note that the Wannier functions are not eigenstates of the Hamiltonian.

2.0 1.5 1.0 0.5 0.0 0.5 1.0 1.5 2.0
Lattice Position

|w
0(

x)
|

V0 = 0
V0 = 5ER

V0 = 30ER

Lattice

Figure 2.2: Wannier function corresponding to ground band of optical lattice at different depths.
The absolute value is taken to show the overlap on the next nearest site.

19



2.7.2 2D Checkerboard Band Structure

In the Rubidium-1 experiment, we work with a 2D optical lattice that is formed by a single

beam that is folded to retro-reflect upon itself twice (sketch in Chapter 3). The relative po-

larization allows for control of an in plane and out of plane lattice, resulting in a dynamically

controllable lattice with 2 sites per unit cell. The total potential seen by the atoms is given by

([34]):

I = I2D + Iz = 2 cos(2kx− 2θxy − 2φxy) + 2 cos(2ky + 2φxy) + 4︸ ︷︷ ︸
I2D

+

4 (cos(kx− θz − φz) + cos(ky + φz))
2︸ ︷︷ ︸

Iz

,

(2.37)

where the φ and θ variables define the relative phase due to the polarization.

For the case of θxy = θz = φxy = φz = 0, we have a pure checkerboard lattice where

the out of plane lattice has twice the spatial period and the double well structure is seen. Figure

2.3 shows an example lattice where the double well structure appears due to the sum of the two

lattices. The phase θ controls the location of the λ (out of plane) lattice, and Fig. 2.3 shows

that this allows for control of the geometry from something that has double wells to a dimerized

lattice.

Just as with the 1D case, we now seek solutions of the Schrödinger equation Hψ = Eψ

and start by expanding in plane waves ψ =
∑
eiq·rφq(r). In the checkerboard case, the lattice

can be written in terms of the reciprocal lattice vectors, b̂1(2) = x̂ ± ŷ. The in plane and out of

20



In Plane Lattice

Out of Plane Lattice

= 0
= /2

Ixy + . 25Iz
 =0
= /2

Figure 2.3: Checkerboard optical lattice, where the white (black) indicates maximum (minimum)
intensities. The total intensity (bottom left) is the sum of the in-plane (top) and out of plane
(middle) intensities. The red curves on the right indicate a different θ, changing the location of
the out of plane lattice from the λ/2 minimum to the barrier (bottom right is offset for visual aid).

21



plane potentials can be rewritten as:

Vxy(x, y) = −V2D

8

(
ei(b1+b2)·r + e−i(b1+b2)·r + ei(b2−b1)·r + e−i(b2−b1)·r + 4

)
(2.38)

Vz = −Vz
16

(4 + 2eib2·r + 2e−ib2·r + 2eib1·r + 2e−ib1·r + ei(b1+b2)·r + e−i(b1+b2)·r

+ei(b2−b1)·r + e−i(b2−b1)·r)

(2.39)

Expanding the bloch function in terms of plane wave components φn =
∑

m,n amne
imb1+inb2 ,

the 2D potential couples terms that differ in reciprocal vectors by either 1 b1(2) or any combination

of ±b1 ± b2. An example Hamiltonian is illustrated in Fig. 2.4 in the form of a matrix plot, illus-

trating the sparseness. This example is using a cutoff of 2 plane waves, resulting in a (2n + 1)2

= 25 x 25 matrix, where each element denotes a plane wave with momentum nb1 + mb2. The

diagonal elements represent the kinetic energy terms, (p+ 2nk + 2mk)2.

22



0 5 10 15 20
0

5

10

15

20

Figure 2.4: Matrix plot of the elements of a Hamiltonian matrix for a cutoff of 2 plane waves.
The red points on the diagonals represent the kinetic energy, whereas the off-diagonal elements
represent coupling due to the lattice.

Figure 2.5 shows the two lowest bands for a 3 ER lattice with 0.5 ER out of plane lattice.

Under these conditions, the q = 0 gap is approximately 3600 Hz, and the total width of the ground

band is 1100 Hz (defined as max-min). Figure 2.6 shows the effect of increasing the out of plane

lattice while keeping the same in plane lattice. Increasing the tilt increases the gap along the

23



edges while simultaneously flattening the lowest bands.

0

1

2

3

E (kHz)
M
X

0

150

300

450

600

750

900

1050

Figure 2.5: Left: Lowest two bands for a 3 ER lattice with 0.5 ER out of plane lattice. Right:
projection of lowest band onto 2D plane. Γ, K and X points are labeled.

24



M X

E

VZ = 0 ER

VZ = 0.5ER

Figure 2.6: Lowest 4 bands for a 3Er lattice for the case of 0 Er out of plane lattice (black) and a
0.5 Er out of plane lattice (red).

2.8 Bose-Hubbard Model

We now revisit the Hamiltonian for atoms in an optical lattice. Taking the 1D situation

with delta function interactions that are dominated by s-wave scattering, and an external lattice

potential, Eq. 2.13 can be written as

H =

∫
dx Ψ̂†(x)

(
−h̄2

2m
∂2
x + V0 sin2(kLx)

)
Ψ̂(x)

+

∫
dx

g

2
Ψ̂†(x)Ψ̂†(x))Ψ̂(x)Ψ̂(x)

(2.40)

25



If the lattice potential is strong enough, we can make the tight binding approximation, where

the field operators are expanded in a Wannier orbital for each band : Ψ̂†(x) =
∑

iwi(x)b†i .

Substituting into the above Hamiltonian results in the Bose-Hubbard Hamiltonian [35]:

H = −J
∑
<ij>

b†ibj +
U

2

∑
i

ni(ni − 1) (2.41)

The first term allows for hopping/tunneling of an atom from a site i to a site j. The < i, j >

notation indicates that the model accounts only for nearest neighbor tunneling, which has terms

for an atom going from site i to i ± 1. The second term is the interaction energy, which is the

energy associated with having multiple atoms per site. The tunneling matrix element, J , is given

as

J =

∫
dx w(x− xj)

(
h̄2

2m

d2

dx2
+ V (x)

)
w(x− xi) (2.42)

In the tight binding approximation, we can approximate the sinusoidal potential with a harmonic

oscillator, and the on site wannier functions will tend towards the Hermite-Gauss solutions for

the harmonic oscillator. Plugging in those basis states into the above equation, which converges

at deeper lattice depths, the tunneling can be expressed as

J =
4√
π
ER

(
V0

ER

)3/4

exp

(
−2

√
V0

ER

)
(2.43)

This was seen in the previous sections on Optical Lattices, where the numerical band struc-

tures flattened out with a higher potential. In this model, a lower J corresponds to less energy

needed to tunnel between sites and the wave function delocalizes over many lattice sites, and

is known as the superfluid phase. In the limit that J gets large, the energy required to tunnel

26



becomes too large, so the atoms will localize to individual sites which is known as the Mott

Insulating phase [36].

We can additionally calculate the on-site interaction as

U =
4πh̄2as
m

∫
dx|wn(x)|4 ≈

(
V0

Er

)3/4

. (2.44)

The ratio, J/U , can be controlled by changing the lattice depth, and can be used to observe

the superfluid to mott-insulator phase transition [37, 38]. In optical lattices, the superfluid to

mott-insulator (MI) transition can be seen through the coherence of time of flight images. In the

superfluid phase, there is phase coherence in the system and snapping off the lattice will project

the atoms onto plane wave components and diffraction peaks are seen. When the MI phase is

approached, phase coherence is lost across sites and in the time of flight image, we just see a

gaussian thermal cloud after time of flight. An example in our lab is shown below

27



0.0 ER 4.2 ER 8.4 ER

12.6 ER 16.8 ER 21.0 ER

Figure 2.7: Time of flight images for atoms loaded into a variable lattice depth, held for 5ms
then snapped off. As the depth is increased, we see an initial diffraction pattern due to the
phase coherence. As we go past the transition, we note that the phase coherence is lost and
the mott-insulator state is obtained. The presence of the additional diffraction peaks are due to
experimental imperfections in the checkerboard lattice alignment.

2.8.1 Checkerboard Tight Binding

The Bose-Hubbard model from the previous section can be further extended to the 2D

checkerboard lattice when we restrict ourselves to the lowest two bands of the system. The

checkerboard lattice consists of double wells, which are offset in energy by ∆ and there are two

sites per unit cell, which we denote the a, b sites. Following the same convention as the Bose-

Hubbard model (assuming non-interacting, U = 0) we can write the first order tight binding limit

28



as

H = −J
∑
i

a†i (bi+1 + bi+2 + bi+3 + bi+4) + ∆
∑
i

a†iai. (2.45)

To calculate the band structure, we perform a Fourier transform, ai =
∑

q aqe
iqx and the

Hamiltonian is:

H = ∆
∑
q

a†qaq − J
∑
q

a†qbq(4 cos(qx/2) cos(qy/2)) + H.C., (2.46)

This is a 2x2 matrix in the a, b basis, which can be diagonalized to give the tight binding energy

bands:

E =
∆

2
±
√

∆2 + 64J2 cos2(qx/2) cos2(qy/2)

2
. (2.47)

The tight binding model can be interpreted as a fit to the exact band structure. To first order,

we can see that the ground band does not capture the curvature along the band edge, as shown in

Fig. 2.8. As a reference, the exact calculation in the figure is calculated with a total lattice depth

of 5ER, where 0.5ER (≈ 1.7 kHz) is out of plane, which results in a q = 0 gap splitting of 3.2

kHz and an approximate value of J being 80 Hz (taken as 1/4 the total width of the ground band).

The resulting best fit for the tight binding gives J ≈ 210 Hz and ∆ ≈ 2.4 kHz.

29



M X

0

50

100

150

200

250

300

350
E(

Hz
)

Exact
First order
Second order

Figure 2.8: Comparing Tight Binding to the full plane wave band structure using only nearest
neighbor tunneling. The first order tight binding model appears flat along the band edge and
does not completely capture the band, whereas the inclusion of next nearest neighbor tunneling
captures the band dynamics very well.

One of the key features of the first order tight binding is that the bands are flat along the

edge of the Brillouin zone. In order to see a more qualitative match, we can add in the next

nearest neighbor tunneling. This Ja(Jb) term allows for tunneling between an a (b) site and any

of the 4 next nearest neighbors, e.g. HNNN = Ja
∑4

j=1 a
†
iaj . After Fourier transforming the

Hamiltonian into momentum space, we end up with a 2x2 matrix for each q:

H2 =

∆− 4Ja cos((qx + qy)/2) cos((qx − qy)/2) −4J cos(qx/2) cos(qy/2)

−4J cos(qx/2) cos(qy/2) −4Jb cos((qx + qy)/2) cos((qx − qy)/2)


(2.48)

30



The resulting band structure has fourier components along the diagonal axes (qx ± qy), which

results in the curvature along the band edge that is seen through the exact diagonalization in the

plane wave basis. Figure 2.8 shows that even including the second term is enough to approximate

the system.

2.9 Artificial Magnetic Fields on a Lattice

The Hubbard model was initially introduced to describe electrons hopping on a lattice. If

we treat our Bose-Hubbard system as a charged particle and add in a magnetic field, which is

given by the curl of a vector potential, B = ∇× A, then when the charged particle tunnels from

one site to another, it will acquire an Aharanov-Bohm type phase given by φim,n = −eAim,n. In

the tight binding language, we refer to this phase as the Peierls phase, [39]. The non-interacting

tight-binding Hamiltonian can then be expressed as

H = −J
∑
n,m

(
eiφ

x
n,ma†n+1,man,m + eiφ

y
n,ma†n,m+1an,m + H.C

)
(2.49)

A particle that tunnels along a closed loop of one lattice plaquette will acquire a net phase of

Φ = φxn,m + φyn,m+1 − φxn+1,m − φyn,m = 2πα, (2.50)

where the - sign comes from tunneling in the conjugate direction. This phase, Φ, is the flux per

unit cell and in the electron/lattice case, the quantity α represents the flux in units of the magnetic

flux quanta, h/e.

This phase is assumed to come from the charged particle, but if we extend this to a neutral

31



particle, then there would be no interaction with the vector potential. The phenomena of an

artificial magnetic fields arise when we artificially generate complex tunneling coefficients in

the Hamiltonian, as has been demonstrated using raman transitions, or in the case of this thesis,

lattice shaking. These techniques result in an atom that tunnels between sites acquiring a phase,

so while there is no charge, it acts as if it were a charged particle moving in a magnetic field.

2.9.1 Harper-Hofstadter Hamiltonian

Since we are dealing with an equivalent vector potential, we are, in principle granted the

freedom of a gauge choice. The observables, such as energy spectrum and average momentum,

are gauge independent, so we can select what suits the system. A typical choice is the Landau

gauge, where we only allow for complex tunneling in one direction, namely φn,m = (−Φn, 0),

which has a magnetic field acting along z. In this gauge, the Hamiltonian is

H − J
∑
n,m

e−iΦna†n+1,man,m + a†n,m+1an,m + H.C. (2.51)

The single particle energy spectrum of this model as a function of the flux, α, results in the

Hofstadter butterfly [22]. When the flux is a rational fraction, α = p/q, the resulting energy

structure results in q different subbands, and there is a magnetic brillouin zone, which is q times

larger than the lattice cell.

For our artificial field lattice system, the fact that we can turn on and off the gauge field

means that the system is in fact not gauge independent. Projective measurements that arise when

snapping off the modulation generating the gauge field, depend on the effective gauge of the

experiment.

32



2.10 Periodically Driven Optical Lattices

In order to simulate Floquet physics in an optical lattice, there are various techniques used,

including lattice shaking and lattice amplitude modulation. Lattice shaking is performed by pe-

riodically displacing the lattice, either through a frequency shift in the lattice light or through

a physical displacement of a mirror (typically using a piezoelectric device), whereas amplitude

modulation is just changing the lattice power using AOMs or other electronics. While the two

techniques result in different effective Hamiltonians, the coupling between lattice bands differs

in each scenario.

In these sections, I describe some theoretical treatments of Floquet problems, which are

highlighted through different example systems.

2.10.1 Tight Binding Floquet Treatment

As previously discussed, cold atoms in an optical lattice can be treated with the use of a

tight-binding model, where we start with a kinetic Hamiltonian: H = −J
∑
〈i,j〉 a

†
iaj . In the

tight-binding Floquet treatment, we start with this kinetic Hamiltonian and add in an additional,

time-periodic term. We then find a unitary transformation that can be applied, which rotates out

the time dependent component.

Time Dependent Force

Before considering the Floquet treatment of performing a unitary transformation, we start

with a simple, intuitive picture in this tight-binding model. In the presence of a time dependent

force, F = K sin(ωt), q changes based on the semi-classical equation of motion, q̇ = F . At a

given time, the quasimomentum can be expressed as q = q0 +F sin(ωt)/ω [40]. If the frequency

33



is fast enough that we can approximate the dynamics with the time averaged dynamics, and

recalling that the tight-binding dispersion was given by E = −2J cos(q), we obtain an effective

Hamiltonian

Heff = −2J
1

T

∫ T

0

cos(q0 +
K

w
sin(ωt))dt = −2JJ0(K0) cos(q), (2.52)

where Jn is the nth order Bessel function, and a scaled driving strength, K0 = K/ω has been

introduced. This scaled hopping amplitude is typically referred to as the effective tunneling,

Jeff ≡ JJ0(K0). The effective band structure is plotted for different drive strengths in Fig. 2.9.

When the drive strength reaches 2.4, the band flattens as a result of the Bessel function being 0,

and beyond that point, J0 is negative, so the band inverts, making q = 0 the maximum point.

The fact that the effective time averaged dispersion can become flat is known as dynamic

localization (coherent destruction of tunneling) and has been shown to occur in cold atom systems

[41, 42, 43]. Due to the lower effective J , this Floquet shaking has been used to realize the Mott

insulator to superfluid transition [44].

Tilt Modulation

The next case to consider is a double well system, where there are two sites, A and B,

with tunneling J and an offset ∆ (taken to be on B). We take ∆ >> J and assume there is no

tunneling between A and B sites. The Hamiltonian is

H0 = −J
∑
i

a†ibi + b†iai + ∆
∑
i

b†ibi, (2.53)

34



0 2 4 6
K0

0.4

0.2

0.0

0.2

0.4

0.6

0.8

1.0

0

3 2 1 0 1 2 3
q

E(
q)

K0=0
K0=1.2
K0=2.4
K0=3.6

Figure 2.9: Left: J0(K0). Right: Modification of the effective band structure for different K0.
For certain ranges of the drive strength, the band inverts and the minima of the band structure is
at the edge of the band.

35



which can be rewritten in matrix form as

H0 =
∑
i

(a†ib
†
i )

 0 −J

−J ∆


ai
bi

 (2.54)

Diagonalizing the system gives 2 energies separated by

δE = E+ − E− =
√

4J2 + ∆2. (2.55)

We now allow for modulation of A sites, such that the Hamiltonian is time dependent, with

frequency ω,

H(t) = H0 + α cos(ωt)a†iai. (2.56)

To analyze this system, we start by moving into a rotated frame to eliminate the total

potential (V (t) = ∆ b†b+ α cos(ωt)a†a), which can be done by using the following rotation:

R = exp

(
iα

ω
sin(ωt)a†iai + i∆tb†ibi

)
(2.57)

Under this transformation, the new Hamiltonian is given by, Ht = RHR† − iR∂tR† which will

become

Ht = −J exp
(
−i∆t+

α

ω
sin(ωt)

)
a†ibi + H.C. (2.58)

In the case of resonant driving, ω = ∆, using the HFE, Eq.2.10, it can be seen that the transformed

Hamiltonian only has an H [0] contribution, and the effective Hamiltonian is given by (after using

36



the Jacobi-Anger expansion for Bessel functions)

Heff = −JJ1

(α
ω

)
a†ibi + H.C. (2.59)

Just like in the previous example, there is now an effective tunneling that takes the form

of a Bessel function. When we take ∆ = nω, which is an n-photon resonance, the effective

tunneling is given by the nth order Bessel function, Jn
(
α
ω

)
. The presence of the drive increases

the effective tunneling (at low enough α/ω) and can be interpreted as allowing tunneling through

absorption of a drive photon. Note that this is the lowest order term in the expansion, which

reintroduces coupling between A and B sites. Higher orders in the expansion can be calculated,

and are actually necessary in some systems to introduce topological gaps in energy bands [45].

2.10.2 Trotter Method

The previous section introduced two interesting aspects of Floquet engineering in tight

binding systems, dynamic localization and renormalization of tunneling, which were solved an-

alytically. In this section, I will describe the Trotter Method for calculating band structures in

the Floquet picture. I apply this to an optical lattice subject to amplitude modulation, where the

intensity of light seen by the atoms is periodically modulated.

For an amplitude modulated lattice, the Hamiltonian is given by

HAM =
p2

2m
+ (1 + α sin(ωt)) VLAT = H0 +H ′(t), (2.60)

where H0 is the same static Hamiltonian as previously discussed, and there is now a time de-

37



pendent contribution H ′(t). In section 2.7.2, H0 was solved in the plane wave basis, and is the

starting point for this method. The time evolution operator can be expressed as

U(T ) = exp[− i
h̄

∫ T

0

(H0 +H ′(t′)) dt′]. (2.61)

To numerically calculate this, we break the integral into a discretized product of N steps, where

T = N(∆t), U(T ) = UNUN−1 . . . U2U1, with Ui = exp[−i(∆t)H(ti)/h̄]. The Trotter-Suzuki

expansion allows us to expand this as [46]

exp[−i(∆t)H(ti)/h̄] ≈ exp[− i

2h̄
(∆t)H ′(ti)] exp[− i

h̄
(∆t)H0]

exp[− i

2h̄
(∆t)H ′(ti)].

(2.62)

To calculate the quasi-energy spectrum, this time evolution is calculated in the plane wave

basis and then diagonalized, which can then be used to extract the quasi-energy. One thing to

note is that since all the energies overlap into a Brillouin zone structure, there are no ”lowest

energy states”, and the labeling of states is not as straightforward as they were in the static case.

In order to label, we work under the assumption that in the Trotter picture, the time dependent

contribution is sufficiently small, so that it only slightly perturbs the initial states. We can label

the Floquet bands by using the overlap of the static eigenvectors with the Floquet ones,

|ψfi 〉 = max
j
|〈ψfj |ψsi 〉|. (2.63)

Note that at sufficiently strong drives, the concept of a particular Floquet state arising from an

undriven state loses meaning as the states become heavily mixed.

38



A sample band structure plot is shown in Fig. 2.10, which shows that the two static bands

are dressed (the second band is shifted down by the drive frequency, ω), and the amplitude mod-

ulation couples the two bands, resulting in the structure shown in 2.10(b). This figure represents

a cut along qy = 0, and it is interesting to note that due to the symmetry of the problem, the full

2D band structure exhibits a ring structure, which is of interest in condensed matter fields [47].

1 0 1
qx

0

250

500

750

1000

1250

1500

1750

2000

E i

(a)

1 0 1
qx

0

250

500

750

1000

1250

1500

1750

2000
i

(b)

1 0 1
qx

0

500

1000

1500

2000

2500

i

(c)

Figure 2.10: Band structure from trotter calculation for α = 0.12, ω = 2800 Hz. (a) shows the
static bands, with the first excited band (blue) offset by ω. (b) shows the Floquet band structure
for stated with maximum overlap with the two lowest bands in the system in the weakly coupled
limit. (c) shows the full quasi-energy spectrum for all diagonalized basis states. The lowest two
bands are indicated.

2.10.3 Extended Basis

I now present another technique for calculating the Floquet band structures, called the

extended basis method. In order to describe this theoretical treatment, I apply it to the case

39



of a periodically displaced (shaken) lattice. I note that this problem can be solved using the

Trotter method as well, but at higher drive amplitudes, the overlap and sorting of bands becomes

problematic and computationally expensive.

For this example, the Hamiltonian can be written as

HS =
p2

2m
+ VLAT(x+ δx sin(ωt), y + δy sin(ωt + φ), t), (2.64)

where VLAT(x, y) is the time dependent lattice potential.

We now introduce the co-moving frame through the unitary transformation U = exp[−iδ ·

p/h̄]. Treating ψ0 as the untransformed wavefunction, the following unitary transformation

moves us into a co-moving frame, ψcm = Uψ0. Plugging into the Schrödinger equation iψ̇ = Hψ

one obtains the following Hamiltonian for the co-moving frame

Hcm =
p2

2m
+ VLAT (x, y)− ω (δx sin(ωt)px + δy sin(ωt+ φ)py) (2.65)

We begin by recalling the eigenvalue problem for Floquet modes, Eq. 2.3 and define the

quasi-momentum operator to be

Q̂ ≡ Ĥ(t)− ih̄∂t, (2.66)

such that Q̂|φnm〉 = εnm|φnm〉, where the states and energies now are defined with two indices,

n originating from the initial Hamiltonian, and m, related to the basis of time dependent periodic

functions. In general, if |un〉 are solutions of the time independent Hamiltonian (such as a Bloch

wave in the static lattice), then the form of the combined state in the extended basis can be

40



expressed as:

|φn,m〉 = eimωt|un〉. (2.67)

By denoting the |φn.m〉, as |n,m〉, the matrix elements of this operator can be expressed as

〈n′m′|Q|nm〉 = 〈n′m′|H|n,m〉+ h̄mωtδnn′δmm′ , (2.68)

where the δ functions arise due to the orthogonality of the states for n and m. The first matrix

element, 〈n′m′|H|nm〉 can be further decomposed

〈n′m′|H|nm〉 =

∫ T

0

dt〈n′|H|n〉ei(m′−m)ωt = 〈n′|Hm′−m|n〉. (2.69)

The expression Hm′−m represents the fourier components of the co-moving Hamiltonian (2.65).

Intuitively, it can be thought of as describing an ‘m photon’ transition, where we couple two n

states that differ in energy by mω. Since our shaken Hamiltonian has a DC (ω = 0) component,

as well as a component at ω, the calculation will involve just two fourier components, and the

quasi-momentum operator can be represented as a tri-diagonal block matrix.

Q =


H0 − h̄ω ⊗ 1̂ H−1 0

H1 H0 H−1

0 H1 H0 + h̄ω ⊗ 1̂

 . (2.70)

The calculation of this matrix requires truncation of the time basis (how many m states),

and there are some edge effects related to the truncation. For our parameter region, we tend to

find that using mmax = 5 (11 total time “blocks”) results in convergence of the eigenvalues. The

41



parameter that defines the coupling between the bands, α is given as the Fourier component of

the co-moving Hamiltonian (2.65), which is the same as the dimensionless parameter K0 from

the tight binding picture. Figure 2.11 shows an example of the extended basis calculation in

the case of a frequency that is less than the resonance condition between the two bands. In this

calculation, the number of time basis states was taken to be 2, such that the calculation includes

[−2ω,−ω, 0, ω, 2ω]. As shown in the figure, the bands repeat themselves, which can get intricate

as the number of time states is increased. In order to isolate the relevant bands, the ground band

is taken to be the one closest to the zero energy (where there is no frequency offsets and no

coupling included). We note that there are additional levels repelling each other, even with states

differing by more than ω which indicates that a tight binding model at high drive will not capture

all relevant levels and a more complete description is necessary to understand the entire system.

42



1 0 1
qx

E0 2

E0

E1 2

E0

E1

E0 +
E1

E0 + 2
E1 +

1 0 1
qx

Figure 2.11: Sample bands for extended matrix calculation in the case of no coupling (left) and
a coupling strength, α=2.2 (right), using a truncation of 2 time states. This shows the general
repeating structure and how the coupling can act to couple different bands.

43



Chapter 3: Experiment

In this section, I will discuss the experimental apparatus designed for our optical lattice ex-

periments. We employ standard magnetic and optical trapping techniques to trap an atomic cloud

at sub-doppler temperatures. Bose-Einstein Condensation is achieved through RF evaporation of

magnetically trapped atoms, followed by evaporative cooling in an optical dipole trap

Details about the experimental apparatus beyond the brief overview found here can be

found in the PhD thesis of Roger Brown [48].

3.1 Vacuum System

3.1.1 Atom Source and Oven

Our vacuum system for the experiment consists of an oven chamber and a UHV cham-

ber. The oven chamber contains a 5g rubidium source within a 1.33” conflat bellows, which is

connected to the ”bright wall”, which is a 4” nipple kept at a higher temperature, designed to

prevent migration of the Rb metal to that port of the oven. While operating, the oven is kept at

85 C, and the bright wall is kept at 105 C. The bright wall is connected to a standard spherical

octagon chamber, which is attached to a Varian ion pump, which pumps at 125 L/s to achieve

pressures in the low 10−10 Torr. The atomic beam is directed from the source to a pinch off tube

44



within the bright wall through a Uniblitz shutter, which is mounted to another pinch-off tube with

a Swage-lok fitting, allowing for the atomic beam to be mechanically shuttered with an external

TTL pulse. The final piece of the oven system is the cold cup, which is a copper piece thermally

contacted through a special designed copper feed-through to a peltier unit on the outside of the

vacuum chamber. The hot side of the peltier unit is water cooled and flushed with nitrogen, while

the cold cup is kept at a temperature close to -18 C, which collects any of the atomic beam that is

not collected through the pinch off tubes. A hole in the cold cup allows for a collimated atomic

beam to pass through the cup and towards the main chamber. A pneumatic gate valve is used

to separate the oven from the Zeeman slower, which directs atoms from the source to the UHV

chamber.

3.1.2 UHV Side

The UHV chamber is another Spherical Octagon with 8 ports and 2 custom bucket win-

dows, designed to allow for closer placement of magnetic coils and optics and is attached to a

Varian ion pump (55 L/s), and an ion gauge. The top of the UHV chamber contains a titanium

cartridge which is used for titanium sublimation pumping. In order to avoid coating the win-

dows on the vacuum side, hand adjustable flags are installed that are rotated prior to pumping the

titanium, which cover and protect the windows.

3.2 Magnetic Fields

The experimental apparatus uses a magnetic quadrupole trap, which is generated from a

pair of anti-Helmholtz coils along the vertical axis. At the center of the trap, the magnetic gradient

45



is 1.05 G/A·cm while additional shim coils along all three directions provide a bias field and

allows for the center of the quadrupole trap to be shifted in space. To allow for fast switching,

the power supply for the quad coils, held at a constant 15V, capable of 440 A, is controlled using

a water cooled MOSFET bank, capable of switching the current in ≈ 1 ms. The water cooling

is interlocked to the current supply through flow switches, ensuring that the current supply is not

on while the water is off.

3.3 Laser Systems

3.3.1 Cooling Light

Our experiment relies on the cooling and trapping of 87Rb, specifically on the D2 line. The

relevant level structure is shown in Fig. 3.1. We initially lock our master laser to the crossover

peak between the |F = 2〉 → |F ′ = 2〉 and |F = 2〉 → |F ′ = 3〉 using saturated absorption

spectroscopy. The cooling light addressing the |F = 2〉 → |F ′ = 3〉 transition is beatnote

locked to the master source. The master is also used to beatnote lock the repump laser to the

|F = 1〉 → |F ′ = 2〉 transition. The cooling and repump locking systems are both controlled

by voltage controlled oscillators, which can be dynamically adjusted during an experiment and

the intensity of the light being sent to the experiment are controlled through AOMs. Additional

mechanical shutters allow for total blocking of light to the experiment.

3.4 Microwave and RF

In the presence of a magnetic field, the hyperfine states will split into 2F +1 zeeman levels,

meaning that the 5S1/2 level splits into 8 different zeeman levels, 3 from the F = 1 and 5 from

46



F = 1

F = 2

F' = 1

F' = 0

F' = 2

F' = 3

5S1/2

5P3/2

M
a
s
te

r

R
e
p
u
m

p

C
o
o
li
n
g
/I
m

a
g
in

g

S
lo

w
e
r
 C

o
o
li
n
g

S
lo

w
e
r
 R

e
p
u
m

p

Figure 3.1: 87Rb levels used for cooling and trapping.

the F = 2. The initial magnetic trapping of the BEC results in a pure |F = 1,mF = −1〉

condensate. The splitting between the F = 1 and F = 2 without the zeeman sublevels is ≈ 6.8

GHz, a frequency in the microwave regime. The use of microwaves allows for control over the

internal state of our atoms.

In order to reach a frequency of 6.8 GHz, we start with a Rhode and Schwartz frequency

synthesizer, which generates 13 dBm of 3.438 GHz microwaves. The signal is frequency doubled

and then mixed with an RF signal generated from a DDS device and sent to an amplifier, which

then passes through a stub tuner and then to a homemade copper horn. A TTL switch before the

amplifier allows for fast switching and control during experimental sequences.

The frequency seen by the atoms is controlled by the output of the RF signal from the DDS.

An example microwave spectroscopy is shown for the |F = 1,mF = −1〉 → |F = 2,mF = −2〉

transition with a 155 µs pulse. These spectra are fit to the standard Rabi formula (3.1) for a two

level system, where the only free parameters are the center frequency and the coupling strength.

47



47.780 47.785 47.790 47.795 47.800 47.805
W Frequency (MHz)

0.0

0.2

0.4

0.6

0.8

1.0

Fr
ac

tio
na

l P
op

ul
at

io
n 

in
 |2

,-2
> 

 W spectrum

Figure 3.2: Typical µW spectrum for the |F = 1,mF = −1〉 → |F = 2,mF = −2〉 transition.
The frequency is the value sent from the DDS device to be mixed with our microwave source.

Pop(ω; Ω, ω0, Tpulse) =
Ω2

Ω2 + (ω − ω0)2 sin2
(
πTpulse

√
Ω2 + (ω − ω0)2

)
(3.1)

3.5 Imaging

The experiment has two different imaging options. The first is a Point Grey Flea3, which

is set up to image perpendicular to direction of gravity. The second camera, a Princeton Images

PIXIS, is set below the atoms, allowing for analysis of the x-y plane after a time of flight (typically

27 ms).

Our images are obtained through standard absorption imaging. We take 3 sequential im-

48



ages: probe with atoms (ATOMS), probe with no atoms (NO ATOMS) and then a background

(BG), which has no probe light. Using Beers law, which states the absorbance scales as I(z) =

I0 exp(−ODz), we obtain the optical depth as

OD = − ln

(
ATOMS-BG

NO ATOMS-BG

)
(3.2)

3.6 Computer Control

The experimental sequences are programmed through Setlist, which is an in-house LabView-

based program that is responsible for generating the desired ramps and triggers. Analog/digital

outputs are generated through National Instruments DAQ devices, which communicate to our

Setlist through the National Instruments DAQmx driver. Calculations of ramp values are done on

the backend, either through LabView or MatLab.

3.7 Optical Lattice

The optical lattice is formed in a bow tie configuration, resulting in a 4-beam lattice as

shown in Fig. 3.3, which is a single beam retroreflected back on itself. The light is produced by

a tunable MBR Ti-Saph laser, which is pumped with a 10 W Verdi (532 nm) and fiber coupled

from a separate optical table to the experimental apparatus. An adjustable waveplate allows for

adjustment of powers between the 2D optical lattice and the vertical optical lattice.

The bowtie configuration is used due to its phase stability, despite having more than 3

(2D+1) beams. Any phase noise will act on all four beams simultaneously, resulting in a total

displacement of the lattice, with no change to the lattice unit cell structure. Our setup uses two

49



Pockels cells, which are electro-acoustic modulators (EOMs) that change the polarization of light

based on the voltage that is applied to the crystal. For our lattice setup, light in the x-y plane will

result in a 2D lattice with spacing λ/2, whereas light out of the plane (z direction) results in a

2D lattice with spacing λ. The total lattice seen by the atoms will be the combination of these

two lattices, as discussed in section 2.7.2. Computer control of the voltage on the input Pockels

cell sets the ratio of in-plane and out of plane light. The voltage applied to the 2D Pockels cell

translates the position of the out of plane lattice, allowing for different lattice geometries and

double well structures.

In order to precisely control the intensity of the lattice beams, a small amount of light is

picked off after the fiber and sent to a photodiode. The photodiode signal is sent to the DDS de-

vice, which has digital feed back circuitry to stabilize the lattice AOM at the desired (dynamically

adjustable) amplitude.

50



Amplitude 
Lock 
Feedback

Lattice AOM

INPUT PC

2D PC

HIGH 
VOLTAGE

HIGH 
VOLTAGE

Figure 3.3: Experimental configuration of the bowtie optical lattice.

3.7.1 Optical Lattice Calibration

In order to calibrate the depth of our optical lattice, we perform a standard Raman-Nath

diffraction pulse [49]. The Raman-Nath approximation is valid for the case when the lattice

pulse time, t, is shorter than the approximate harmonic oscillation period for the lattice, which

is a valid assumption for our lattice depths of up to 50ER and pulse times of 2µs. After pulsing

the lattice of depth V0 for a time t, the relative population of atoms that are transferred to the 1st

51



order diffraction peak can be approximated as

P1/P0 =
J0(V0t/2)

J1(V0t/2)
(3.3)

Figure 3.4: Example image of lattice diffraction as well as a cross-section indicating the relative
populations in the first order. In this image, the imbalance between the cross-sections along the
diagonals indicates that there are slightly different depths along the x and y directions. This can
be due to slight alignment imperfections or power imbalances.

Figure 3.4 shows a sample image of lattice diffraction as well as the line profile of the

atomic peaks after diffracting the lattice on for 2µs. The resulting integrated peak magnitudes

give the relative populations (the R fit value gives the relative population in the first order). By

solving this value for Eq. 3.3, we find that this particular run gave a lattice depth of approximately

72ER. By combining this with the signal on the photodiode, this gives a reference for lattice

intensity as a function of voltage, which we can control for experimental runs.

3.7.2 Tilt Calibration

In order to calibrate the offset between the checkerboard sublattices, we follow the tech-

nique described in [34]. For a checkerboard lattice with an offset ∆ between neighboring sites,

the potential is ramped on in a time that is adiabatic with respect to vibrational spacing, but not

52



adiabatic with respect to tunneling in the lattice. For our experiment, we ramp on the lattice to a

depth of 30ER in 200µs, which loads into both sublattices, even though the ground state for large

tilts has atoms in one sublattice. The wavefunction is therefore a combination of atoms occupy-

ing both sublattices and the different sublattices acquire a phase ∆t. By holding the atoms in the

lattice for a variable time, the phase of the wavefunction in the sublattices oscillates, giving rise to

oscillations in the amplitude of the diffraction peaks. We fit the visibility between the diffraction

peaks to a damped sine wave, from which we can determine the tilt between neighboring sites.

The tilt calibration procedure is typically done with large lattice depths of 30ER, which re-

sults in a frequency fast enough to measure before thermalization occurs and interaction induced

phase coherence is lost. In order to extrapolate the measured tilts, we assume a linear scaling

of tilt with the lattice intensity, which is an approximation, as single particle tunneling and mean

field effects in the lattice can result in slightly different tilts (we find that these higher order effects

are small, and the tilts agree with values obtained through parametric spectroscopy techniques,

to within experimental accuracy). We also note that this is for a BEC prepared at q = 0 meaning

that in the band picture, this gives the q = 0 gap and does not contain any information about the

bandwidths or gaps at non-zero q.

3.8 Tripod Piezo Mounts

The primary physics in this thesis relies on periodically shaking, or modulating the position

of the optical lattice. Given the bowtie lattice geometry is a single retro-reflected beam that is

phase stable, using a scheme such as modulating an AOM to adjust the frequency of one arm

would not work. The shaking performed relies on physically modulating the position of the

53



lattice by moving the mirrors as shown in Fig. 3.5.

Figure 3.5: Schematic showing the position of the two tripod piezos. The linear combination of
the two displacements results in an arbitrary displacement of the lattice. The arrows indicate the
lattice motion resulting from shaking the associated piezo mirror.

54



3.8.1 Lattice Displacement

Flat Mirror

Figure 3.6: Geometry of flat mirror displacement. When the piezo moves a distance d, the initial
optical path, as shown in red, gets an additional phase from the optical path l1 − l2.

When moving the flat mirror (M2), the path length can be seen in Fig. 3.6. After moving a

distance d, the total path length is given by ∆` = `1 − `2, where

`1 =
2d

cos(π/8)
(3.4a)

`2 = 2d sin(π/8) tan(π/8) (3.4b)

The total path length after a displacement, d, of the mirror is given by

∆` = 2d cos(π/8). (3.5)

55



In the four beam lattice, k1 is unaffected, k2 and k3 acquire a phase φ = k∆` and k4

acquires 2φ. The total intensity (recalling k1/2 = −k3/4 and |ki| = k) is given by:

I(x, y) = |x̂
(
eik2x+iφ + e−ik2x+iφ

)
+ ŷ
(
eik1y + e−ik1y+2iφ

)
|2 = 4+2 cos(2kx)+2 cos(2ky+2φ)

(3.6)

In terms of the mirror displacement, the intensity is

I(x, y)− I0 = 2 cos(2kx) + 2 cos(2k(y + 2d cos(π/8))) (3.7)

Periodically shaking the flat mirror translates the lattice along the y direction (direction of

k1/k4 lattice) and the physical displacement is 2d cos(π/8).

Retro Mirror

The geometry for the retro mirror (M1) is a bit simpler. Since the beam is normally incident,

when the mirror is moved by a distance d, the path difference is 2d and the phase is φ = 2kd. In

this case, k1 and k2 are unaffected, while k3 and k4 each acquire an additional phase φ. The total

intensity is

I(x, y) = |x̂
(
eik2x + e−ik2x+iφ

)
+ ŷ

(
eik1y + e−ik1y+iφ

)
|2. (3.8)

After some algebra, this becomes

I(x, y)− I0 = 2 cos(2kx− φ) + 2 cos(2ky − φ)

→ (cos(2k(x− d)) + cos(2k(y − d)))

(3.9)

This shows that displacing the mirror will result in the lattice moving along the (x+y)

direction.

56



1.0 0.5 0.0 0.5 1.0
X (a.u)

1.00

0.75

0.50

0.25

0.00

0.25

0.50

0.75

1.00

Y 
(a

.u
)

Lattice Shaking Trajectories
d1 = 1, d2 = 0
d1 = 0, d2 = 1
d2 = 2d1, = 5

4
d1 = d2, =

Figure 3.7: Shaking trajectories for different mirror driving parameter combinations.

Combined Shaking

Given the geometry of the lattice, we are unable to decouple the mirror shaking into two

orthogonal dimensions. In this section we look at the form of the shaking needed to produce

a circular trajectory. When the flat mirror is moved a distance d1, the lattice is translated a

distance δ1 = 2d1 cos(π/2) along the y direction. The retro mirror, when translated a distance

d2, results in a lattice displacement of δ2 = d2 along x and y. If both are periodically varied,

di = di cos(ωt + φi), the total displacement of the lattice can be expressed as (setting φ1 = 0

without loss of generality)

x(t)

y(t)

 =

δ1 cos(ωt) + δ2 cos(ωt+ φ)

δ1 cos(ωt)

 (3.10)

The position of the lattice in 2D can be arbitrarily set, and Fig. 3.7 shows a set of different

57



parameters that can be used to create different 2D trajectories. Circular shaking, which is partic-

ularly of interest, requires a relative phase of 1.25π between the two shaking of the two mirrors

with a ratio of the displacement amplitudes of δ2/δ1 =
√

2.

3.8.2 Design

To perform the shaking, we designed a translatable mirror, which consists of three individ-

ually controlled piezo actuators in a tripod configuration. This design was the focus of the 2018

paper, [50], and was inspired based on interaction with the Esslinger group. The use of three

piezos in the tripod configuration instead of a single one allows for control of steering errors,

which can become important when dealing either with sensitive systems, such as cold atoms, or

in setups where there are long optical paths, and slight angular deviations can result in large phys-

ical displacements of the beam. It was found by experience that the use of a single piezo mounted

to a mirror would result in a tilt as well as a displacement, resulting in undesired steering.

The system design can be seen in Fig. 3.8. Our design consists of a steel “bullet” which has

a mass of 310g, designed to be a mechanical insulator and remove low frequency resonances. The

choice of steel in lieu of potentially heavier materials such as tungsten, was made as a compromise

between density, cost and ease of machining. The backside of the steel support is cut such that

it can be placed into a standard 1-inch optical mount, a requirement for fine-tuning the beam

direction. We include a thin 1 mm, ceramic disk onto the face of the steel bullet, to act as an

electrical insulator between the piezo and the steel. The piezos are Noliac NAC2012, which are

3 mm x 3 mm in size, have a capacitance of 65nF, and a maximum displacement of 3.3µm. The

three piezos are epoxied to the ceramic insulator in a tripod configuration with a standard Loctite

58



Figure 3.8: Design of the tripod mirror system. (a) shows the heavy steel ”bullet” design that fits
into a standard mirror mount. There is a small ceramic layer in between the steel and piezo. (b)
shows the technical specs of the steel bullet.

epoxy from any hardware shop and the other side of the piezos were then epoxied to the mirror.

In an effort to make sure that the epoxy was fully cured, the epoxying was done over different

days. Each epoxy bond was given over 24 hours to cure, in order to make sure there was no free

movement, which could introduce unwanted noise into the device.

The independent control of the piezo displacements needed to control the beam steering

was provided by an electronic circuit that allows for the relative tuning of the three applied volt-

ages based on a trimpot. By adjusting the trimpots, the relative fractional displacement can be

adjusted, which allows for fine control of the reflected beam direction as the mirror is displaced.

59



The circuit is shown in Fig. 3.9.

High voltage driver

20 V/V
Input

LM

10kΩ

20kΩ

10kΩ

LM7171

2kΩ

2kΩ

LM7171

2kΩ

2kΩ

LM7171

2kΩ

2kΩ

LM7171

2kΩ

2kΩ

1Ω

1Ω

1Ω

1Ω

LF356

10kΩ

20kΩ

10kΩ

LM7171

2kΩ

2kΩ

LM7171

2kΩ

2kΩ

LM7171

2kΩ

2kΩ

LM7171

2kΩ

2kΩ

1Ω

1Ω

1Ω

1Ω

LF356

PZ1 PZ2 PZ3

A B

Figure 3.9: Electronics board for controlling the three tripod piezos. After the High Voltage
Amplifier, which gives a 20 V/V amplification, the fractional signal to each of the individual
piezos is controlled by the value of individual trimpots on the board.

3.8.3 Piezo Characterization

Frequency Response

Given that the typical energy scales of cold atom systems are in the kHz range, it is desirable

60



that the piezo design works in this range, with up to 20kHz response bandwidth being our goal.

To test the frequency response of the mirror displacement, which reveals mechanical resonances

that may be present in the system, we used a Bode 100 Network Analyzer to measure the response

to a given input signal. The setup used for measurement is shown in Fig. 3.10. The tripod mirror

was placed in a Michelson interferometer, which has arms of approximately 10 cm, short enough

to avoid any steering. The other arm of the interferometer had a separate piezo mirror, which was

used solely for feedback stabilization of the low frequency fringe drifts of the interferometer.

Figure 3.10: Experimental setup and mechanical response for mechanical resonances. The inset
in the lower panel indicates how flat the response is up to 40 kHz.

The response is remarkably flat for frequencies up through 40 kHz, which we deemed

suitable for the cold atom experiments that are performed later in this thesis.

61



Displacement

Characterization of the displacement of the tripod configuration was done using the same

interferometer setup as shown in Fig. 3.10, just with the phase locking mirror set at a fixed

position. By driving the tripod arm with a triangular voltage signal, we observed interference

fringes on the APD and counting how many of these fringes occur in a period determines the

physical displacement of the mirror An example measurement is shown in Fig. 3.11.

Figure 3.11: By driving the tripod mirror with a periodic triangular voltage, interference fringes
pass on the detector. This example shows the signal produced by driving the piezo at 1kHz, with a
maximum voltage of 160V. This particular measurement shows 12 fringes pass in a ramp, which
corresponds to a displacement of 2.34 microns.

Based on our experimental determination, we found that the maximum displacement is

approximately 2 µm at frequencies up to 8 kHz. To maximize the lifetime of the mirror, we

added a software control in the system such that the product of the voltage and the frequency

do not exceed a certain threshold. The displacement measurements over various frequencies are

shown in Fig. 3.12.

62



Figure 3.12: Measurement of piezo displacement vs drive voltage for various frequencies. We
observe that all frequencies in this range can achieve a linear displacement up to about 2 microns.

Steering

We performed test experiments to check that the steering of the tripod mirror was in fact

an improvement on a single piezo mirror. To check the static beam steering response, we place a

beam profiler at a distance of 12.4 m away from the testing mirror, which allows us to measure

small displacements from the steering angle. A DC voltage was applied and at each voltage, a

set of 10 images with 1s exposure are taken and averaged. The beam profiler performs a fit of the

centers using a 2D gaussian profile, and the extracted position allows us to determine the steering

angle. The results, shown in Fig. 3.13, show that just adding the tripod piezo design, with equal

voltages, can improve the steering by a factor of 20, and by optimizing the relative ratios of the

voltages driving the piezos, steering can be reduced further to a factor of 60.

The final test for calibrating the system is to determine the performance in the actual cold

atom experiment when driven with an AC signal. As discussed in section 3.7.1, lattice diffraction

provides a measure of the lattice depth seen by the atoms. The diffraction of the lattice depends

63



Figure 3.13: Measurement of DC steering angle vs drive voltage. (a) shows the experimental
setup for measuring the center of the optical beam. (b) indicates the improvement from fine-
tuning the relative voltages through the feed forward system, (c) A comparison of the two designs
with respect to just using a single piezo mirror.

64



on the co-propogating beam’s overlap with the atomic cloud. In the case of a large steering error

from the piezo, we expect that the lattice intensity seen by the atoms would decrease, as the beam

steers away from the position of the cloud. Using the atoms as a test for the steering, we re-

checked the DC voltage and compared the resulting lattice intensity, as found from the first order

diffraction. At the DC voltage, the single piezo quickly lost diffraction signal, with the lattice

depth decreasing to half the original value at half of the maximum displacement.

The tripod design with the feed forward activated was impressively stable up to the max-

imum displacement, as shown in Fig. 3.14. To test the steering for AC driving, we periodically

drove the piezo mirrors at a given frequency and performed a 2µs lattice pulse at various drive

phases. As shown in the figure, at an example drive of 2kHz, the lattice depth remained flat,

within error, throughout a full driving cycle of the piezo, which is indicative that the steering is

not an issue for the cold atom system.

Figure 3.14: Determination of lattice steering on the lattice by measuring intensity seen by atoms.
Left panel shows the results of applying a static DC voltage to the retro reflected mirror. The
performance for a single piezo is notably bad, while the tripod design is significantly better,
showing the same lattice intensity across the whole range of displacements. The right figure
shows the results of lattice diffraction at various phases of the shaking at a frequency of 2kHz.
The resulting lattice intensities show a very slight modulation slightly out of phase with the piezo
drive. The deviation across the piezo drive cycle does not exceed the error, so we conclude the
frequency dependent steering effects were small.

65



Chapter 4: Parametric Heating in Periodically Driven Optical Lattice

Heating in periodically driven systems in an interesting problem. In principle, periodic

driving provides an infinite source of energy that can be transferred to the system. The most likely

result is that the system heats towards an infinite temperature [14]. Remarkably, there are situa-

tions in which the drive induced heating does not immediately lead to the system reaching infinite

temperature, such as bichromatic drives [51] or prethermal and disordered states [17, 52, 53, 54].

The detailed understanding of heating in interacting many-body systems is an outstanding the-

oretical question. In this section, I address our studies of one such heating mechanism in these

periodically driven systems, which is the emergence of parametric instabilities.

4.1 Theory

The basis for theoretically understanding the emergence of the parametric heating mecha-

nism was discussed in [55]. The high level approach is to treat the time dependent problem at

a mean field level and then calculate time dependent fluctuations on top of the time dependent

mean field solution, using a “Floquet-Bogoliubov” treatment. The starting point is a tight bind-

ing Bose-Hubbard Hamiltonian with a periodically driven term [56]. We present the theoretical

approach using the example of a 1D optical lattice, but the results are easily extended to the 2D

case we studied experimentally. The Hamiltonian is given by:

66



H = −J
∑
<n,m>

a†nam + a†man +
U

2

∑
m

nm(nm − 1) +K cos(ωt)
∑
m

mnm, (4.1)

where J is the tunneling energy, U is the on-site interaction energy and K is the periodic drive

strength.

In the mean-field limit of a condensate, we replace the quantum field on site n, an, with a

constant classical field, αn/
√
N , representing the condensate, with N total atoms. Substituting

this approximation into the Heisenberg equations of motion, one obtains the equations for the

classical amplitudes:

i∂tαn = −J(αn+1 + αn−1) + g|αn|2αn +K cos(ωt)nαn, (4.2)

where g = Uρ is an effective interaction term with ρ being the number of atoms per site, and the

amplitudes obey the normalization
∑

i |αi|2 = N . As discussed in Section 2.1, the drive results

in a renormalized hopping amplitude, and the initial wave function a0
n(t) will be the solution of

the effective GPE

− JJ0(K/ω)(a0
n+1 + a0

n−1) + U |a0
n|2a0

n = µan. (4.3)

In this effective equation of motion, the time dependent drive term has been eliminated, but the

tunneling is renormalized to an effective Jeff = JJ0(K/ω)

The ground state condensate wave function will be a bloch wave with either q0 = 0 or

q0 = π, where the momentum of the ground state depends on the value of the effective tunneling.

If J0(K/ω) > 0, then the ground state momentum will be q0 = 0, but when J0(K/ω) < 0, the

67



band structure is inverted and the ground state wave function is the bloch wave at q0 = π.

Given the initial condensate wave function, one is interested in the stability of the fluctua-

tions, δan on top of the condensate a0
n:

an(t) = a0
n(t)(1 + δan(t)). (4.4)

This mode is plugged into the equation of motion (4.2) and linearized by ignoring terms propor-

tional to |δan|2, which results in coupled equations for δan(t) and δa∗n(t):

i∂tδan(t) =
−J
a0
n(t)

(
a0
n+1δan+1 + a0

n−1δan−1

)
+K cos(ωt)nδan + 2U |a0

n|2δan

+U |a0
n|2δa∗n − i

∂ta
0
n

a0
n

δan

(4.5)

as well as the complex conjugate. Recalling the initial state was given as a0
n = eiq0n, one can

remove the periodic driving term by performing a time dependent gauge transformation,

an(t)→ an(t)e−inα sin(ωt), (4.6)

where α = K/ω. The transformed amplitudes with the additional fluctuation are substituted into

Eq. 4.2

i ˙δan = −J
(
eizδan+1 + e−izδan−1

)
+ 2Uρδan + Uρδa∗n, (4.7)

where I introduce ρ ≡ |a0
n|2 and z ≡ q0 − α sin(ωt). We now perform the standard Bogoliubov

68



transformation on the fluctuation term δan

δan =
∑
q

uq(t)e
iqn + v∗q (t)e

−iqn, (4.8)

which results in a Bogoloubov-de Gennes (BdG) equation:

i∂t

uq
vq

 =

E(q, t) + g g

−g −E(q, t)− g


uq
vq

 , (4.9)

where the effective, time averaged ‘energy’ is given byE(q, t) = 4J sin(q/2) sin (q/2 + q0 − α sin(ωt))

and the interaction parameter g ≡ Uρ. The effective time averaged E(q, t) can be expressed as

−2J cos(q + q0 + α sin(ωt)) + 2J cos(q0 + α sin(ωt)), which is equivalent to the difference

between the Bogoliubov energy on top of the BEC and the energy of a BEC at q0. The matrix

appearing in this expression contains the periodic time dependence of the system. To study the

dynamics, the initial system is time evolved one period to determine the propagator matrix, U .

Following Floquet theory, the quasi-energies of the system, εq are related to the eigenvalues, λq

of the propagator through the relation λq = e−iεqT .

Considering the case where the quasi energy contains both a real and imaginary part εq =

rq + isq, the time evolution of the system will be e−irqt esqt, which displays oscillatory behavior

in the real part, but can exponentially grow (decay) in the case where sq is positive (negative).

These Lyapunov exponents, where the quasi energy contains a positive imaginary component, are

indicative of a dynamical instability, where there is an exponential growth of an excitation/mode

at a given q.

Since each mode can have an imaginary part that contributes to the dynamical instability, a

69



good measure of the system response is the growth rate of the most unstable mode,

Γmum = max
q
sq. (4.10)

The reasoning for this selection is that in the long time limit, the exponential behavior will be

dominated by the largest value, as shown in [55]. In realistic systems, this model is an oversim-

plification due to the interactions between the growing modes, which mix in other modes as time

evolves.

Following the techniques presented in [55, 57], we can map the BdG equations, 4.9, into a

parametric oscillator. To do so, we express the time averaged energy as

Eavg(q) =
√

4|Jeff| sin2(q/2)(4|Jeff| sin2(q/2) + 2g), (4.11)

which is obtained by time averaging the effective GPE, Eq. 4.3, with Jeff = JJ0(K0). As noted

in [55], this time averaged energy corresponds to the zeroth order Floquet quasienergy. The full

transformation can be broken into two parts, the first diagonalizing the time averaged component:

uq
vq

 =

cosh(θq) sinh(θq)

sinh(θq) cosh(θq)


u′q
v′q

 (4.12)

where cosh(2θq) =
(
4|Jeff| sin2(q/2) + g

)
/Eavg(q) and sinh(2θq) = g/Eavg(q). The second

70



transformation, which transforms into a frame rotating at 2Eavg, is given by

ũ′q
ṽ′q

 =

e2iEavg(q)t 0

0 1


u′q
v′q

 . (4.13)

The equation of motion can then be expressed as

i∂t

ũ′q
ṽ′q

 =

Eavg(q) +Wq(t) +
g

Eavg

 0 hq(t)e
−2iEavg(q)t

−hq(t)e2iEavg(q)t 0



ũ′q
ṽ′q

 (4.14)

which are equivalent to the equations describing a parametric oscillator, with Wq(t) a diagonal

matrix that time averages to 0 and hq(t) is a real-valued function which can be expressed as a

fourier series using the Jacobi-Anger identity, and the relation, J`(x) + J−`(x) = 0:

hq(t) = 4J sin2(q/2)
∞∑

`=−∞

J2`(K/ω)e2i`ωt, (4.15)

which contains only even harmonics of the drive frequency.

Using the analogy of a parametric oscillator, we can map the oscillator’s eigenfrequency,

ω0 to the time averaged dispersion, Eavg and use the substitution ωpo = 2ω (po indicating para-

metric oscillator). From the classical results, the oscillator will display a parametric instability

for ω = 2ω0. Given that each q has a different average E, each mode can be treated as an in-

dependent parametric oscillator. Additionally, because the Fourier expansion of hq contains an

infinite number of harmonics, the resonance will occur once any even harmonic, 2`ω, becomes

equal to 2E(q). For sufficiently low interaction strengths, one can approximate the function hq

71



to the second harmonic (` = 1 Fourier component). The last term in 4.14 can be expressed as

8JJ2(K/ω) sin2(q/2)g

Eavg

 0 cos(2ωt)−2iEavg(q)t

− cos(2ωt)e2iEavg(q)t 0

 (4.16)

From the classical parametric oscillator [58], and the substitutions above, the instability rate can

be expressed as

sq =
4JJ2(K/ω) sin2(q/2)g

Eavg

√
1−

(
(ω − Eavg)Eavg

4JJ2(K/ω) sin2(q/2)g

)2

, (4.17)

In the long time limit, the dynamics will be governed by the most unstable mode which is the

highest value of sq.

At first order, the emergence of a parametric instability can be seen in the case where

Eavg ≈ ω. In this scenario, our condensate effectively absorbs a drive photon and excites a

resonant Bogoliubov mode, which grows exponentially with a rate given by:

Γq =
4JJ2(K/ω) sin2(q/2)g

ω
(4.18)

At second order, the instability can occur in a region around the resonance point and the

maximally unstable mode has a rate given by

Γ =
|J2(K/ω)| g
|J0(K/ω)| ω

(√
g2 + ω2 − g

)
. (4.19)

In this treatment of the 1D shaken optical lattice, there is a cutoff in the drive frequency.

When ω becomes greater than Eavg, the q = π point will be maximally unstable, but an excitation

72



would violate energy conservation. In a system with weak confinement in directions perpendic-

ular to the lattice, a 1D lattice would consist of an unbounded energy spectrum in the transverse

degree of freedom. When the drive frequency exceeds the maximum excitation energy, the sys-

tem can become unstable at the maximum dispersion, namely Eavg(q = π), and the difference

will be absorbed into the transverse direction. This was the subject of [59], which measured

parametric instabilities in a driven 1D system.

4.1.1 Driving in 2D

Our experimental system allows us to realize a 2D optical lattice with arbitrary 2D driving

trajectories. In addition, the system is a lattice of tubes, meaning there is a transverse degree of

freedom along the z direction. In this section, we follow the methodology from [55, 60] and go

through the key results for the 2D system.

Similar to the 1D case, we consider the time dependent GPE on a square lattice in the tight

binding limit:

i∂tam,n = −J(am,n+1 + am,n−1 + am−1,n + am+1,n) + U |an,m|2am,n

+am,nK(m sin(ωt) + εn sin(ωt+ φ))

(4.20)

where ε, φ represent the shake trajectory. When ε = 0, we consider shaking along the x direction

only, and when ε = 1, a phase of 0 (π/2) results in diagonal (circular) shaking. The treatment

of the mode functions, is the same as in Eq. 4.9, but each drive will have a different dispersion,

E(q, t) and a different Fourier decomposition of hq. For the following, we consider just the case

when the normalized drive amplitude, K0 ≡ K/ω < 2.4, which corresponds to the effective

73



tunneling, Jeff being positive.

X Shaking

For the case of shaking along one single direction, ε = 0, the BDG equation gives the

following values for E(q, t) and hq(t)

E(q, t) = 4J
(
sin(qx/2) sin(qx/2−K sin(ωt)/ω) + sin2(qy/2)

)
, (4.21)

hq(t