Quantum correlations of ultracold atoms in optical lattices
byShuming Li
B.A., Nanjing University, 2004
M.S., Chinese Academy of Science, 2007

A thesis submitted to the
Faculty of the Graduate School of the
University of Colorado in partial fulfillmentof the requirements for the degree ofDoctor of Philosophy
Department of Physics
2014This thesis entitled:
Quantum correlations of ultracold atoms in optical lattices
written by Shuming Li
has been approved for the Department of Physics

Prof. Ana Maria Rey

Prof. John Bohn

Date

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

iii
Li, Shuming (Ph.D., Physics)
Quantum correlations of ultracold atoms in optical lattices
Thesis directed by Prof. Ana Maria Rey

An optical lattice is the periodic potential that atoms experience via the ac-Stark shift when
they are illuminated by counter-propagating laser beams that form a standing-wave pattern. Optical lattices have been widely used as a versatile platform in cooling, trapping, controlling atoms,
and also for the study of a variety of problems in physics. The many-body states of ultracold
atoms in optical lattices can be characterized by the quantum-correlations encoded in time-of-flightimages. In this thesis, we mainly discuss the use of correlations function as a natural framework
for characterizing quantum states in optical lattices.
The outline of the thesis is as follows. Chapter 1 gives a brief introduction to optical latticepotentials and ways for describing particles moving in a periodic potential. Chapter 2 explains theimportance of correlations and discusses common methods to detect them in cold atom experiments.
Chapter 3 presents work done to study the many-body Schrödinger equation in a quasiperiodic potential and discusses its connection with the Kolmogorov-Arnold-Moser (KAM) problem
of classical mechanics. We posed a possible visualization of such a connection in experimentally
accessible many-body observables. These observables are useful probes for the three characteristicphases of the problem: the metallic, Anderson, and band-insulator phases. In addition, they exhibit
fingerprints of nonlinear phenomena such as bifurcations and devil’s staircases. Our numericaltreatment is complemented with a perturbative analysis that provides insight on the underlying
physics. The perturbative-theory approach is particularly useful in illuminating the distinction
between the Anderson-insulator and the band-insulator phases in terms of paired sets of dimerizedstates.
Chapter 4 discusses several theoretical procedures developed to understand a recent experiment on macroscopic quantum self-trapping (ST) performed in a 2D optical lattice. Mean fieldiv
and truncated-Wigner-approximation (TWA) calculations are performed trying to reproduce theexperimental observations. The discrepancy between the theory and the experiment lead to the
hypothesis of a new type of ST caused by strong correlations. We analyze toy models to supportit.

Dedication

To my husband and my son.

vi

Acknowledgements

I want to thank all the people who made this dissertation possible. First and foremost I
would like to thank my advisor Ana Maria Rey for giving me the opportunity to join her group.
For the past six years, she gave me not only very good advice and detailed guidance, but also theencouragement I needed to feel that the work I was doing was important and worthwhile. It has
been a great pleasure to work with and learn from such an extraordinary person.
I also want to thank Prof. Indubala I. Satija from George Mason University and Prof. Charles
W. Clark from the Joint Quantum Institute and NIST in Maryland. Without their great ideas andcritical thinking, the work described in Chapter 3 would have been impossible.
I am also very thankful to Prof. David Weiss, as well as Aaron Reinhard, Jean-Felix Riou and
Laura Zundel from his group. They did the self-trapping experiment that is described in Chapter 4
while we did the theoretical counterpart. After a lot of discussions in Skype meetings and hundredsof emails back and forth, we finally brought forward a new mechanism of self-trapping.
I also owe special gratitude to everyone including former members in the group for the beneficial discussions, Javier Von Stecher (Seagate), Michael Foss-Feig (JQI), Chester Rubbo (CedarCrestone), Kaden Hazzard (Rice university), Salvatore R. Manmana (University of Göttingen),Andrew Koller, Bihui Zhu, and all the new members in the group. Especially, I want to thankSalvatore R. Manmana for his help on DMRG calculations in this thesis. Also, thanks for theuseful discussions with Murray Holland, Dana Anderson, Jinx Cooper, Dominic Meiser and Ryan
Wilson.

vii
At last, I owe my deepest thanks to my family — my grandmother, my parents, my sister, myson, and my husband. My parents are helping me by taking care of my son when I am working on mythesis. Without their support, the complete of this thesis would not have been possible.

viii

Contents

Chapter
1 Optical lattices
1.1

1.2

1.3

1

Light-matter interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.1.1

AC-Stark shift . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.1.2

Dissipative interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5

1.1.3

Gaussian beams . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.1.4

1D, 2D, and 3D optical lattices . . . . . . . . . . . . . . . . . . . . . . . . . .

6Single particle properties . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.2.1

Band structure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.2.2

Wannier functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.2.3

External potentials . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

1.2.4

Semiclassical dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

Bose-Hubbard model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
1.3.1

The superfluid to Mott-insulator transition . . . . . . . . . . . . . . . . . . . 15

2 Detection of correlations

18

2.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.2

Coherence, interference and correlations . . . . . . . . . . . . . . . . . . . . . . . . . 19
2.2.1

Coherence . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19

2.2.2

Interference . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

ix
2.2.3
2.3

Correlations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

Time-of-flight technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
2.3.1

Experimental evidence of superfluid to Mott-insulator phase transition . . . . 24

2.4

Noise correlations and the Hanbury-Brown-Twiss spectroscopy . . . . . . . . . . . . 25

2.5

Cold atoms in one dimension . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.5.1

Interaction strength γ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 30

2.5.2

Two-particle correlations g (2) . . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3 Quantum correlations in bichromatic lattices

33

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

3.2

Ultracold atoms in a 1D bichromatic optical lattice . . . . . . . . . . . . . . . . . . . 35

3.3

Our model system: The Harper equation, many-body observables and self-duality

. 36

3.3.1

The Harper equation. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.3.2

Many-body observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37

3.3.3

Self-duality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

3.4

Localization transition as a KAM-Cantori transition . . . . . . . . . . . . . . . . . . 40

3.5

Fingerprints of nonlinear phenomena in many-body observables . . . . . . . . . . . . 44

3.6

3.5.1

Fragmented Fermi sea . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

3.5.2

Bifurcations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

3.5.3

Devil’s staircases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55

4 Self-trapping dynamics in a 2D optical lattice

58

4.1

Introduction to self-trapping . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59

4.2

MQST: an overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

4.3

MQST in a 2D lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
4.3.1

Mean field: Coupled GPEs . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

4.3.2

General behavior: Three different regimes . . . . . . . . . . . . . . . . . . . . 68

x
4.3.3
4.4

4.5

4.6

Comparison with experiment . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

Beyond mean-field model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
4.4.1

Comparison between theory and experiment: phase fluctuations . . . . . . . . 77

4.4.2

Number fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84

4.4.3

Phase and number fluctuations . . . . . . . . . . . . . . . . . . . . . . . . . . 84

Beyond the aTWA—two coupled tubes . . . . . . . . . . . . . . . . . . . . . . . . . . 86
4.5.1

The DMRG method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

4.5.2

Short-time dynamics

4.5.3

Long-time dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

4.5.4

Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

Conclusion. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 95

Bibliography

97

Appendix
A Time-of-flight image in a single-particle expansion

103

B Mapping from position-to-momentum space104

C Perturbation theory and dimerized states

105

D Variational method for 1D optical lattice

109

xi

TablesTable
4.1

Set of parameters and different initial conditions used to study the dynamics. Ntotal
is the total particle number in the ladder. U is the on-site interaction strength, Jk
is the nearest-neighbor tunneling matrix element within each tube, and J⊥ is thetunneling matrix element between tubes.

4.2. . . . . . . . . . . . . . . . . . . . . . . . 88

Fitted values of η for different set of parameters used to investigate the dynamics: at
short times the t-DMRG dynamics are well captured by the aTWA when η is used
as a fitting parameter.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

xii

FiguresFigure1.1x)|
AC Stark shift induced by the atom-laser interaction. The energy shifts are ~ |Ω(~
4δ

2

2

x)|
for the |gi state (and −~ |Ω(~for the |ei state). When δ > 0 (the laser is blue
4δ

detuned), the minima of the lattice potential correspond to the places with minimum
laser intensity. In contrast, when δ < 0 (the laser is red detuned), the minima of thelattice potential correspond to the place with maximum laser intensity. . . . . . . . .

4

1.2

Optical lattice in different dimensions. Courtesy of Michael Foss-Feig. . . . . . . . .

7

1.3Band structure of an optical lattice with different lattice depths. . . . . . . . . . . .

9

1.4Band structure: orange curves—exact solution from Mathieu equation, black-dashed1.5

curves—from the tight-binding approximation. . . . . . . . . . . . . . . . . . . . . . 10
R
Comparisons of dx|w(x)|4 , for the red solid line, w(x) is approximated by theGaussian function ψ0 (x), and for the blue dashed line w(x) is calculated numericallyby its definition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.1

Schematic of the Hanbury-Brown and Twiss experiment, showing two independent
detectors 1 and 2 that record the intensities of radiations from two sources a and b
simultaneously. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

2.2

Schematic 3D interference-pattern (in the superfluid regime) with measured absorptionimages taken along x and y directions after some time-of-flight. From [7]. . . . . . . 25

xiii
2.3

Time-of-flight images. The figure shows the interference pattern for different latticedepths, (a) 0ER , (b) 3ER , (c) 7ER , (d) 10ER , (e) 13ER , (f) 14ER , (g) 16ER , and (h)20ER . As the lattice depth is increased, a quantum phase-transition from superfluid
to Mott-insulator occurs (The superfluid-Mott transition happens around 12-13 ER ).
The ballistic-expansion time was 15ms. The narrow peaks developed during the expansion demonstrate its phase coherence across the lattice at the moment of release.
Peaks appear at Reciprocal lattice vectors, with an envelope function |w̃(~q)|2 . From[7]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26

2.4

Illustration of the HBT interferometry and the origin of quantum correlations. (a),After a time-of-flight period, the absorption image of the atomic cloud is detected by
a CCD camera. Two pixels of the camera P1 and P2 are highlighted, each of which
records the atoms in a column along its line of sight. (b), When two atoms initiallytrapped at lattice sites i and j are released and detected at P1 and P2 simultaneously, the two indistinguishable paths, illustrated by solid and dashed lines, interfere
constructively for bosons or destructively for fermions. (c), The resulting probabilityof simultaneously finding an atom at each detector is modulated sinusoidally as afunction of d in the unit of l (black curve). The multiple wave generalization to aregular array of six sources with the same spacing is shown in green. a.u., arbitrary
units. From [34]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

2.5

Noise correlations of a bosonic Mott insulator released from a 3D optical lattice. (a)Single shot image of the two-dimensional density-distribution of a Mott-insulating
atomic cloud released from a 3D optical lattice. (b) Horizontal cross section throughthe center of the image in (a), and Gaussian fit to the average over many independentimages. (c) Spatial noise correlations obtained by analyzing the same set of images.
(d) Horizontal profile through the center of the pattern in (c). Those peaks are
separated by integer multiples of l = Mhta , reflecting the underlying lattice-structure.
From [34] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

xiv
2.6

Noise correlation of a fermionic band insulator. Each plot has a similar meaning as
in Fig. 2.5. Instead of the correlation “bunching” peaks observed in Fig. 2.5, thecorrelation function of fermioninc quantum-gas shows an antibunching effect. From[35]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29

2.7

Local pair correlation function g (2) ≡ g (2) (0) as a function of interaction strength γ.
The solid line shows the theory prediction from [38] . The plot is taken from [39]. . . 32

3.1

Return maps in real space for (a) and (c) λ = 10(red), 1(green), 0.5(blue) and the corresponding return map in momentum space (b) and (d) λ = 0.1(red), 1(green), 2(blue).
The filling factors are ν = 0.25 for the upper panels, and ν = 0.5 for the lower panels.
The insets are averaged over 50 random phases.

3.2. . . . . . . . . . . . . . . . . . . . 42

Return maps for the irrational filling factor ν = 1 − σ. Panel (a) is in real spacewith λ = 10(red), 1(green), 0.5(blue) and panel (b) is in momentum space with thedual values λ = 0.1(red), 1(green), 2(blue). The insets are averaged over 50 random
phases.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

3.3

Various boundaries of the fragmented Fermi sea. . . . . . . . . . . . . . . . . . . . . 45

3.4

Momentum distribution for (a)λ = 0.1, (b)λ = 1, (c)λ = 2. All figures are averaged
over 50 random phases to mimic a realistic experimental case. We can see that the
sharp features of these distributions survive phase averaging. . . . . . . . . . . . . . 47

3.5

Momentum distribution for fixed Q(k = 1045) , FM = 4181 as a function of filling
factors vs. the disorder parameter. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

3.6

Momentum distribution vs. the disorder parameter. (a) ν = 0.5, and (b) ν =FM −1 /FM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

xv
3.7

Momentum distributions as a function of the disorder parameter. (a) filling factor
ν = 0.5. Here we see the formation of two tongues that merge at the critical value
λ = 1. (b) ν = σ. When ν is irrational, the system is a band insulator and theArnold tongue-like structure disappears. We have checked for the existence of thesimilar structure in position space when converting λ to 1/λ. . . . . . . . . . . . . . 48

3.8

Bifurcation of local density at Fibonacci sites. The central curve (red triangle)
corresponds to quarter filling while the lowest (black blocks) and the topmost (bluedisks) curves, respectively, correspond to ν = 1/5, and ν = 1/3. Quarter filling isa special case where the local density at Fibonacci sites is 0 or 1 as λ → ∞. For
ν > 0.25, the Fibonacci sites are filled, while for ν < 0.25 Fibonacci sites are empty
as λ → ∞. We have checked for the existence of the same bifurcation phenomena in
the quasi-momentum distribution, but with the weak and the strong coupling limits
reversed ( λ → 1/λ).

3.9. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

Local density distribution as the disorder parameter and filling factor are varied for(a): Fibonacci site j = 377 and (b): normal site j = 100. We have checked theself-duality behavior in the quasi-momentum distribution with the weak and thestrong-coupling limits reversed (λ → 1/λ).

. . . . . . . . . . . . . . . . . . . . . . . 51

3.10 Noise correlation for (a)λ = 0.1, (b)λ = 1, (c)λ = 2. All plots display averaged
quantities over 50 random phases. We can see that the structure in the noise correlations survives phase averaging, as they do in the momentum distributions of Fig.
3.4. ∆(Q = 0) is not displayed in those plots. . . . . . . . . . . . . . . . . . . . . . . 53
3.11 This figure shows the evolution of ∆(Qn ) ( n = 0, 1) as a function of the filling
factor for different λ values. In the limit λ  1, ∆(Qn ) exhibits steps occurring at(m)

the filling factors νju (m = 1, 2, . . . ) (see text). In the limit λ  1, it acquires asinusoidal profile. The inset shows the self-similar nature of the steps at λ = 1. . . . 54

xvi
3.12 This figure shows a measure of the gaps in ∆(0) vs. λ for different system sizes, FM .
We deem that a gap occurs if |∆(0)(Np + 1) − ∆(0)(Np )| >  (we choose  = 10−15here). We can see that with increasing M , the measure of the gaps approaches a
step function with the step position at λ = 1.
4.1

. . . . . . . . . . . . . . . . . . . . . 55

Schematic of the experimental setup. A BEC of 87 Rb atoms was initially preparedin the crossed-dipole-harmonic trap. A 2D optical lattice along the x-y direction
was then adiabatically ramped on to create an array of quasi-1D tubes. The crossed
dipole-trap was then turned off, and the expansion dynamics in the presence of the2D periodic potential were investigated by direct imaging. The line of sight was at 45
degrees from the lattice directions (45 degrees from x and y). The 2D atom density
distribution was recorded after various expansion times. . . . . . . . . . . . . . . . . 64

4.2

The acceleration of xrms (ẍrms ) obtained from the mean-field calculations for all
lattice depths used in the experiment: (a) 7.25 ER , (b) 9.25 ER , (c) 11 ER and (d)13 ER . From the sign of ẍrms , the evolution can be separated into three different
regimes indicated by a solid green line and a dashed blue line. The three regimes
are discussed in detail in the text. MQST is signaled in the xrms as a negative
acceleration. It starts at the solid green line and stops at ∼ tc (indicated by thedashed blue line) when the acceleration curves cross 0 from below. The values ofU ρ(tc )/Ja3 for the four depths under consideration at tc are {0.4, 3.0, 3.0, 3.0}. Those−1
correspond to ρ(tc )/J = {60, 350, 320, 300} µm−3 ER, respectively.

4.3. . . . . . . . . 67

Time evolutions of: the kinetic energy along the lattice direction (Ekx , dashed orangeline), the kinetic energy along z (Ekz , dotted black line), and the interaction energy(Eint , solid red line). Here energy is in units of ER . This plot is computed using themean-field model for a 7.25 ER lattice. The vertical lines are at the same positions
as in Fig. 4.2 (a). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

xvii
4.4

Expansion dynamics predicted by mean-field calculations for the lowest lattice-depth(7.25 ER ): Panel (a) displays the evolution of the density profiles at z = 0, after
integrating along one transverse direction. Panel (b) displays the same density profiles, but after integrating along a direction that is 45 degrees between the x andy axes (as done in the experiment). Panel (c) shows transverse density profiles at
z = 0, during the MQST regime at t = 8.7 ms. The insets in Panel (a) and (b) are
intersecting profiles at t = 8.7 ms marked by white lines in the 3D density plots andhave the same units as in the 3D plots. . . . . . . . . . . . . . . . . . . . . . . . . . . 69

4.5

Panel (a) is a schematic picture of the tunneling in 2D optical lattices. In 2D, nosite is self-trapped in all directions, so atoms are not fully frozen. The arrows show
the directions that atoms may tunnel, while the crosses mean that tunneling along
those directions is forbidden by MQST. Panel (b) is a density (ρ at z = 0) contour
plot of a square fort-like barrier developed during the MQST regime.

4.6. . . . . . . . 72

The evolution of (a) σn2 and (b) Rz for 7.25ER from mean field calculations andthe experiment. Each black point is the average of 10 experimental measurements.
In Panel (a), the error bars represent random uncertainty, but an overall systematic
uncertainty of 0.5µm associated with the imaging resolution is not included. To
illustrate the dependence of the mean-field dynamics on the initial width, we show
the evolution of σn2 for different initial conditions for the 7.25ER lattice. Each line
corresponds to a particular initial condition, subject to the constraint of matchingaxial expansion rates at long times. In Panel (b), error bars are not shown, since weare mainly interested in the long-time dynamics when Rz is very large and the error
bars are negligible. Even with different initial conditions, the expansion rates alongz are almost the same, therefore in Panel (b) all lines are almost on top of each other. 74

xviii
4.7

Time evolution of (a) σn2 and (b) Rz for 13 ER from mean-field calculations and theexperiment. As in Fig. 4.6, the solid and dashed lines are simulation results, andeach corresponds to a particular initial condition. The long-time dynamics of σn2
are sensitive to the initial width of the cloud also at 13 ER . The black points areexperimental data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

4.8

Dimensionless coupling strength at t = 0 as a function of tube location along they axis, if the tubes are approximated to be independent. Here V0 = 13 ER . The
approximation is never fully warranted in the regime of this work, where the tubes
are manifestly not independent. Still, we suspect that this γ gives a sense of theinitial axial correlations in the relatively deep lattices. . . . . . . . . . . . . . . . . . 75

4.9Transverse density profiles at z = 0 during the MQST regime at t = 8.7 ms for the
13 ER lattice obtained from the aTWA results. . . . . . . . . . . . . . . . . . . . . . 78

4.10 Mean-field, aTWA, and experimental results for the evolution of σn2 for all lattice
depths used in the experiment. The lattice depths are (a) 7.25 ER , (b) 9.25 ER ,(c) 11 ER and (d) 13 ER . At low lattice depth, a small η seems to account for theobserved behavior at short times. The dashed blue line is the result of a mean-fieldcalculation of the dynamics with no random phase between the tubes (η = 0). The
dotted purple line is for η = 0.2, the dashed-dotted green lines are for (a) η = 0.4
and (b) η = 0.5, and the solid orange line is for η = 1. The dashed-vertical red
lines indicate the transition from MQST to BE predicted by the aTWA. The yellow
shadow regions indicate the tc inferred from the experimental data. . . . . . . . . . . 80
4.11 Evolution of Rz for all lattice depths used in the experiment. The dots show experimental data, and the solid and dashed lines show the predictions obtained from theaTWA for different η. The lattice depths are (a) 7.25 ER , (b) 9.25 ER , (c) 11 ER
and (d) 13 ER . Rz is obtained by fitting the density along z to a Thomas-Fermi
profile. The theoretical curves use the same parameters as the ones shown in Fig.
4.10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81

xix
4.12 ẍrms calculated using the aTWA for different lattice depths and for the corresponding
optimal η, (a) 7.25 ER , η = 0.4 and (b) 13 ER , η = 1. The boundaries between
MQST and BE are indicated by the dashed blue lines. In panel (a), the boundary
between EXP and MQST is indicated by the solid-vertical green line. In panel (b),because ẍrms is almost 0 at all times, it can be interpreted that the self-trapping
starts at t = 0 ms. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.13 Transverse density profiles: The solid and dashed lines are obtained from the aTWA
results for the V = 7.25 ER lattice and for η = 0.4. The profiles are taken at z = 0
and viewed transversely at an angle 45o from the lattice axis and normalized to1. The solid blue lines indicate the theoretical profiles before convolution, and thedashed black lines after convolution. The red points are averaged experimental dataof 10 measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82
4.14 Transverse density profiles: The solid and dashed lines are obtained from the aTWA
results for the V = 13 ER lattice and for η = 1. The profiles are taken at z = 0 andviewed transversely at an angle 45o from the lattice axis and normalized to 1. The
solid blue lines indicate the theoretical profiles before convolution and the dashedblack lines after convolution. The red points are averaged experimental data of 10
measurements. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
4.15 Time evolution of σn2 computed for a 7.25 ER lattice, and for η = 0. The dashed
red line shows the result when number fluctuations are neglected, and the solid blue√line shows the result when number fluctuations (≈ 1/ Nmn ) are included. The plot
clearly shows that number fluctuations suppress the expansion, but only by a small
amount. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

xx
4.16 aTWA results for the time evolution of σn2 for a 7.25 ER lattice, and for η = 0.35.
The dashed red line shows the result when only phase fluctuations are included, and√the solid blue line shows the result when both number fluctuations (≈ 1/ Nmn ) andphase fluctuations are included. The plot shows that number fluctuations slightly
suppress the initial expansion rate.

. . . . . . . . . . . . . . . . . . . . . . . . . . . 85

4.17 Schematic configuration of the two-tube system. . . . . . . . . . . . . . . . . . . . . 87
1
4.18 NN̈|t=0 (in units of Jk2 ) as a function of NN1 (t=0), J = 0.05, U Ntotal = 30. The solid
total
total

black line shows the GPE limit, and the dash-dotted gray line the hard-core boson
limit [see Eqn. (4.11)]. The red squares, magenta upper-triangles, blue diamonds,1
and orange disks show the initial curvature NN̈
|t=0 extracted from the t-DMRG
total

results at U = 1, U = 2, U = 3 and U = 5, respectively. The dashed black lines
show the aTWA solutions using η as a fitting parameter. The values of η are listed
in Table 4.2 . This plot is computed using the “case 2” parameters shown in Table
4.2.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

4.19 Evolutions of U N1 (in units of Jk ) and N1 as functions of time (in units of 1/Jk )calculated using the DMRG and GPE methods. In panel (a) and (b), the valuesof U are 1 and 5, respectively, while U Ntotal is fixed to 30. In panel (c) and (d),the values of U are 2 and 5, respectively, while Ntotal is fixed to 30. In each panel,
the solid line is the t-DMRG solution, and the dashed black line is the mean-fieldcalculation with no random phase (GPE). The dotted gray line is for (a) η = 0.3,(b) η = 0.7, (c) η = 0.4 and (d) η = 0.53.
(2)

4.20 gl. . . . . . . . . . . . . . . . . . . . . . . . 93

as a function of time in different regimes. The correlations are extracted from(2)

the t-DMRG results for the parameters describing case 2. Panel (a) shows g1 ,(2) L PL P
2
2i hn̂1,i i − n1,i , and Panel (b) shows g2 , N 2
i hn̂2,i i − n2,i . In each panel, fromN2
1

2

top to bottom, the lines correspond to U = 1, U = 2, U = 3, and U = 5, respectively. 93

xxi
B.1 Mapping from j to k , denoted by K(j) used to determine momentum-space observables from position-space observables and vice versa. Larger points are Fibonacci
sites. Here we choose FM = 87, and K(j) is the relation that maps j to k. . . . . . . 104
C.1 Panel (a) shows numerically that at the two red points, m1 = FM −1 and m2 = FM −2 ,Λm vanishes, so that degenerate perturbation theory is required. Panel (b) shows
the relations of Lm , Bm and m (see text) for φ = 3π
4 . The relations between Lm andm are very sensitive to the value of φ while the relations between Bm and m are notsensitive at all. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 107

Chapter 1

Optical lattices

An optical lattice is created by the interference of counter-propagating laser beams, whichgive rise to a spatially periodic intensity-pattern. The intensity pattern generates a potential forneutral atoms via the a.c. Stark shift of the atomic energy levels. Ultracold atoms in an opticallattice are analogous to valence electrons in a solid crystal, where the optical lattice potential mimics
the crystal lattice in a solid, and the atoms loaded in the lattice mimic the valence electrons. In
recent years, optical lattices have been widely used as a platform for studying various problems insolid state physics [1].
In contrast to solid-state materials where impurities and structural defects are inevitable,
optical lattices are rigid, free of defects, and fully controllable. The two important parametersthat characterize an optical lattice are the depth of the lattice potential wells, V0 , and the lattice
constant, a. The depth of the lattice can be tuned by changing the intensity of the laser, while
the lattice constant can be tuned by changing the wavelength of the laser or by changing therelative angle between the two laser beams. By interfering laser beams in different ways, atoms
can be confined in different types of geometries [2, 3]. By controlling the intensity of the trapping
lasers, the dimensionality of the system can be varied from 3D to 0D [4, 5, 6, 7]. All these features
make optical lattices a versatile arena for the investigation of nonlinear phenomena, nonequilibrium
dynamics, quantum phase transitions, quantum information, and so on [1, 8].
In this chapter, we introduce the basic theory of optical lattices and discuss the singleparticle physics of atoms in such a periodic potential. We follow with a brief introduction to the2Bose-Hubbard model.

1.1Light-matter interactions
Light-matter interactions are what generate the optical lattice potential. In the dipole ap-

proximation, the light-matter Hamiltonian is given by [9],~Hdip = −d~ · E,

(1.1)

~ is the electric field.
where d~ is the dipole moment of the atom, and E
This dipole interaction can be considerably simplified if the external field is weak and nearlyon resonance with an atomic transition. In this case, only the ground state and the near-resonanceexcited level have appreciable probability of occupation, and thus the complicated multilevel atomic
system can be treated as a two level atom with states |gi and |ei.
We denote the energy difference between the |ei and |gi states as ~ωa , the laser frequency as
ωo , and the detuning of the laser from the atomic transition as δ = ωo − ωa . If both the atom and~ can be replaced by,
the laser field are treated quantum mechanically, Er~ = ~eEi~ωo hΞ(~x)â + Ξ∗ (~x)â† ,
2

(1.2)

with ~e the unit polarization vector of the field, and Ξ(~x) the field mode evaluated at the atomic
position ~x (for example, for a 1D optical lattice potential, Ξ(x) = √2V cos(kx)). â is the annihilation
operator of photons.
The dipole moment can be written in matrix form as,ê · d~ = µge |gihe| + µeg |eihg|,(1.3)~ ≡ µ. The total Hamiltonian of the system can thus be written as,with µeg = ~e · he|d|gi

 r

i
~ωa  †
1
~ωo h
†
†
H=
âe âe − âg âg + ~ωo â â ++Ξ(~x)â + Ξ∗ (~x)â† (µ|eihg| + µ∗ |gihe|) , (1.4)
2220where âg and âe are, respectively, the atomic annihilation operators of states |gi and |ei.

3
When δ is much smaller than the laser frequency, |δ|  ωo , the rapidly oscillating terms,â|gihe| and â† |eihg| can be neglected, i.e., we use the rotating-wave approximation. Now only thenear-resonant processes need to be considered. In the rotating frame of the laser and by treating
the laser field classically, the Hamiltonian can then be reduced to ~Ω(~x)~δ  †~Ω∗ (~x) †âe âe − â†g âg +a†e ag +ag ae ,
22
2
qo hN iΞ(~x)µ (hN i = hâ† âi).
where Ω(~x) is the Rabi frequency given by ~Ω(~x) = −2 ~ω2
0
H≈−

(1.5)

|Ω(x)|2 is proportional to the laser intensity I(~x) = ~ωo chN i|Ξ(~x)|2 , and~2 |Ω(~x)|2 =1.1.1

2I(~x) 2
|µ| .
0 c

(1.6)

AC-Stark shift

If the detuning is large compared to the Rabi frequency, i.e., |δ|  Ω, the atom-laser interaction gives rise to energy shifts in states |gi and |ei. Using second-order perturbation theory, theenergy shifts are∆E = ±~|Ω(~x)|2
,
4δ

(1.7)

with the plus sign for the |gi state and the minus sign for the |ei state. The ac-Stark energy shift
2x)|of the atoms, ~ |Ω(~
4δ , gives rise to the optical lattice potential in the ground state. The sign of thepotential depends on the sign of the detuning. See Fig.1.1.
When δ > 0 (the laser is blue detuned), the dipole interaction repels atoms out of the field,
and the potential minima correspond to the places with minimum intensity. In contrast, when
δ < 0 (the laser is red detuned), the dipole interaction attracts atoms into the light field, and thepotential minima are therefore found at positions with maximum laser intensity. Therefore, atoms
accumulate at the bright parts in a red-detuned lattice and at the dark parts in a blue-detuned
lattice.

4

δ>0

δ<0
e

g
δ



e

(a)

 2
Ω( x)
4δ

−δ

−

 2
Ω( x)
4δ

g

(b)
2

x)|
Figure 1.1: AC Stark shift induced by the atom-laser interaction. The energy shifts are ~ |Ω(~
for
4δ
2

x)|
the |gi state (and −~ |Ω(~for the |ei state). When δ > 0 (the laser is blue detuned), the minima
4δ
of the lattice potential correspond to the places with minimum laser intensity. In contrast, when
δ < 0 (the laser is red detuned), the minima of the lattice potential correspond to the place withmaximum laser intensity.

5
1.1.2

Dissipative interactions

In real atoms, besides the atomic transitions induced by the laser, atoms in the excitedstates also spontaneously emit photons and decay to the ground state. This process destroys theinternal coherence between the |gi and |ei states of the atoms. This effect can be taken into
account phenomenologically by attributing to the excited state a complex energy. Denoting the
spontaneous-decay rate of the excited state into the ground state as Γ, when δ  Γ the energy ofthe perturbed ground state becomes a complex quantity that can be written as [10, 11],E =V (~x) =Γsc =~ |Ω(~x)|2= V (~x) + iΓsc (~x),
4 δ − iΓ~ |Ω(~x)|2 δ|Ω(~x)|2≈~,
4 δ 2 + Γ24δ
~ |Ω(~x)|2 Γ|Ω(~x)|2 Γ≈~.
4 δ 2 + Γ2
4δ 2

(1.8)
(1.9)
(1.10)The real part of the energy corresponds to the optical potential.
The damping rate Γ can be determined from the theory of spontaneous emission as [10]Γ=ω03
|he|µ|gi|2 .
3π0 ~c3(1.11)The general expressions for the dipole potential and the scattering rate can be rewritten interms of Γ and I(~x),
V (~x) =Γsc (~x) =3πc2 Γ
I(~x).
2ω03 δ
 
3πc2 Γ 2
I(~x).
2~ω03 δ(1.12)(1.13)The dipole potential scales as I/δ, whereas the scattering rate scales as I/δ 2 . Therefore optical
dipole traps usually use large detunings and high intensities to keep the scattering rate as low as
possible at a certain potential depth [10].

6
1.1.3

Gaussian beams

A typical laser beam has a Gaussian intensity-profile. For example, the intensity of the laser
propagating along the z direction has the form,I(r, z) =2P
exp[−2r2 /w2 (z)],πw2 (z)where P is the total power of the laser beam, r =(1.14)

px2 + y 2 is the radial distance from the beam’s

center, and w(z) is the 1/e2 radius. For a beam with wavelength λ, the variation of the spot sizeq
2 , where w is called the beam waist, and z = πw 2 /λat point z is given by w(z) = w0 1 + z 2 /zR0R0is called the Rayleigh length. Around the intensity maximum, the ac-Stark shift for a red-detuned
laser beam is approximated by the harmonic potential
Vdip (r, z) ≈ −V0 [1 − 2(r/w0 )2 − (z/zR )2 ].

(1.15)

Here V0 is the trap depth, and it is linearly proportional to the laser power. V0 typically ranges froma few kHz up to 1 MHz, i.e., from the nanokelvin to the microkelvin regime. The harmonic conp
finement is characterized by radial (ωr ) and axial (ωz ) trapping frequencies, with ωr = 4V0 /M w02q
2 , where M is the mass of the atom.
and ωz = 2V0 /M zR
1.1.4

1D, 2D, and 3D optical lattices

The simplest lattice-potential is a one-dimensional lattice formed by two counter-propagating
beams. For a Gaussian beam, the typical value of zR is significantly larger than w0 , so around theintensity maximum, the trapping potential has the simplified form of2

2

Vdip (r, z) ≈ −V0 e−2r /w0 sin2 (kz),(1.16)

where k = 2π
λ is the wave vector of the laser light, and V0 is the maximum depth of the latticepotential. The 1D optical potential generates an array of disk-shaped potentials (pancakes).
If two additional counter-propagating beams are added perpendicularly to the original 1D
lattice (the polarization vectors of the second-set of beams are also orthogonal to the first set of

7

Figure 1.2: Optical lattice in different dimensions. Courtesy of Michael Foss-Feig.

beams), the resulting trapping potential has the form,
2

2

2

2

2

2

Vdip (x, y, z) ≈ −Vx e−2(y +z )/wx sin2 (kx) − Vy e−2(x +z )/wy sin2 (ky)




z2
x2
z2
y2
2
≈ −Vx 1 − 2 2 − 2 2 sin (kx) − Vy 1 − 2 2 − 2 2 sin2 (ky),
wx
wx
wy
wy

(1.17)

where wx and wy are, respectively, the beam waists along x and y. The above 2D optical potential
generates an array of elongated 1D traps (tubes).
If a third pair of counter-propagating beams are added to the 2D optical lattice, the resulting
trapping-potential has the form in the case where any additional potential from the Gaussian shapeof the laser intensity is neglected
Vdip (x, y, z) ≈ −Vx sin2 (kx) − Vy sin2 (ky) − Vz sin2 (kz).
Schematic configurations of optical lattices in 1D, 2D, and 3D are shown in Fig. 1.2.

(1.18)

8

1.2

Single particle properties

1.2.1

Band structure

Considering a single atom in an 1D optical lattice potential, the single-particle Hamiltonian
is
H0 =

p2
+ V0 sin2 (kx).
2M(1.19)For convenience, the depth of the lattice potential V0 is generally specified by the recoil energy
ER = ~2 k 2 /2M , which is the kinetic energy gained by an atom with M by absorbing one photon.
The typical value of a lattice depth ranges from several recoil up to hundreds.
From condensed matter theory, it is known that in such a periodic potential, the eigenstates(n)

are the so-called Bloch functions φq (x) [12], whereH0 φq(n) (x) = Eq(n) φ(n)q (x).

(1.20)(n)Bloch functions have the form of a plane wave eiqx times a periodic function uq (x) with
(n)

(n)

uq (x) = uq (x + a). Here a is the lattice spacing and q is a quantum number that characterizes
the translational symmetry of the periodic potential. ~q is called quasimomentum, and it plays asimilar role in a periodic potential as momentum does in free space. Each Bloch function has a
(n)

(n)

specific value of q and eigenvalue Eq , where n is an integer, called the band index. If φq (x) is(n)an eigenfunction of H0 , so is φq+Gm (x), where Gm = m 2π
a is called the reciprocal lattice vector
with m an arbitrary integer. In general, the wave vector q is confined to the first Brillouin zone,(n)i.e., q ∈ (−π/a, π/a]. If we plot Eq

as a function of q, we will get a series of smooth curves

from q ∈ (−π/a, π/a]. Each curve corresponds to a specific band with band index n. Larger ncorresponds to higher band with larger energy.
(n)

(n)

In the case of a sinusoidal potential, the Bloch functions φq (x) and their eigenvalues Eq(n)can be solved analytically. The Bloch functions φq (x) can be expressed in terms of Mathieu(n)functions, and the eigenvalues Eq

be written in terms of Mathieu characteristic values. Fig. 1.3

shows the band structure of an optical lattice for different lattice depths. For V0 = 0, particles

8

8

6

6

4
2

4
2

0
-1.0 -0.5

0.0
qaΠ

0.5

0
-1.0 -0.5

1.0

V0=4 ER
8
6
4
2
0-1.0 -0.5

0.0
qaΠ

0.5

1.0

0.5

1.0

V0=20 ER
15
EqER

EqER

9

V0=1 ER
EqER

EqER

V0=0 ER

0.0
qaΠ

0.5

1.0

10
5
0
-1.0 -0.5

0.0
qaΠ

Figure 1.3: Band structure of an optical lattice with different lattice depths.

move in free space, and the spectrum is quadratic in q. As the potential depth is increased, theband structure appears. For small V0 , a gap opens at the edge of the first Brillouin zone. As thedepth increases, the band gap enlarges, the bandwidth shrinks, and more and more gaps open at
the edge of the first Brillouin zone. For very deep lattices, the spectrum is almost degenerate in q
and exhibits a dependence on n similar to the one of a particle in a harmonic trap.

1.2.2

Wannier functions

In many cases, it is more convenient to use another complete set of orthonormal wave functions, which are called Wannier functions, to describe atoms in a optical lattice. As Fourier transforms of Bloch functions, Wannier functions constitute an equivalent basis to describe the periodic
system, i.e.,wn (x − xi ) =(n)Replacing φq (x) with √1NPie

iqxi w

1 X −iqxi (n)√eφq (x).
N q(1.21)n (x − xi ) in equation (1.20) and using the orthonormal

4

4

3

3

2

2
1

1
0
-1.0 -0.5

0.0
qaΠ

0.5

0
-1.0 -0.5

1.0

V0=4 ER

0.0
qaΠ

0.5

1.0

0.5

1.0

V0=20 ER

4

8

3

6
EqER

EqER

10

V0=1 ER

EqER

EqER

V0=0 ER

2

4
2

1
0
-1.0 -0.5

0.0
qaΠ

0.5

1.0

0
-1.0 -0.5

0.0
qaΠ

Figure 1.4: Band structure: orange curves—exact solution from Mathieu equation, black-dashedcurves—from the tight-binding approximation.

properties of Wannier functions, we get
Eq(n)

=

X
Zi

≈

iqxi

e

Z

dx wn∗ (x)H0 wn (x − xi )dx wn∗ (x)H0 wn (x) + eiqa

= n − 2J (n) cos(qa) + · · · .

Zdx wn∗ (x)H0 wn (x − a) + e−iqa

Zdx wn∗ (x)H0 wn (x + a) + · · ·(1.22)The last line is called the tight-binding approximation.1 The tight-binding approximation is widelyused for atoms in deep lattices.
As we mentioned in the last section, the band structure of an atom in a sinusoidal potential
can be solved exactly. Therefore, we can get the value of J and n in Eqn. (1.22) from the
bandwidths and the band gaps of the exact band-structure. On the other hand, the validity ofthe tight-binding approximation can be checked by comparing the exact band-structure with the1
The terminology of “tight binding” was first used in condensed matter physics for calculating the electronic band
structure of tightly bound electrons in solids. When an electron is tightly bound to the atom to which it belongs andthe effect of surrounding atoms on the electrons is limited, the Wannier functions are replaced by atomic orbitals of
the free atom to which the electron belongs. In the tight-binding approximation, usually only the on-site energy En
and first order correction due to the nearest-neighboring tunneling are included, other terms are neglected.

11band-structure in the tight-binding approximation, as shown in Fig. 1.4. We can see that the tightbinding approximation is a good approximation when the lattice depth is deep enough. Provided
atoms only occupy the lowest band, the tight-biding approximation is good for a lattice depth
larger than 4ER .

1.2.3

External potentials

The tight-binding approximation is commonly used to solve the problem of a particle in aperiodic potential when an external field is also applied.
When the lattice is deep enough and the external perturbation is not strong enough to induceinter-band transitions, one may only expand the moving-particle wave function in the Wannier basesin the first band, i.e.,ψ(x, t) =Xzi (t)w0 (x − xi ),(1.23)iwhere i is the lattice site index, and xi is the position of the ith lattice-site. The Schrödinger
equation for a particle that moves in the periodic potential of a 1D lattice Vlat plus a perturbative
external-potential Vext (x) is∂~2 ∂ 2i~ ψ(x, t) = −+ Vlat (x) ψ(t) + Vext (x)ψ(x, t).
∂t
2M ∂x2(1.24)

By substituting Eqn. (1.23) in the Schrödinger equation, we get the following equations of
motion,i~X∂zi (t) = −Jij zj (t) + Vij (x)zj (t),∂t(1.25)

j

where
Z
Jij
Vij

dx w0∗ (x − xi )~2 ∂ 2−+ Vlat (x) w0 (x − xj ),2M ∂x2= −Z=
dx w0∗ (x − xi )Vext (x)w0 (x − xj ).

(1.26)(1.27)

If the external perturbation is a slowly varying function, it can be assumed constant insideeach individual lattice-site. And the Wannier function w0 (x − xi ) can be replaced with the delta

12function δ(x − xi ), and then Vij = Vext (xi )δij . Moreover, since the lattice is deep enough, only thetunneling to nearest-neighbors (for j = i ± 1) is relevant. The equations of motion reduce toi~∂zi (t) = −J [zi−1 (t) + zi+1 (t)] + 0 zi (t) + Vi zi (t),∂t(1.28)with
~2 ∂ 2+ Vlat (x) w0 (x + a),= −−2m ∂x2
Z~2 ∂ 2∗+ Vlat (x) w0 (x),=dx w0 (x) −2m ∂x2Z
dx w0∗ (x − xi )Vext (x)w0 (x − xi ) ≈ Vext (xi ).
=
Z

J
0
Vi

dx w0∗ (x)(1.29)(1.30)(1.31)

Here J is the tunneling matrix element between nearest-neighboring sites, 0 is the on-siteenergy-shift, and Vi is the external potential at lattice-site i. Equation (1.28) is called the discrete
Schrödinger equation (DSE).

1.2.4

Semiclassical dynamics

The wave-packet solutions of the Schrödinger equation behave like classical particles obeying
the classical Hamiltonian equations,
ẋ =

∂H
,
∂p

ṗ = −∂H.
∂x(1.32)

Bloch’s theorem introduces a wave-vector q that turns out to play the same role in thegeneral problem of motion in a periodic potential as the momentum p/~ does in free space. So
the semiclassical equations 2 of a particle with momentum ~q in the first band can be written as[12, 13],(0)ẋ = v (0) (q) =

1 dEq
,
~ dq

ṗ = ~q̇ = −dVext (x).
dx(1.33)

The first equation simply states that the velocity of a semiclassical atom is the group-velocityof the underlying wave-packet. The second equation states that the rate-of-change of an atom’s
2

The semiclassical model is a partial classical limit in the sense that the externally applied fields are treatedclassically, but the periodic potential is not [12].

13quasi-momentum is given only by the external force, but not by the periodic field of the lattice.
In sum, the semiclassical equations-of-motion describe how, in the absence of collisions and in thepresence of an external potential, the position x and the wave vector q of a particle evolve entirely
in terms of the band structure of the lattice.
If we compare the acceleration predicted by the model with the conventional Newtonian
equation, M ẍ = −dVext (x)/dx, we can associate an effective-mass m∗ induced by the presence ofthe lattice. This is given by,1
1 d2 (0)
~2
=
E
=
.
q
m∗
~2 dq 22Ja2 cos(qa)(1.34)

2~
Close to the bottom of the band, q → 0, m∗ → 2Ja
2 . Since the tunneling decreases ex-

ponentially with the lattice depth [11], the effective mass grows exponentially. In the presence ofinteractions, the large effective-mass manifests itself in a substantial enhancement of the interaction
to kinetic energy ratio, in comparison to the free particle case. Therefore atoms in optical lattices
can easily reach the strongly interacting regime, in which the interaction energy of the atoms at agiven density dominates over their characteristic quantum kinetic-energy.

1.3

Bose-Hubbard model
The simplest nontrivial model that gives an approximate description of the physics of inter-acting bosons in a lattice is the Bose-Hubbard model [14]. It is closely related to the Hubbard
model [15], which originated in solid-state physics as an approximate model to describe the motionof electrons in a crystal lattice. The Bose-Hubbard Hamiltonian (for bosons) exhibits a quantum
phase-transition from a superfluid state to a Mott-insulator state. Following the proposal by Jaksch
et al in 1998 [16], the Bose-Hubbard model has been realized experimentally using ultracold atomic
systems in various dimensions [7, 17, 18, 19].
The Bose-Hubbard model can be derived using second quantization which describes interacting atoms in an external trapping-potential. When the lattice is deep enough, and the external
perturbation is too weak to induce inter-band transitions, in the language of second quantization,

14
the field operator can be expanded in terms of Wannier orbitals in the first band and their associatedannihilation-operators, i.e.,Ψ̂ =Xâi w0 (~x − ~xi ),(1.35)iwhere i is the lattice site index, and ~xi is the position of the ith lattice site. The Bose-Hubbard
model can be derived by substituting Eqn. (1.35)


Z
~2 2
d~x Ψ̂ (~x) −∇ + Vlat (~x) Ψ̂(~x) + d~x Vext (~x)Ψ̂† (~x)Ψ̂(~x)
2m
Z
Z
1
0
†
† 0
0
0
d~xd~x Ψ̂ (~x)Ψ̂ (~x )Uint (~x − ~x )Ψ̂(~x)Ψ̂(~x ) − d~x µ(~x)Ψ̂† (~x)Ψ̂(~x).
2

Z
Ĥ =
+

†

(1.36)As mentioned in the last section, Vlat (~x) is the periodic lattice potential, Vext (~x) denotes any
additional slowly varying external-potential, and µ is the chemical potential. Uint (~x − ~x0 ) is thetwo-body interatomic potential. When the atomic distance is much larger than the effective extent
of the Uint , the atomic interaction due to s-wave scattering can be well approximated by a contact
interaction, i.e.,
Uint (~x − ~x0 ) ≈ gδ(~x − ~x0 ),g=

4πas ~2,M(1.37)

where as is the s-wave scattering length.
In the tight-binding approximation, we get the Hamiltonian in the second quantization form
H = −JX †X
UX † †âi âi âi âi ,âi âj +(Vi − µ)â†i âi +2(1.38)

i

i

hi,ji

with J and Vi as in Equation (1.29, 1.31) andU=

4πas ~2MZ

dx|w0 (x)|4 .

(1.39)

The first term in the Hamiltonian proportional to J is a measure of the kinetic energy ofthe system. Next-to-nearest neighbor tunneling amplitudes are typically two orders of magnitude
smaller than nearest-neighbor ones, and to a good approximation they can be neglected [11]. The
second term in the Hamiltonian takes into account the energy offset between lattice sites due toa slow-varying external potential Vext (x), and chemical potential µ. The third term accounts forinteratomic interactions. The parameter U measures the strength of the repulsion between two15aÙ wHxL¤4âx
ER
4
3
2
1
10

20

30

40

50

60

70

V0

R
Figure 1.5: Comparisons of dx|w(x)|4 , for the red solid line, w(x) is approximated by the Gaussianfunction ψ0 (x), and for the blue dashed line w(x) is calculated numerically by its definition.

atoms at the same lattice-site. While J decreases exponentially with lattice depth V0 , U increases
D/4

as a power law, V0, where D is the dimensionality of the lattice.

The integral of U is not sensitive to the oscillatory tails of the Wannier orbitals, thereforeRwhen calculating dx|w0 (x)|4 , the Wannier functions can be approximated by the eigenfunctionsq 2 1/4 − x2 2xof the harmonic trap, i.e., ψn (x) = √21n n! π0
e 2x0 Hn xx0 with x0 = M~ωx the characteristic length of an harmonic potential. In Fig. (1.5), we compare the numerical values betweenRR
dx |w0 (x)|4 and dx |ψ0 (x)|4 . We can see that these two curves are very close to each other,which means that Gaussian functions are good approximations of Wannier orbitals in estimating
U.

1.3.1

The superfluid to Mott-insulator transition

Here we skip the details about the phase diagram of the Bose-Hubbard model, but only
introduce the superfluid and Mott-insulator states, since we will use them a lot in this thesis.
At zero temperature, the physics described by the Bose-Hubbard Hamiltonian can be divided into two different regimes in the homogeneous case when Vi = 0. The first regime is the
kinetic-energy-dominated regime, when tunneling overwhelms the repulsion, and the atoms exhibit

16
superfluid properties. Another regime is the interaction-dominated regime when J is much smaller
than U , and the atoms behave like an insulator, the so-called Mott-insulator. The transition between the two regimes is a consequence of the competition between the kinetic energy, which tries
to delocalize the particles, and the interaction energy, which tries to localize them and penalizes
multiple-occupied lattice sites.
In the superfluid regime, the kinetic energy term dominates the Hamiltonian, and the systembehaves as a weakly interacting Bose gas. Quantum fluctuations can be neglected, and the system
can be described by a macroscopic wave function. For a translationally invariant lattice, the groundstate consists of all atoms in the zero-quasi-momentum mode. In this case,1  † N
|ΨSF i = √
b̂0
|0i,
N!

(1.40)

where N is the total number of atoms, |0i is the state of vacuum, and b̂†q = √1L† iqxjis thej âj ePcreation operator of an atom with quasi-momentum q, and L is the total number of lattice sites.
With increasing interatomic interactions, the average kinetic-energy required for an atom to
hop from one site to the next becomes insufficient to overcome the interaction-energy cost. Atoms
tend to get localized at individual lattice-sites, and fluctuations of the atom-number at a given
lattice-site are penalized. In this strongly interacting limit, if the number of atoms is commensuratewith the number of lattice sites, the ground state of the system enters the Mott-insulator phasecharacterized by localized atomic wave functions and a fixed number of atoms per site, i.e.,|ΨM I i =Y 1  † n√âj |0i,n!
j(1.41)

where n = N/L is the filling factor or mean number of particles per site. We have assumed thatthe number of particles is commensurate with lattice sites, i.e., that n is an integer.
For fermionic atoms with two possible spin-states | ↑i, | ↓i confined to the lowest band in 1D
optical-lattices, the system can be described by the single-band Hubbard model [1],H = −J

X †X
X
ĉi,σ ĉi+1,σ + ĉ†i+1,σ ĉi,σ + U
n̂i,↑ n̂i,↓ +Vi n̂i,σ .
i,σ

i

i,σ

(1.42)

17
Similar to the Bose-Hubbard model, the Hubbard model of fermonic atoms in an opticallattice describes a quantum phase-transition from a superfluid state to a fermionic Mott-insulator
state [20].
For spin-polarized fermionic atoms, because of the Pauli-exclusion principle, the s-wave interaction between two atoms is zero, therefore the low-energy physics is described by the Hamiltonianof free atoms moving in an optical lattice in the presence of an additional external potential,
H = −J X
X †ĉi ĉi+1 + ĉ†i+1 ĉi +Vi n̂i .
i(1.43)i

In the case of Vi = 0, the atomic eigenstates are the Bloch states. As the Puali principle does
not allow two fermionic particles to occupy the same Bloch state, atoms occupy the Bloch states
one by one, starting from the bottom of the lowest-energy band to the top. The full occupation of
the lowest band and the energy gap to the empty second band characterizes the fermionic bandinsulating state.

Chapter 2

Detection of correlations

2.1

Introduction
Correlation functions are a natural framework for characterizing quantum states. This chap-ter starts with a general discussion about the basic concepts of coherency, interference, and correlation. We then introduce the amplitude interferometry such as Young’s double slit experiment,
and Hanbury Brown and Twiss (HBT) intensity-interferometry. Next the investigation of quantum
correlations and their detection in ultracold gases in an optical lattice is discussed. For example,
absorption images after time of flight (TOF) provide a useful tool in determining correlation functions. The reason is that after the atoms are released from the trap and allowed to expand until thesize of the cloud becomes much larger than the original separation between the atoms, the momentum distribution of the trapped gas gets mapped into the density distribution in TOF images after
the expansion. The latter displays an interference pattern in the superfluid regime, and disappears
in the Mott-insulator phase. Moreover, it has been pointed out that measurements of shotnoise in
TOF images after the expansion (the results of which are different from one experimental run toanother) contain information about the quantum statistics of the particles and can help to further
characterize the quantum phases. These ideas were first theoretically developed by Altman et al.
in 2004 [21].

19

2.2

Coherence, interference and correlations

2.2.1

Coherence

A Bose-Einstein condensate (BEC) is a state of matter that emerges in a degenerate Bosegas, and it is usually described as a macroscopic matter-wave field. The first study of BEC goesback to nighty years ago when it was first predicted for photons by S. Bose in 1924 and for atoms
by A. Einstein in 1925. The experimental realization of BEC in trapped atomic clouds [22, 23]
opened up the exploration of quantum phenomena on a macroscopic scale.
In atomic BECs, most atoms occupy the same ground state, so quantum fluctuations can beneglected to a good approximation. Therefore the field operator can be replaced by a c-number,Ψ̂(~x) → Ψ(~x). The function Ψ(~x) is often called the “condensate wave function” or “order parameter”, and it evolves according to the time-dependent GPE [24] as
i~


∂Ψ
~2 2
=−
∇ Ψ + Vext + g|Ψ|2 Ψ.
∂t
2M

(2.1)

2

Here Vext is the external potential, and g = 4π~M as with M the atomic mass and as thescattering length.
Phase coherence is one of the most essential features exhibited by BECs. Assuming a condensate confined by a box with infinitely hard walls, the wave function vanishes at the wall, and in
the interior of the box the wave function has its bulk value. The distance ξ over which the wave
function arises from 0 to its bulk values can be estimated from the GP equation, since far awayfrom the wall the wave function is governed by competitions between the kinetic energy ∼ ~2 /2M ξ 2and the interaction ∼ nU (U as in Eqn. (1.39)). These two energies are equal when
~2= nU,
2M ξ 2(2.2)

1or ξ 2 = 8πna
. ξ is referred to as the coherence length. Since it describes the distance over whichsthe wave function tends to its bulk value when subjected to a localized perturbation, it is also
referred to as the healing length [24].

20
2.2.2

Interference

Interference is a generic property of coherent waves. One of the simplest, yet most profound,interference experiments was implemented by Thomas Young in 1803: The Young’s double-slit
experiment. In the double-slit setup, a monochromatic light source with wavelength λ goes throughtwo small closely spaced holes, separated by a distance d. Two split waves interfere and form aninterference pattern on the screen. The resulting interference-pattern is composed of bright bands
due to constructive interference, interlaced with dark bands due to destructive interference.
A powerful generalization of Young’s interferometer is the diffraction grating. A grating is
an array of N regularly space lines that act as N slits. The diffracted light-signal through a grid
results in a constructive interference-pattern with N 2 enhancement in intensity. This device iswidely used from spectroscopy to telecommunications. For example, regular arrays of atoms in
crystals provide as a natural grating for X-rays and generate a diffraction pattern known as Bragg
diffraction, which are used to study the structure of crystalline solids.
In quantum mechanics, all objects have wave-like properties. For instance, electrons can beused instead of light in the Young’s double-slit experiment. When a large number of electrons isaccumulated, something like regular fringes begin to appear in the detecting plane. Historically,
the wave nature of particles was first put forward by De Broglie in 1924 and confirmed three years
later for electrons with the observation of electron diffraction independently by Clinton Davisson
and George Thomson.
A revolutionary step towards the understanding of macroscopic matter-wave arose from theexperimental-realization of BEC in atomic vapors. BECs along with quantum degenerate Fermi
gases, have provided a new platform and scope to the subject of interferometry of matter waves.
For example, interference between two BECs were observed [25], demonstrating that BECs are“laser-like” matter-waves.
In the discussion so far, interference reflects first-order correlations associated with the amplitude of two coherent waves. Interestingly, interference phenomena can also be observed with

21
incoherent sources provided one probes the effects beyond first-order correlations. In other words,
many-particle correlations can exhibit interference effects even when single-particle correlations do
not. This was first discovered by Hanbury Brown and Twiss more than fifty years ago [26]. The
HBT experiment shows that coherence between sources is not a prerequisite for interference effects.
The simplest picture of HBT interferometry, which illustrates the basic idea behind HBT [26],is to consider two distant point-sources of light, a and b, of the same frequency and with random~ apart. One measures the light falling in two independent detectors, separated
phase at a distance R
by a distance d, as shown in Fig. 2.1. The distance L between the source and the detector is muchlarger than R or d. Treating electromagnetic radiations as classical waves, the radiation field fromtwo sources at point ~r isABeik|~r−~ra |+iφa +eik|~r−~rb |+iφb ,|~r − ~ra ||~r − ~rb |(2.3)

with φa and φb being the random and uncorrelated phases. After averaging the phases, one obtains~ apart to be
the intensity-correlation function from two sources at a distance R
hI1 I2 i = hI1 ihI2 i +2 2 2|A ||B | cos [k(r1a − r2a − r1b + r2b )] ,
L2

(2.4)

where r1a (r1b ) is the distance from source a (b) to detector D1. Same for r2a (r2b ).
The above equation shows a signal that varies as a function of detector separation d on acharacteristic length-scale d = λ/θ, where λ is the wavelength of the light, and θ = R/L is theangular size of the sources as seen from the detectors. Thus, by varying the separation of thedetectors, one learns the apparent angle between the two sources, and with a knowledge of theindividual wave-vectors, the physical size of the source.
To summarize, while the double slit with a single detector is an example of amplitude interferometry, HBT interferometry is an example of intensity interferometry that measures the correlation
function hI1 I2 i of two intensities at two different detectors at a given time.

22
R	
  
a	
  

1	
  

d	
  

b	
  

2	
  

Figure 2.1: Schematic of the Hanbury-Brown and Twiss experiment, showing two independent
detectors 1 and 2 that record the intensities of radiations from two sources a and b simultaneously.

2.2.3

Correlations

We now give the definition of two important correlation-functions in cold atoms and discuss
their measurements and applications. The correlation functions arehn(~x)i = hΨ̂† (~x)Ψ̂(~x)i, andG(2) (~x, ~x0 ) ≡ hΨ̂† (~x)Ψ̂† (~x0 )Ψ̂(~x0 )Ψ̂(~x)i,

(2.5)
(2.6)

where Ψ̂(~x) is the wave-field operator, hn(~x)i is a first-order correlation-function, and G(2) (~x, ~x0 )is a second-order correlation- function. As we discussed above, the first-order correlation-functions
can be measured by using amplitude interferometry, while G(2) (~x, ~x0 ) can be measured by HBTinterferometry. In cold atoms, both interferometries take advantage of the time-of-flight technique.

2.3

Time-of-flight technique
Consider a cloud of atoms loaded in an optical lattice. It is of interest to probe the in situ

correlations of this system. However, imaging limitations make it difficult to resolve two latticesites spaced by less than one µm apart. Although recently there has been great progress in thedevelopment of high-resolution imaging techniques that could resolve single lattice-sites [27, 28, 29],

23
until now most cold-atom experiments still probe the system by using time-of-flight images taken on
the expanded atomic-cloud after turning off the lattice and other external confinement-potentials.
The expansion itself is a rather complex process. In general, interactions between particles
during the expansion can generate additional noise and correlations in the system. However, since
the system’s density sufficiently decreases after a short time of expansion, to a good approximation,
one can neglect collisional processes and assume that upon release the atoms freely evolve.
The key point behind TOF images relies on the fact that the spatial distribution after expansion reflects the momentum distribution of the atoms in the trap. There are two prerequisites
for this to be true. First, the time of flight must be long enough, so that the initial size of the atom
cloud in the trap can be neglected. Second, during the expansion process, interactions between
atoms do not modify the initial momentum-distribution. For an atom of mass M moving with anaverage velocity x/t, it’s characteristic momentum is ~Q(x) = M x/t. Appendix A discusses detailsof the TOF images for a single particle.
We first consider a BEC released from a trapping potential. After a time of flight t, theresulting density distribution yields a three-dimensional density-distribution n(~x). If interactions
can be neglected during the TOF, the average density-distribution is related to the in-trap quantum
state via the relation [1]â(~x)TOF ≈

M
~t

3/2
â(~q)trap ,

(2.7)

~t
where ~q and ~x are related by the ballistic-expansion condition ~x = M~q, and â(~q) is the field operator

for destroying a particle with momentum ~q.
Now, assume that initially a BEC is trapped in a 3D simple-cubic optical-lattice. The field~ of localized Wannieroperator Ψ̂(~r) can be expressed as a sum over the destruction operator â(R)P~ ~ . Therefore â(~q)trap is given bystates by using the expansion Ψ̂(~r) = R~ w(~r − R)âRZâ(~q)trap =e−i~q·~r Ψ̂(~r)d3~r = w̃(~q)

X~

e−i~q·R âR~ ,

(2.8)

~
R

where w̃(~q) is the Fourier transform of the Wannier function w(~r). The density distribution observed

24in the TOF image is,n̂(~x)TOF =

M
~t

3

|w̃(~q)|2 G(~q),(2.9)

where the coherence properties of the many-body state are characterized by,G(~q) =

X

~ ~0
~ R
~ 0 ),
ei~q·(R−R ) G(1) (R,(2.10)

~ R
~0
R,

~ R
~ 0 ) = hâ† â ~ 0 i.
which is the Fourier transform of the one-particle density matrix G(1) (R,~ R
R

Below, we compute various correlation-functions extractable from the TOF images [30].

2.3.1

Experimental evidence of superfluid to Mott-insulator phase transition

The Superfluid to Mott-insulator transition (tuned by increasing the lattice depths) has beenprobed in a gas of ultracold bosonic atoms by using TOF images [7].
As explained in Chapter 1, bosonic atoms in optical lattices are described by the BoseHubbard model, which includes nearest-neighbor tunneling and on-site atomic interactions in a
lattice potential. In the superfluid regime, the ground-state energy is minimized if all N atoms
occupy the same identical lowest-Bloch-state (the single-particle wave function for each atom isspread out over the entire lattice with L lattice sites). The many-body ground state is then givenby (for a homogeneous system) as
N
L
X
1
1
or |ΨSF i = √  √âj  |0i,
L j=1
N!


1
|ΨSF i = √N!
where b̂†q = √1L

 N
b̂†0
|0i

(2.11)† iqxj.
j=1 âj e

PL

Then the expected value of hn̂(x)i in the TOF image is given byhn̂(x)iTOF ∝ N |w̃(q)|2 δq,G ,(2.12)

which predicts a high-contrast three-dimensional interference-pattern where sharp peaks arise at
the zero momentum and reciprocal lattice vectors. In real experiments, what is measured is a 2DR
absorption image by integrating n(~x) along one lattice direction such as dx hn̂(~x)iTOF . See Figs.
2.2 and 2.3.

25
In the Mott-insulator limit, to minimize the interaction energy, atoms tend to avoid eachother. The many-body ground state is then a product of local Fock-states (also known as number
states). For a commensurate filling of n atoms per site (n is an integer) in the homogenous case,the many-body ground state isLY
1  † n√|ΨMI i =âj |0i.
n!
j=1(2.13)

Figure 2.2: Schematic 3D interference-pattern (in the superfluid regime) with measured absorptionimages taken along x and y directions after some time-of-flight. From [7].

Then the expected value of hn̂(x)iTOF in the TOF image is given byhn̂(~x)iTOF ∝ hn̂(~q)itrap = N |w̃(~q)|2 ,(2.14)

which predicts an approximated Gaussian background and no visible interference pattern.

2.4

Noise correlations and the Hanbury-Brown-Twiss spectroscopy
Although the strongly correlated properties of the Mott insulator are hidden in the TOF

image of n(~x), they are revealed by higher-order correlations in TOF images. We now discuss how
to detect noise correlations (density-density correlations) in the expanding atomic-cloud [21]. These
correlations are characterized by [21, 1]hn̂(~x)n̂(~x0 )i = hn̂(~x)ihn̂(~x0 )ig (2) (~x, ~x0 ) + δ(~x − ~x0 )hn̂(~x).i(2.15)

26

Figure 2.3: Time-of-flight images. The figure shows the interference pattern for different latticedepths, (a) 0ER , (b) 3ER , (c) 7ER , (d) 10ER , (e) 13ER , (f) 14ER , (g) 16ER , and (h) 20ER .
As the lattice depth is increased, a quantum phase-transition from superfluid to Mott-insulator
occurs (The superfluid-Mott transition happens around 12-13 ER ). The ballistic-expansion timewas 15ms. The narrow peaks developed during the expansion demonstrate its phase coherence
across the lattice at the moment of release. Peaks appear at Reciprocal lattice vectors, with anenvelope function |w̃(~q)|2 . From [7].

27Here g (2) (~x, ~x0 ) is the normalized pair-distribution-function. δ(~x−~x0 )hn̂(~x)i is a self-correlationterm.
In a bosonic Mott-insulating or fermionic band-insulating state [31, 32, 33] in a threedimensional simple-cubic lattice, each lattice site is occupied by a fixed atom number nR~ . Such aquantum gas is released from the lattice potential and freely expands. After a flight time t, thedensity-density correlation-function is given by2hn̂(~x)n̂(~x0 )i = |w̃(M~x/~t)|2 |w̃(M~x0 /~t)|2 N 2 × 1 ±

1
N2

X

0 ~
.ei(~x−~x )·R(M/~t) nR~ (2.16)~R

The plus sign corresponds to the case of bosonic particles, and the minus sign to the caseof fermionic particles. The different sign reveals that noise correlations contain information on thequantum statistics. Equation (2.16) shows peaks (bunching) or dips (antibunching) in the density~ of thedensity correlations for bosons or fermions whenever ~q − ~q0 is a reciprocal lattice vector G
underlying lattice. The peaks (dips) become sharper with increasing N.
The different behavior exhibited by bosons and fermions can be understood in the following
way. Consider the case of two atoms localized at two sites i and j, as indicated in Fig. 2.4 (b).
There are two indistinguishable ways for particles to reach two detectors at position P 1 and P 2.
The two paths constructively interfere for identical bosonic-particles, giving rise to bunching (peaksin the correlation function), and destructively interfere for fermions, giving rise to anti-bunching(dips in the correlation function).
The bunching for a bosonic Mott insulator and the anti bunching for a fermionic band insulator have been observed experimentally [34, 35, 18]. In these experiments, several independentimages of the desired quantum-state were recorded after the atoms were released from the opticaltrap and let to expand for a finite time. These individually recorded images differ from each other
only in the shot noise. Then by analyzing the same set of images, the spatial noise correlation(2)function gex is obtained ashn(~x + ~b/2)n(~x − ~b/2)id2 ~x.
hn(~x + ~b/2)ihn(~x − ~b/2)id2 ~x

R

(2)
gexp
=R

(2.17)

28

Figure 2.4: Illustration of the HBT interferometry and the origin of quantum correlations. (a),After a time-of-flight period, the absorption image of the atomic cloud is detected by a CCD
camera. Two pixels of the camera P1 and P2 are highlighted, each of which records the atoms ina column along its line of sight. (b), When two atoms initially trapped at lattice sites i and j arereleased and detected at P1 and P2 simultaneously, the two indistinguishable paths, illustrated bysolid and dashed lines, interfere constructively for bosons or destructively for fermions. (c), The
resulting probability of simultaneously finding an atom at each detector is modulated sinusoidally
as a function of d in the unit of l (black curve). The multiple wave generalization to a regular arrayof six sources with the same spacing is shown in green. a.u., arbitrary units. From [34].

29

Figure 2.5: Noise correlations of a bosonic Mott insulator released from a 3D optical lattice. (a)Single shot image of the two-dimensional density-distribution of a Mott-insulating atomic cloud
released from a 3D optical lattice. (b) Horizontal cross section through the center of the image in(a), and Gaussian fit to the average over many independent images. (c) Spatial noise correlations
obtained by analyzing the same set of images. (d) Horizontal profile through the center of thepattern in (c). Those peaks are separated by integer multiples of l = Mhta , reflecting the underlying
lattice-structure. From [34]Figure 2.6: Noise correlation of a fermionic band insulator. Each plot has a similar meaning as inFig. 2.5. Instead of the correlation “bunching” peaks observed in Fig. 2.5, the correlation functionof fermioninc quantum-gas shows an antibunching effect. From [35].

30
(2)

As shown in Fig. 2.5 and 2.6, the Mott-insulating state exhibits long-range order in gex .
Since the contribution from the averaged density n(~x) is smoothed out in g (2) (~x, ~x0 ) [See Equation(2.15)], in the superfluid regime, one expects g (2) (~x, ~x0 ) = 1. In contrast, the Mott-insulator state
exhibits long-range order in the pair correlation function, manifesting as an array of peaks withspacing l = ht/aM in the noise correlation function.

2.5

Cold atoms in one dimension

2.5.1

Interaction strength γ

Correlation effects strongly depend on dimensionality. The experimental realization of cold
atoms confined in waveguides [5, 6] provides a way of reaching the strong-interaction regime indilute gases.
In 1963, Elliot Lieb and Werner Liniger [36] calculated the exact solutions for 1D Bose gaseswith δ(z) interactions (Lieb-Liniger gas). In 1998, Olshanii extended the discussion to atoms ina highly elongated cigar-shaped potential [37]. The interaction strength of a 1D gas is usuallycharacterized by the dimensionless parameter γ = Eint /Ekin , where Ekin is the quantum kineticenergy needed to correlate particles by localizing them with respect to each other on the order of the2~mean inter particle distance d, Ekin = 2M
. Eint = g1D n1D is the average mean-field interaction,d2

with g1D the strength of the 1D interaction, and n1D = 1/d the 1D density.
Assume that a gas of cold atoms is loaded in a strong 2D optical lattice in the x, y directions,which tightly confines the atoms in the transverse direction. In the regime where the excitationenergy ~ω⊥ of motion in the x,y direction
is larger than the chemical potential, only the lowest
2vibrational state φ0 (x, y) =along the z direction.

1

exp − x +y
2a2
⊥
√
πa⊥

2

is accessible. At the same time, atoms can move freelyp
In a realistic system, the harmonic-oscillator length a⊥ = ~/µω⊥(where µ = M/2 is the reduced mass) in the transverse direction is usually larger than the 3D
s-wave scattering-length as . Then by integrating the 3D interaction-potential over the transverse
1

Considering two atoms in a harmonic confining-potential, the atomic motion can be separated into a center-ofmass motion and a relative motion. The relative wave functions are of the form φ0 (x, y)e−ikz .

31direction, gR

d3 x|φ0 (x, y)|2 δ(~x), one can get the resulting 1D effective interaction g1 δ(x) with g1 =g|φ0 (0, 0)|2 = 2~ω⊥ as in the weak confining limit as  a⊥ . Consequently
γ=

2as
.
n1D a2⊥(2.18)

In contrast to the 3D situation, the strength of interaction γ scales inversely with the density
n1D . Therefore, in 1D, interactions dominate in the low-density limit. For γ  1, the interacting
bosons behave like non interacting fermions — the so-called Tonks Girardeau (TG) gas. γ =5.5, was reached experimentally [5] and clear signatures of the fermionization of the gas were
demonstrated. The momentum distribution of TG gas with a superimposed lattice-potential alongz was also observed for ultracold atoms [6].

2.5.2

Two-particle correlations g (2)An intuitive understanding of the gradual transition from a weakly interacting gas and TG
regimes with increasing value of γ can be gained by considering the two-particle pair correlations2g (2) (x) ≡ hn̂(x)n̂(0)i−δ(x)hn̂(x)i. It is a measure of the probability of observing two particles separated
hn̂(x)ihn̂(0)i
by x simultaneously. For an ideal BEC in three dimensions, g (2) (~x) = 1 at arbitrary ~x. For the
Lieb-Liniger gas, the local value of the pair correlation g (2) (0) can be calculated exactly [38], which√√
gives g (2) (0) = 1 − 2 γ/π + · · · in the limit of γ  1 and g (2) (0) = (2π/ 3γ)2 → 0 when γ  1.
Experimentally, g (2) (0) was determined by Kinoshita et al, [39] using photoassociation. As shownin Fig. 2.7, the experimental results are in good agreement with theory over a wide range ofinteraction strengths.

2

Those measurements were in situ, but not after time of flight.

32

Figure 2.7: Local pair correlation function g (2) ≡ g (2) (0) as a function of interaction strength γ.
The solid line shows the theory prediction from [38] . The plot is taken from [39].

Chapter 3

Quantum correlations in bichromatic lattices

Correlation functions can be used as powerful tools to characterize complex phases. In thischapter, we study the many-body Schrödinger equation in a quasi-periodic potential and calculate
the momentum distribution and noise correlations. We discuss the connection between them with
the Kolmogorov-Arnold-Moser (KAM) problem of classical mechanics. We propose a possible
visualization of such a connection in experimentally accessible many-body observables. Thoseobservables are useful probes for the three characteristic phases predicted in a bichromatic lattice:the metallic, Anderson, and band-insulator phases. In addition, we show that those correlation
functions exhibit fingerprints of nonlinear phenomena such as bifurcations and devil’s staircases.

3.1

Introduction
Systems with competing lengths have fascinated physicists as well as mathematicians becauseof their exotic fractal-characteristics [40, 41]. Such systems with two competing periodicities,
commonly known as almost-periodic or quasi-periodic, occur very commonly in nature. The mostcommonly studied example in quantum physics is the single-particle Schrödinger equation in thepresence of a quasi-periodic (QP) potential.
Motivated by the experimental realization of two-color lattices, we revisit the problem of
the Schrödinger equation in a quasi-periodic potential and its relationship with the KolmogorovArnold-Moser (KAM) [42] problem of classical mechanics. Our particular focus is the treatment
at the many-body level of those effects that have been extensively studied in QP systems at the34single-particle level. By focusing on experimentally accessible observables, such as the momentum
distribution and density-density correlations, we demonstrate the possibility of experimental visualization of the relationship between the metal-insulator transition and the KAM-Cantori transition[43] in many body systems. These observables are found to exhibit fingerprints of various paradigmsof nonlinear systems such as Arnold tongue-like structures, bifurcations, and devil’s staircases.
The many-body systems that we treat here are ensembles of ultracold spin-polarized fermionic
atoms confined in one-dimensional bichromatic optical-lattices. We treat cases in which the twolattices have incommensurate periodicities, and thus constitute a quasi-periodic potential for atomic
motion [44]. Since the metal-insulator transition in single-particle QP systems is associated witha localization transition of extended single-particle states, a natural probe of such transitions in
many-particle systems is the density distribution and density-density correlations. At present, thedensity distribution of atoms confined in a lattice is not an easily accessible experimental observable;
instead, the quasi-momentum distribution and its corresponding correlation-functions can be easily
measured after the lattice potential is suddenly removed.
Using the momentum-position duality of our basic model [45, 46], we show that TOF images
encode local-density information and can be used to identify the possible phases of ultracold atomicgases, e.g., metallic phases and various types of insulators. Many of our results are explained
analytically with perturbation theory, although the complete picture is obtained by exact numericalcalculations, particularly near the metal Anderson-insulator transition.
This chapter is organized as follows. First, we give an overview of experiments on ultracold
atoms in QP optical lattices. Then we discuss the basic Hamiltonian, and the experimental observables. After that, we describe the relationship between the metal Anderson-insulator transition
and KAM-Cantori transition, and follow it by a discussion of possible experimental manifestationsof effectively nonlinear behavior in these many-body systems. Finally, we provide a summary of
our results and state our conclusions.

35

3.2

Ultracold atoms in a 1D bichromatic optical lattice
The systems of interest here are gases of ultracold spin-polarized fermionic atoms trappedin the lowest band of a transversal 2D optical-lattice. For many systems of current experimental
interest, where atomic interactions are short-range, spin-polarized fermions are effectively noninteracting because of Pauli blocking. The 2D lattice depth is made strong enough to freeze themotion of the atoms transversally, creating an array of independent 1D tubes [5, 47]. Along theaxis of the 1D tubes, additional optical lattices can be imposed, as has been done in Refs. [48, 6].
We discuss cases in which two such lattices are imposed, with incommensurate periods [8] (whenone of these lattices is much stronger than the other, we refer to it as the primary lattice). The
combined lattices therefore generate an effective QP potential along the axial direction. Such an
experiment has been recently implemented for bosonic atoms in Ref. [49]. Our treatment considers√
cases where the ratio of the two lattice constants is equal to the “golden mean”, σ = ( 5 − 1)/2,
which is one of the best-studied examples in QP systems [46].
In most experiments, atoms are first trapped and cooled to quantum degeneracy. They aresubsequently loaded into the lattice by adiabatically turning on the lattice laser beams. At the endof each experimental sequence, atoms are probed by using TOF images. These are obtained after
releasing the atoms by turning off all the confinement potentials. The atomic cloud expands and isthen photographed after it enters the ballistic regime. Assuming that the atoms are noninteractingfrom the time of release, properties of the initial state can be inferred from the spatial images[21]: the column density distribution image reflects the initial quasi-momentum distribution, andthe density fluctuations, namely the noise correlations, reflect the quasi-momentum fluctuations.
These quantities, which will be defined below —see Eqns. (3.5,3.6)—, have been shown to be
successful diagnostic tools for characterizing quantum phases and have been recently measured inbosonic quasi-periodic systems [50].

36

3.3

Our model system: The Harper equation, many-body observables andself-duality

3.3.1

The Harper equation

If the intensity of the secondary lattice is much weaker than that of the primary lattice, the low
energy physics of the fermionic system can be well described by the tight-binding Hamiltonian[8]:H = −JX(ĉ†j ĉj+1 + ĉ†j+1 ĉj ) +jX

Vj n̂j ,(3.1)

jwhere ĉj is the fermionic annihilation operator at the lattice site j, and J is the hopping energy
between adjacent sites. The main effect of the QP potential is to modulate the on-site potential.
It is accounted for by the term Vj = 2V0 cos(2πσj + φ). The parameter V0 is proportional to theintensity of the lasers used to create the secondary lattice [44], σ is the ratio between the wave√
vectors of the two lattices that we choose to take the value σ = ( 5 − 1)/2, and φ is a phase factor
that is determined by the absolute registration of the two lattices.
To model real experimental-conditions, averaging over φ is required. This averaging takes
into account, on one hand, the phase fluctuations from one preparation to another. Those arise
because of the difficulty of locking the position of the cloud over several shots. On the other hand,the fact that typical experimental setups generally consist of an assembly of one-dimensional tubes
with different lengths and phases among them.
(m)

For a single atom, the eigenfunctions ψjand eigenenergies E (m) of the Hamiltonian in Eqn.

(3.1) satisfy:


(m)
(m)
(m)
(m)
− ψj+1 + ψj−1 + 2λ cos(2πσj + φ)ψj = (m) ψj ,(3.2)

where λ = V0 /J, (m) = E (m) /J, and (m) ≤ (m+1) . Equation (3.2) is known as the Harper
equation, a paradigm in the study of 1D quasi-periodic (QP) systems [46]. For irrational σ, themodel exhibits a transition from extended-to-localized states at λc = 1. Below criticality, all thestates are extended Bloch-like states characteristic of a periodic potential. Above criticality, theHarper model becomes equivalent to a corresponding Anderson model, the spectrum is a pure point

37
spectrum, and all states are exponentially localized. At criticality, the spectrum is a Cantor set
and the gaps form a devil’s staircase of measure unity [40].
In our numerical studies, σ is approximated by the ratio of two Fibonacci-numbers FM −1 /FM ,(F1 = F0 = 1, Fi+1 = Fi + Fi−1 ) that describe the best rational approximant by continued fractionexpansion of σ. For this rational approximation, the unit cell has length FM , and the single-particle
spectrum consists of FM bands and FM − 1 gaps. The gaps occur at Qn /2, ±(π − |Qn |/2) with
Qn = ±(2π/a)hnσi reciprocal lattice vectors constrained in the interval aQn ∈ (−π, π]. Herehnσi = nσ (mod 1), and n an integer. We denote by Np the number of atoms, Nl = FM is thenumber of lattice sites, and the filling factor ν = Np /Nl ranges from 0 to 1.

3.3.2

Many-body observables

An ensemble of spin-polarized fermions at zero temperature are “stacked up” into the singleparticle eigenstates of increasing energy, with one particle per quantum-state. The energy of thehighest occupied-level, which depends on the filling factor ν, is the Fermi energy, EF . Since at
the critical point all the single-particle wave-functions become localized, at the many-body level
polarized fermions also exhibit a transition from metal to insulator at λc . However, in addition tothese two phases, the fragmentation of the single-particle spectrum in a series of bands and gaps
introduces additional band-insulator phases when the Fermi energy lies within a gap. The most
relevant insulating phases occur at the irrational filling factors: ν = FM −1 /FM and ν = FM −2 /FM(which respectively correspond to ν = σ, 1 − σ), associated with the leading band gaps. In theband insulator phases the many body system is an insulator, irrespective of the value of λ, as theextended-single-particle states are occupied.
Since the metal-insulator transition is clearly signaled by the onset of localization of extended
single-particle states, a natural probe of this transition is the many-body density profile, ρj , andthe density-density correlations, ∆(j1 , j2 ), which can be written in terms of single particle wave

38
functions as:ρj= hn̂j i =

Np
X

(m) 2

ψj

, and

(3.3)

m=1

∆(j1 , j2 ) = hn̂j1 n̂j2 i − hn̂j1 ihn̂j2 i =

Np
X

(m) 2

ψj1

δj1 ,j2 −NpX2(m) ∗ (m)
ψj1 ψj2

.

(3.4)

m=1

m=1

Here n̂j = ĉ†j ĉj , and we have used Wick’s theorem to evaluate ∆(j1 , j2 ). However, in general,such local observables are hard to measure experimentally because of the lack of addressability of
individual lattice- sites for typical lattice-spacing. Instead, TOF images access nonlocal observables
such as the quasi-momentum distribution, n̂(Q) and the quasi-momentum fluctuations, ∆(Q, Q0 ),which are given byn̂(Q) =

1 X iQa(i−j) †eĉi ĉj , and
Nl

(3.5)

i,j

∆(Q, Q0 ) = hn̂(Q)n̂(Q0 )i − hn̂(Q)ihn̂(Q0 )i,(3.6)

where Q is the quasi momentum that can be expressed in terms of the index k (k is an integer) ask
Q(k) = 2π
a Nl . Q(k) is constrained to the interval Qa ∈ (−π, π] .

Introducing ηk , the Fourier transform of ψk ,(m)

ηk

1 X −i 2πkj(m)=√e Nl ψj ,
Nl j

(3.7)

where the observables nk (also denoted sometimes as n[Q(k)]) and ∆(Q) ≡ ∆(Q, 0) can be written
as
nk =

Np
X

(m) ∗ (m)
ηk

ηk

,

(3.8)

m=1

∆[Q(k)] =

Np
X
m=1

3.3.3

(m) 2
η0
δk,0 −Np
X2(m) ∗ (m)
ηk η0

.

(3.9)

m=1

Self-duality

The self-duality of the Harper equation corresponds to the property that single-particle eigenP
states ψj and their corresponding Fourier-transformed eigenstates, zn ≡ √1N e−inφ ψj e−ij(2πσn+θ)l

39
satisfy the same equation with the roles of J and V0 interchanged [45, 46], i.e.,− (zn−1 + zn+1 ) +2
cos(2πσn + θ)zn = − zn ..
λλ(3.10)The self-duality of the Harper equation implies that if ψ is an extended-eigenstate, its selfdual eigenstate z is localized and vice-versa. ηk is connected to zn by the relation zn = ηk e−inφwith k = Mod[nFM −1 , FM ] when θ is set to 0.
Thouless has developed a formula [51] that relates the exponent for the decay of a localizedstate to the integrated density of states. Denoting γ1 and γ2 , respectively, as the decay exponentof the localized states for ψ and z, the Thouless relation can be written as,2(E − E 0 )ρ(E 0 ).
JZ
2(E − E 0 )γ2 (E) =dE 0 lnρ(E 0 ).
V0Z

γ1 (E) =dE 0 ln(3.11)(3.12)

The above equations imply that,
γ1 (E) = γ2 (E) + ln

V0
J


.

(3.13)

Since γi (i = 1, 2) is related to the integrated density of states, γi ≥ 0. Consequently, if ψ islocalized, then z must be extended, and therefore γ2 must be zero and γ1 = ln(λ). The exponentγ, also known as the Lyapunov exponent, describes the exponential growth of the renormalized
couplings in the renormalization-group equations [52, 53].
Self-duality is a key to obtaining experimental information on local quantities from measurements. In this context, it is important to understand the relationship between the index j of ψjthat satisfies the Harper equation (3.2) and its corresponding index k in ηk that satisfies the dual
equation (B.1). We obtain this relationship numerically by diagonalizing the Harper equation fora given λ, labeling the states in increasing order in energy, and then finding the corresponding
momentum-space dual by repeating the same procedure, but with λ replaced by 1/λ. An explicit example of the mapping is provided in Appendix B. Note that Fibonacci sites in real space
are mapped to Fibonacci sites in momentum space up to a common displacement. This shift isdependent on the phase factor φ.

40

3.4

Localization transition as a KAM-Cantori transition
As discussed by Ostlund et al. [52, 53], there is a close parallel between the KAM problem

of classical mechanics and the quantum problem of a particle in a QP potential.
Firstly, both systems share the mathematical difficulty of having higher-order terms withsmall denominators when the QP potential in the quantum problem or the nonintegrable term in
the classical problem are treated perturbatively. In the classical system, Kolmogorov, Arnold andMoser (KAM) [42] solved this “small divisor” problem and demonstrated that most of the invariant
tori in the phase space are not destroyed by a sufficiently weak nonintegrable perturbation, andthe phase space trajectories of these weakly perturbed nonintegrable systems are topologicallyequivalent to those of an integrable system. Outside the perturbative regime, extensive studies
have shown that the KAM tori break, becoming an invariant cantor-set, known as Cantori [43].
The transition from KAM tori to Cantori has been studied extensively using numerical and by
nonperturbative renormalization-group (RG) methods both in the dissipative systems, such as
circle map [54], and in Hamiltonian systems such as the area-preserving standard maps [55]. The
KAM-Cantori transitions resemble second-order phase transitions of statistical mechanics because
the critical point can be described by a scaling theory with universal critical exponents.
In the context of QP potential, Dinaburg and Sinai [56] showed that for a sufficiently weak
QP potential, a large part of the spectrum remains continuous, and the wave functions exhibit
a generalized Bloch form of ψn = e2πinκ χ(nσ), where χ(x) is smooth and periodic. Outside theperturbative regime, QP systems have been studied numerically and using various RG schemes.
As discussed by Ostlund et al. [53], the RG methods used to describe the critical point of the
QP trajectories in the circle map are closely related to the RG methods used to study scaling
properties of the spectrum and the wave function of the discrete Schrödinger equation with QP
potential. Numerical studies [52] have demonstrated a close parallel between the KAM tori and thesmooth Hull-function of the extended phase and the Cantori and the discontinuous Hull-functionof the localized phase. Similar to the critical torus, the critical Hull-functions are self-similar and41
are characterized by universal scaling exponents.
The correspondence between the KAM-type problem of classical physics and the quantum
problem of a particle in a QP potential (as described above), has been visualized within a singleparticle picture. Below, we elucidate this correspondence in the many body context, which is
directly relevant for current cold-atom experiments. In other words, instead of looking at thesingle-particle wave function, whose Hull function, χ(nσ), is a smooth and continuous function in
the extended phase but becomes discontinuous for λ > 1, we propose to look at the return mapof the local density of the atomic cloud, ρj vs ρj+1 . This type of return map provides a cleaner
visualization of the KAM-to-Cantori transition at the many-body level.
Figure 3.1 shows the return map for various rational filling factors ν = Np /Nl = 1/4, 1/2 anddifferent values of λ. For λ < 1, the return maps are smooth curves and correspond to the KAM
tori of the extended or metallic phase. For λ > 1, the density profile is a discontinuous function, aCantori. The discreteness of the return map for generic filling factors can be easily understood inthe λ → ∞ limit, where the wave functions are fully localized, and thus the return map can only
take the four possible values of (0, 0), (0, 1), (1, 0), (1, 1). Exactly at the transition point λ = 1, thesmooth curves become disconnected.
Using the duality transformation, similar return maps can be drawn in momentum space,nk vs nk±k0 (k0 = FM −1 ). While nearest-neighbor sites are connected in position space becauseof a finite J, quasi-momentum components separated by the main reciprocal-lattice vectors of
the secondary lattice are connected by a finite V0 . We find the momentum return-maps exhibit
an important advantage compared with the local-density maps related to the fact that they areinsensitive to variations of the phase and retain their pattern when averaged over phase. This is
not the case in the density maps shown in the insets of Fig. 3.1.
We now consider the case where the filling factor is a ratio of two main Fibonacci-numbers(ν = σ, 1 − σ). We refer to these filling factors as irrational filling factors since they approachan irrational number in the thermodynamical limit. As discussed earlier, at irrational filling, thesystem is in a band-insulating phase as the Fermi energy EF lies in the gap, just outside a filled

42

(a)

(b)

(c)

(d)

Figure 3.1: Return maps in real space for (a) and (c) λ = 10(red), 1(green), 0.5(blue) and thecorresponding return map in momentum space (b) and (d) λ = 0.1(red), 1(green), 2(blue). The
filling factors are ν = 0.25 for the upper panels, and ν = 0.5 for the lower panels. The insets areaveraged over 50 random phases.

band.
In contrast to the generic or rational filling factors discussed earlier, the return map forthe irrational-filling case (ν = 1 − σ) is found to remain smooth regardless of the value of λ(See Fig. 3.2). This result may appear somewhat counterintuitive because as λ → ∞, all thesingle-particle wave functions become localized. However, an exception to this simple picture
occurs near irrational filling when the density function exhibits a continuous distribution. Thisdifference between rational and irrational filling is due to the presence of a group of paired states
centered around the most dominant band edges (associated with the dominant gaps). Each of

43

(a)

(b)

Figure 3.2: Return maps for the irrational filling factor ν = 1 − σ. Panel (a) is in real spacewith λ = 10(red), 1(green), 0.5(blue) and panel (b) is in momentum space with the dual values
λ = 0.1(red), 1(green), 2(blue). The insets are averaged over 50 random phases.

44
those states is dimerized, in other words we mean that they localize at two neighboring sites. The
existence of a pair of such dimerized states is the key to understanding the difference between
rational and irrational filling. Although a single dimerized-state causes delocalization near bandfilling (explaining irrational filling case), the existence of a pair can cancel the delocalization effect.
Technical details of this argument are presented in Appendix C, where we show that the smoothcharacter of the map at the irrational fillings, σ, 1 − σ associated with the leading gaps can beunderstood by the breakdown of nondegenerate perturbation theory in the large λ limit.
In summary, the band insulating phase belongs to the KAM phase irrespective of the strengthof the disorder.

3.5

Fingerprints of nonlinear phenomena in many-body observables
We now study the quasi-momentum distribution of polarized fermions, that is directly acces-

sible in the TOF images of ultracold atoms. In the following, we will show how the quasi-momentum
distribution imprints a signature of the KAM-Cantori transition and provides experimental realization of various landmarks of nonlinear systems.

3.5.1

Fragmented Fermi sea

Fermions in the extended phase have metallic properties. For λ = 0, the single particle
eigenstates are fully localized in quasi-momentum space (and thus delocalized in position space).
The quasi momentum distribution n(Q) has a step-like profile: n(Q) = 1 for |Q| 6 QF and n(Q) = 0
for |Q| > QF (QF is the Fermi momentum). For 0 < λ < 1, the single particle eigenstates localized
at k acquire some admixture of other quasi-momentum components ( nk±FM −1 to leading orderin λ). In this regime, the quasi-momentum distribution retains part of the step-like profile but
gets fragmented into additional structures centered at different reciprocal lattice vectors of thesecondary lattice. We call the filled states centered around Q = 0 the main Fermi sea and thosearound the QP-related reciprocal lattice vectors the quasi-Fermi seas.
The fragmentation of the quasi-momentum distribution and the development of QP Fermi

45

Ν
1.0
0.8

n 6
n 2

n 3

0.6
0.4 n 4
0.2
-Π

- Π2

n 1

n 5
Π
2

Π

Qa

Figure 3.3: Various boundaries of the fragmented Fermi sea.

seas can be understood from first-order perturbation theory. To leading order in λ, the perturbation tools discussed in Appendix C can be used to show that the momentum landscape becomesfragmented in six regions, as shown in Fig. 3.3 and given by: 2
λ1
ñ1 (Qa) =,22
2
4 sin Θ sin (|Qa| − Θ)

where Θ =(3.14)ñ2 (Qa) = 1 − ñ1 (π − |Qa|),
 2


λ
1
1
1
ñ3 (Qa) =+,2
4 sin2 Θ sin2 (|Qa| + Θ) sin2 (|Qa| − Θ)(3.15)ñ4 (Qa) = 1 − ñ3 (π − |Qa|),(3.17)ñ5 (Qa) = 0, and(3.18)ñ6 (Qa) = 1,

(3.19)(3.16)FM −2|Q1 a|FM π =
2 . The various regions are delineated by the edges of the QP Fermi seas,

46
shown by white lines in Fig. 3.3 and given in the Q vs ν space byν=aa|Q(k ± σNl )| = |Q(k ± FM −1 )|.
π
π(3.20)The regions ñ1,2 , ñ3,4 and ñ5,6 reflect the characteristic particle-hole symmetry in fermionicsystems.
As disorder is increased, more quasi-Fermi seas become visible, and with the increase of ν they
overlap in a complicated pattern. Exactly at criticality, λ = 1, the pattern evolves into a fractal-likestructure. Beyond the critical point, λ > 1, the fragmented quasi-momentum distribution profile
disappears and instead becomes a smooth function of Q and ν. This behavior is summarized inFig. 3.4, and the corresponding cross sections for fixed Q and ν are displayed in Figs. 3.5 and 3.6.
At irrational filling factors, however, the quasi-momentum distribution profile remains smoothregardless of the value of λ. The smoothness of the momentum distribution at these special filling
factors can be understood using the same reasoning as the one used to understand the smoothcharacter of the return map in position space in terms of nondegenerate perturbation theory. In
this case, however, the states that are coupled are the ones localized at the quasi-momenta QF + and −QF + , (a  1).
The fragmentation of the Fermi sea in the λ < 1 regime, and the smooth profile in both theλ > 1 and band-insulator phases are directly connected via the self-duality property to the discrete
nature of the density return map in the localized phase, and its smooth character in the extended
and band-insulator phases.
An important point to emphasize again, which is crucial for possible experimental observationof the predicted behavior, is the insensitivity of the momentum distribution to variations of thephase φ. We demonstrate such insensitivity by noticing that the quasi-momentum distributions
plotted in Fig. 3.4 are actually averaged over many different values of φ. Hereafter, all the plots inmomentum space are always averaged over 50 random phases.
The emergence of a tongue-like structure (see Fig. 3.7) as the QP potential increases andits disappearance beyond criticality is reminiscent of the mode-locked windows such as Arnold

47

(a)

(b)

(c)

nHQL

Figure 3.4: Momentum distribution for (a)λ = 0.1, (b)λ = 1, (c)λ = 2. All figures are averaged
over 50 random phases to mimic a realistic experimental case. We can see that the sharp featuresof these distributions survive phase averaging.

1.0
0.8
0.6
0.4
0.2
0.0
0.0

Λ=0.1
Λ=1
Λ=2

0.2

0.4

0.6

0.8

1.0

Ν
Figure 3.5: Momentum distribution for fixed Q(k = 1045) , FM = 4181 as a function of filling
factors vs. the disorder parameter.

1.0
0.8
0.6
0.4
0.2
0.0
-Π

Λ=0.1
Λ=1
Λ=2

- Π2

0
Qa
(a)

Figure 3.6:
FM −1 /FM .

nHQL

nHQL

48

Π
2

Π

1.0
0.8
0.6
0.4
0.2
0.0
-Π

Λ=0.1
Λ=1
Λ=2

- Π2

0
Qa

Π
2

Π

(b)

Momentum distribution vs. the disorder parameter. (a) ν = 0.5, and (b) ν =

(a)

(b)

Figure 3.7: Momentum distributions as a function of the disorder parameter. (a) filling factor
ν = 0.5. Here we see the formation of two tongues that merge at the critical value λ = 1. (b)ν = σ. When ν is irrational, the system is a band insulator and the Arnold tongue-like structure
disappears. We have checked for the existence of the similar structure in position space whenconverting λ to 1/λ.

49
tongues in QP circle maps [57]. For λ = 0, the Fermi distribution is a binary distribution with
only two values n(Q) = 0, 1. For finite λ and for generic filling factors, two windows of values
centered around zero and unity become allowed. Their width increases with λ as they can bederived from a perturbative analysis. We call such windows QP tongues. Exactly at the onset ofthe metal-insulator transition, the two distributions overlap, mimicking the Arnold tongue overlap
at the onset to transition. We have checked that the duality of the Harper equation allows one toobserve the formation of analogous structures in position space. The absence of a metal-insulatortransition at the irrational filling factors at which the system is a band insulator is also signaled by
the disappearance of the tongues at these fillings, as shown in Fig. 3.7(b).

3.5.2

Bifurcations

Bifurcations are common features observed in nonlinear dynamical systems that occur whena small smooth change made to specific parameters (the bifurcation parameters) of the systemcauses a sudden ‘qualitative’ or topological change in system behavior. A series of bifurcations canlead the system from order to chaos.
The density distribution provides a nice manifestation of a single bifurcation with λ as abifurcation parameter. When the density at the various Fibonacci sites is plotted as a function of
λ, for a specific filling factor that depends on the value of φ, a bifurcation opens up at λ = 1. In
Fig. 3.8, the existence of a bifurcation at quarter filling (ν = 1/4) when φ = 3π/4 is shown.
A qualitative understanding of the bifurcation can be gained by considering the two limiting
regimes, λ  1 and λ  1. In the weak-coupling limit, λ  1, the local density is uniform anddirectly proportional to the filling factor ν. In contrast, in the strong coupling limit, a given siteremains empty or occupied depending upon whether the on-site potential is greater than or lessthan the Fermi energy.
Fibonacci sites j = Fn have similar on-site energies that oscillate about c ≡ 2λ cos(φ)

50Figure 3.8:
Bifurcation of local density at Fibonacci sites. The central curve (red triangle)
corresponds to quarter filling while the lowest (black blocks) and the topmost (blue disks) curves,
respectively, correspond to ν = 1/5, and ν = 1/3. Quarter filling is a special case where the local
density at Fibonacci sites is 0 or 1 as λ → ∞. For ν > 0.25, the Fibonacci sites are filled, whilefor ν < 0.25 Fibonacci sites are empty as λ → ∞. We have checked for the existence of the same
bifurcation phenomena in the quasi-momentum distribution, but with the weak and the strong
coupling limits reversed ( λ → 1/λ).

(assuming n ≈ M and φ > 2πFM −n /FM ) as

2λ cos(2πσFn + φ) = 2λ cos 2π(−1)n−1 FM −n /FM + φ  < cn is odd > cn is even.

(3.21)

These oscillations can be seen in Fig. B.1 in Appendix B.
Consequently, at the filling factor, νc , at which the Fermi energy matches c (e.g., if φ = 3π/4,
then νc ∼ 0.25) Equation (3.21) implies that as λ goes to infinity, even Fibonacci sites (F2l )become empty and odd Fibonacci sites (F2l+1 ) become occupied. This behavior, combined withthe monotonic increase of the density with λ in the weak coupling limit, qualitatively explains theobserved bifurcation.
The self-duality behavior described by Eq. (3.10) implies the existence of the same bifurcation
phenomena in quasi-momentum distribution, but with the weak and the strong coupling limits

51

(a)

(b)

Figure 3.9: Local density distribution as the disorder parameter and filling factor are varied for (a):Fibonacci site j = 377 and (b): normal site j = 100. We have checked the self-duality behavior in
the quasi-momentum distribution with the weak and the strong-coupling limits reversed (λ → 1/λ).

reversed ( λ → 1/λ).
When the local density is plotted for a generic lattice site as a function of λ and ν, thelandscape shows a similar change in topology. In contrast to Fig. 3.8 where the filling factor isfixed but different curves corresponding to different Fibonacci sites are shown, in Fig. 3.9 thedensity at a specific lattice site is plotted, but both the filling factor and disorder strength arevaried. In this plot, besides the large bifurcation, there are smaller ones which cannot be explainedby studying the two limiting cases. However, in momentum space they are qualitatively understoodby considering cuts at a fixed λ as a function of ν. In these cuts the big jump at small λ can be
identified with the edges of the main Fermi sea in the quasi-momentum profile and the smaller
jumps correspond to edges of the quasi-Fermi seas.

3.5.3

Devil’s staircases

A Devil’s staircase describes a self similar function f (x), with a hierarchy of jumps or steps.
In other words, f (x) exhibits more and more steps as one views the function at a smaller andsmaller length scale. The derivative of f (x) vanishes almost everywhere, meaning there exists a set

52
of points of measure 0 such that for all x outside it, the derivative of f (x) exists and is zero.
Cold atom experiments may provide an opportunity to visualize Devil’s staircases in momentummomentum correlations, namely the noise correlations.
Measurement of noise correlations is an example of Hanbury-Brown-Twiss interferometry(HBTI), which is sensitive to intrinsic quantum noise in intensity correlations. HBTI is emerging
as one of the most important tools for providing information beyond that offered by standard
momentum distribution-based characterization of phase coherence. The noise correlation patternin 1D QP bosonic systems has been studied theoretically [58, 59, 60, 61] and has also been measured
experimentally [50]. Here, however, we focus on the fermionic system.
In the extended phase, noise correlations exhibit a series of plateaus as ν is varied, and thenumber of steps or plateaus increases as the strength of the disorder λ increases. The origin of this(n)step-like structure with jumps occurring at the filling factors νju = |Qn |a/π can be understoodfrom a perturbative argument as follows. For noninteracting fermions, for Q 6= 0 (See Fig. 3.10),PNp(m) (m) 2(m)∆(Q(k)) = − m=1
ηk η0, with ηk being the Fourier transform of the mth single-particle
eigenfunction. For λ = 0, the overlap between any two different Fourier components is always zero(1)

as only the ground (m = 1) state has a zero quasi-momentum component, i.e., η0 = 1. For small(1)

λ  1, first-order perturbation theory yields a single step observed at νju = |Q1 |a/π, as only(m(1) )

η0(1)

with m(1) = νju Nl is nonzero.
(1)

Quantitatively, from first-order perturbation theory, the heights of the steps at νju for λ  1
are given by, 2
λ1
ζ(0) = −2,
2(1 − cos(2πσ))2 2
λ1
ζ(Q1 ) =.
2(1 − cos(2πσ))2(3.22)(1)

where ζ(0) (ζ(Q1 )) is the step height for ∆(0) [∆(Q1 )] at filling factor νju . The minus sign implies
a decrease in the noise as the filling factor ν increases. As λ increases, more and more steps are seen
and can be explained using higher-order perturbation theory. At criticality (see Figs. 3.10–3.11),
the steps acquire a hierarchical structure that resembles a Devil’s staircase and correlates with the

53

(a)

(b)

(c)

Figure 3.10: Noise correlation for (a)λ = 0.1, (b)λ = 1, (c)λ = 2. All plots display averaged
quantities over 50 random phases. We can see that the structure in the noise correlations survives
phase averaging, as they do in the momentum distributions of Fig. 3.4. ∆(Q = 0) is not displayedin those plots.

54

HDH0L-0.16L104
5
4
3
2
1
0
0.24 0.26

DH0L

0.20
0.15

6

0.28

Ν
0.3

0.10

4

0.06

2
Ν
0.42 0.44 0.46 0.48

0.04
0.02

0.05
0.00
0.0

H-DHQ1L-0.0506L105
8

0.08-DHQ1L

0.25

0.2

0.4

0.6

0.8

1.0

0.00
0.0

0.2

0.4

Ν
(a)

0.6

0.8

1.0

Ν
(b)

Figure 3.11: This figure shows the evolution of ∆(Qn ) ( n = 0, 1) as a function of the filling factor
for different λ values. In the limit λ  1, ∆(Qn ) exhibits steps occurring at the filling factors(m)νju (m = 1, 2, . . . ) (see text). In the limit λ  1, it acquires a sinusoidal profile. The inset shows
the self-similar nature of the steps at λ = 1.

55
fractal structure of the energy spectrum [52, 53]. In contrast to the momentum distribution, noisecorrelations do not show significant differences between the rational and irrational filling factors.
In the localized phase, λ > 1, we see the smoothing of the step-structure and the noise2(πνn)
correlation function tends towards a sinusoidal profile ∆(Qn ) ≈ νδn,0 − sin(πn).
2To emphasize the striking difference between the noise patterns in the extended and localized
phases (the step-like vs. smooth profile), in Fig. 3.12, we plot a measure of the gaps developed in∆(0) as a function of λ.

1.0

ò
à
àæ
à
àæ

Measure

0.8
0.4
0.2
à
æ
ò

0.0æàò

à
æ
ò

à
æ
ò

0.4

à
æ
ò

ò
æ
à

ò
æ
à

ò
æ
à

ò
æ
à

ò
æ
à

à

M=15
M=17
M=19

0.6

ò
æ
à

à
æ
à

à
à æò
àà æ
à
àà ææòàà ææææòò
à
æ
à
æ òò
à
à æ
ææ
òòòòòòò
æ
òò
ò æ
ò æ
ò ò æ

0.6

0.8

1.0

1.2

1.4

Λ
Figure 3.12: This figure shows a measure of the gaps in ∆(0) vs. λ for different system sizes,
FM . We deem that a gap occurs if |∆(0)(Np + 1) − ∆(0)(Np )| >  (we choose  = 10−15 here). We
can see that with increasing M , the measure of the gaps approaches a step function with the step
position at λ = 1.

3.6

Summary
We now summarize our key results in this chapter:(1) Extended-localized vs KAM-Cantori phases
At the many-body level, the connection between the metal-insulator transition and the
KAM-Cantori transition is signaled in the return map of the local density between nearestneighbor sites. For a given angle φ, the metallic phase with extended single-particle wave
functions exhibits a smooth return map. We identify this with the preserved-invariant

56
KAM tori in phase space at weak perturbations in the corresponding classical system. The
localized phase displays a discrete return map, which we identify with the remaining tori
or Cantori outside the perturbative regime. This behavior is shown in Fig.3.1. In addition,
for the particular filling factors at which the system becomes a band insulator, the returnmap remains smooth for any value of λ (see Fig. 3.2).
Similar behavior can be observed in the corresponding momentum-distribution return map.
However, while the local-density return map loses these characteristic features after averaging over φ, the momentum-distribution return map remains almost unaffected. The
robustness of the momentum return map to phase variations is ideal for the experimental
visualization of the metal-insulator vs KAM-Cantori connection in cold atoms.
(2) Quasi-fractal structures
The introduction of weak disorder modifies the characteristic step-function–Fermi-sea profile of the quasi-momentum distribution. Additional step-function structures centered at
different reciprocal lattice vectors of the QP lattice appear for λ > 0. We refer to thosestructures as the “quasi-Fermi seas”. With increasing λ and ν, the number and widthof the various quasi-Fermi seas increase the fragmentation of the momentum distribution,
turning it into a complex pattern. At λ = 1 the fragmentation becomes maximal, and themomentum distribution evolves into a smooth profile as the system enters the localized
phase (See figures 3.3-3.6).
A more appropriate analogy of such behavior can be observed by considering the set ofvalues taken by quasi-momentum distribution for a given filling factor. In the absence ofdisorder, this distribution can only have the values of 1 or 0, depending upon whether quasi
momentum is greater or smaller than the Fermi quasi momentum. As the strength of thequasi-periodic lattice increases, two distributions of values develop, centered around 0 and
1 respectively. Their width increases with increasing disorder and they overlap exactly at
criticality (Fig. 3.7). We refer to those windows as Arnold tongue-like structures because

57of their resemblance to Arnold tongues in classical mechanical systems.
(3) Bifurcations
The overlap between various Arnold tongue-like structures at λc , can be linked, using thespace-momentum duality transformation, to the appearance of a bifurcation in the density
profile. The bifurcation occurs at a common, but phase-dependent, filling for the various
Fibonacci sites (Fig. 3.8). At generic lattice sites, the filling factor at which the bifurcation
takes place also depends on the lattice site under consideration and can be observed when
the local density is plotted as a function of λ and ν (Fig. 3.9) .
(4) Devil’s staircases
Noise correlations plotted as functions of the filling factor exhibit step-like structures thatevolve into a Devil’s staircase at the onset to the metal-insulator transition ( Figs. 3.10–3.12).
Systems with competing periodicities stand in between periodic and random systems. The
richness and complexity underlying such systems have been studied extensively [46]. Ultracold
atoms are emerging as a promising candidate to simulate a wide variety of physical phenomena. Here
we have shown they offer opportunities to experimentally realize various paradigms of nonlinear
dynamics.
Our focus here was on spin-polarized fermionic systems, since we wanted to look at the simplest consequences of many-body physics in disordered systems. However, Bose-Einstein condensed
systems may also be used as tools for laboratory investigation of various predictions made for thequasi-periodic systems based on single-particle arguments [44]. For example, it might be possible
to confirm the strong-coupling universality prediction, which establishes that the ratio of the single
particle density at two consecutive Fibonacci sites should be a universal number [62].

Chapter 4

Self-trapping dynamics in a 2D optical lattice

Cold gas experiments are well suited to studying out-of-equilibrium dynamics, such as theevolution of many-body correlations, since experimental time scales are relatively slow yet fasterthan typical relaxation and decoherence rates [63, 28, 64]. In this chapter, we introduce an interesting nonequilibrium phenomenon called “self-trapping” that is caused by a competition between
discreteness introduced by the optical lattice and nonlinear effects due to atomic interactions. In
collaboration with experimental colleagues at Pennsylvania State University [65], we found what
seems to be a new type of Macroscopic Quantum Self-Trapping (MQST) effect induced by quantum
correlations in a gas of ultracold 87 Rb atoms trapped in a two-dimensional optical lattice. The lattice potential created coupled arrays of 1D tubes. A striking inhibition of the transverse expansion
dynamics was observed as the depth of the 2D lattice was increased.
Here we describe in detail various theoretical models developed to understand the observed
behavior. The main conclusion is that the ST phenomena observed in this experiment is very likely
caused by the strong correlations between atoms in the same tube. At first, the pure mean-fieldmodel based on the solution of coupled-nonlinear equations fails to reproduce the experimental
observations. It greatly overestimates the initial expansion rates at short times and predicts aslower expansion rate of the cloud at longer times. It also predicts the formation of a hole surrounded by a steep square fort-like barrier that was not observed in the experiment. An improved
theoretical description based on a simplified Truncated Wigner Approximation (TWA) (See Ref.
[66] and references therein), which adds phase and number fluctuations to the initial conditions,

59
pushes the theoretical results closer to the experimental observations, but fails to quantitatively
reproduce them. To further investigate the role of genuine quantum correlations as a localization
mechanism and to benchmark the validity of the approximated TWA, we study the expansion dynamics in a simpler two-tube configuration using exact numerical time-dependent Density Matrix
Renormalization Group (t-DMRG) methods (See Ref. [67, 68] and references therein).
The outline of this chapter is as follows. In section 4.1, we introduce the concept of selftrapping and review the experiments where the MQST has been observed. In section 4.2, we
review the treatment of MQST in terms of a mean-field model and the application of this method
to prior experiments. In Sec. 4.3, we introduce our experimental setup [69] and explain a mean-field
treatment of the dynamics capable of dealing with the axial expansion along the tubes. We compare
our solutions obtained by numerical evolution of the nonlinear equations with the experimental data.
To overcome the limitations of the mean-field model, in section 4.4 we discuss ways to incorporate
quantum fluctuations using the aTWA methods. In section 4.5, we use t-DMRG to study the
problem of atoms tunneling between two tubes and use it to test the parameter regime in which
the GPE and the aTWA are valid. Those calculations show the relevance of quantum correlations
during the expansion dynamics. Finally, we conclude in Sec. 4.6.

4.1

Introduction to self-trapping
The term “self-trapping” was first introduced by Landau to describe the motion of an electron

in a crystal lattice [70]. According to Landau, the electron gets dressed by the lattice polarizationand deformation caused by its own motion. The resulting polaron can be localized or self-trapped in
the presence of strong electron-lattice interactions. Nowadays, the concept of self-trapped polarons
is applied to various condensed matter systems such as semiconductors, organic molecular crystals,
high-temperature cuprate superconductors and colossal magnetoresistance manganates [71, 72, 73].
Still, the theory of the dynamical coupling of a conduction electron to lattice phonons is a complicated highly nonlinear problem that is difficult to tackle even in one-dimensional situations [74].
Many questions remain open.

60
Macroscopic self-trapping in quantum degenerate Bose gases, where mean-field energy gradients suppress tunneling, presents an interesting set of non equilibrium phenomena, with analogs
nonlinear photonics optics [75, 76] and Josephoson junction arrays [77]. Macroscopic self-trapping
has been studied theoretically in all dimensions [78, 79, 80, 81, 82, 83, 84, 85, 86] and experimentallyin double-well systems and in arrays of 2D pancake-like BECs [87, 88]. In those systems, self trapping is caused by a competition between the discreteness introduced by the lattice and nonlinear
effects. However, in these systems, the lattice potential is rigid (imposed by optical fields), and thenonlinearity has its roots in intrinsic interatomic interactions. The self-trapping mechanism is theninherently a many-particle phenomenon, commonly referred to as MQST.
MQST can be understood as wave-packet localization caused by the interaction-induced
suppression of tunneling in the presence of chemical potential gradients. In those cases, a mean
field description in terms of the so called time-dependent Gross-Pitaevskii equation (GPE) is shown
to be sufficient to describe the experiments.
In contrast to other self-trapping experiments, in the experiment we study, the cloud isallowed to freely expand along the axis of the tubes, so that the atom density decreases with time.
The inhibited expansion persists until the density became too low to sustain MQST.
The observed transverse localization was not reproduced by a pure mean-field analysis. Whenphase fluctuations among the tubes were added to the theory in an attempt to approximate theeffect of correlations, the agreement with the experiment improved for shallow lattices. But for
deep lattices, the calculated suppression of tunneling was less complete than what was observedin the experiment, even with maximal phase fluctuations. In Ref. [65] we proposed a possible
mechanism for this effect that cannot be reproduced by modifying mean-field theory. 1D gases
develop quantum correlations to reduce their mean-field energy. Atoms might be prevented fromtunneling because of the mean field energy cost they would have to pay in the adjacent tube, withwhose atoms they are not properly correlated.
A simple understanding of this mechanism can be gained in the Tonks-Girardeau regime
where because of quantum correlations bosons fermionize [89, 5, 90]. Instead of a macroscopically61
occupied single-particle wave function, the many-body wave function is more like a set of spatially
distinct single-particle wave functions to avoid interactions. In this case it is clear that tunnelingbetween adjacent tubes can be suppressed because of the large interaction-energy cost an atom
would have to pay because of mismatched correlations between tubes, i.e., if the only region thatan atom can access via tunneling (which exponentially decreases with distance) in the adjacenttube overlaps with the wave function of an atom already present there.
The succession of MSQT experiments from the double well BEC (coupled 3D gases with nolattice), to coupled 2D gases in 1D lattices, to the coupled 1D gases in 2D lattices, is somethingof a microcosm of the way cold atom physics has developed. As atoms become more confined andthe coupling strength increases, the GPE ceases to capture all the important physics. Quantum
correlations play an enhanced role, giving rise to new phenomena and posing a challenge to theory.
It is natural to speculate that the new type of MSQT that appears in 2D lattice experiments could
also be affecting transport in 3D lattice experiments.

4.2

MQST: an overview
A dilute bosonic gas forms a Bose Einstein condensate (BEC) below a critical temperature.

The GPE has proven to be a very powerful and successful way to describe the dynamics of BECs.
A BEC confined in a double well potential can experience MQST. This was predicted theoretically [78, 79, 80] and has been directly observed experimentally [87]. MQST was induced in theexperiment by creating a population difference between the left and right well. For small population imbalance, the atoms oscillate back and forth between the wells, as expected for noninteracting
atoms, but beyond a critical population imbalance, the mean-field interactions inhibit tunnelingand the atoms become self-trapped.
At the mean-field level, the double-well system can be described by a two-state model as

 


0
−J
d  ψL   UL NL + EL
  ψL 
i~ 
(4.1)

,
=
dt
0
ψR
−J
UR NR + ER
ψR
where ψL and ψR are “macroscopic” wave functions of the particles in the left and right wells, J is62
the tunneling matrix element, and UL,R describe the interaction energy cost of having two particles
in the left or in the right well respectively. NL,R = |ψL,R |2 are the numbers of particles in the two
wells, and EL,R are the respective zero-point energies. These parameters can be determined by theoverlap integrals of the eigenfunctions of isolated wells.
We set EL + ER = 0 and denote the zero point energy difference between the two levels
as ∆E = ER − EL . Assuming that all particles initially occupy one well, when UL,R NL,R → 0,√
Eqn. (4.1) yields Josephson oscillations of the population with a frequency 4J 2 + ∆E 2 and an
2. When ∆E  J, the system is off-resonant and the atoms oscillate rapidly
amplitude 4J 24J+∆E 2
with a small amplitude around the initial configuration.
When ∆E = 0 and |UL NL − UR NR |  J, the atoms also show reduced-amplitude oscillations. In this case, however, the population imbalance is self-locked to the initial value because of
MQST.
MQST in 1D optical lattices has also been studied theoretically and observed experimentally[82, 88, 83, 81]. Atoms were first loaded in a 1D optical lattice with an additional dipole trap,
forming arrays of 2D pancake-like BECs. Then the dipole-trapping beam along the lattice direction
was suddenly removed. Even though the 1D lattice system is more complicated than the double
well, the physics responsible for the MQST is similar in the two cases. In the 1D lattice MQST
manifests as a dynamical localization of an initially prepared wave packet. In contrast to thenoninteracting case, in which a continuous increase of the width of the wave packet with time isexpected, interactions can stop the expansion. The MQST starts at the edges of the cloud where thedensity gradient is the largest. As the system evolves, atoms form a hole at the center surroundedby immobile steep edges. In Ref. [88], the MQST behavior was probed by directly imaging the atom
spatial distribution. At a critical value of interaction energy, a transition from diffusive dynamics(monotonic expansion of the wave packet with time) to MQST was observed.
The local dynamics of MQST in lattices is complex, but the global dynamics, i.e., evolutionof the root-mean-square (RMS) width of the wave packet can be predicted analytically using avery simple model. A variational ansatz has been used to predict when MQST will occur, and the63
result is solely determined by global properties of the gas [81]. It has been applied to 1D latticesystems in which the lattice splits the BEC into an array of two dimensional “disks” or pancakes[88]. The predictions of the model were found to be in qualitative agreement with the experimental
observations. Details of the variational method are presented in the Appendix D.

4.3

MQST in a 2D latticeDuring the expansion, dynamics took place in an array of quasi-1D tubes, where atoms

expanded along the tubes as they underwent their transverse dynamics. This made it possible toobserve a self-trapping transition “on the fly”, as the overall density steadily dropped because of
the axial expansion.
A BEC of N ∼ 3.5 × 105 87 Rb atoms was initially prepared in a crossed-dipole-harmonic
trap with transverse frequency ω⊥ = 2π × 38 Hz and vertical confinement ω|| = 2π × 94 Hz. It
was then loaded into a 2D square-lattice potential with lattice spacing a = 385 nm. Details of theexperimental setup are given in Ref. [69] and the experimental geometry is depicted in Fig. 4.1.
The blue-detuned 2D lattice was ramped up in time according to I(t) ∝ [1 − (t/τ )]−2 , reaching a
final lattice depth of V0 = 7.25 ER , 9.25 ER , 11 ER or 13 ER and creating an array of 1D tubes. HereI is the lattice intensity, τ = 4.15 ms is the time constant, and ER = h2 /(8M a2 ) is the recoil energy,
where d is the lattice spacing. After the initial preparation, all the harmonic confining potentials
were suddenly turned off so that the atoms could expand in the 2D lattice. The expansion dynamics
were investigated by direct imaging. The line of sight was along a direction 45 degrees between thelattice directions (45 degrees from x and y), and the image of the 2D density-distribution squared
was recorded as a function of time. The main finding was a suppressed expansion rate of the RMS
width of the density profile along the lattice, consistent with MQST. However, the expansion rate
did not agree with the mean-field predictions. Moreover, no signature of MQST behavior predictedby the mean-field theory was observed, such as the formation of steep edges or a hole at the centerof the density profile.

64

crosseddipole trap

Bose gases

Figure 4.1: Schematic of the experimental setup. A BEC of 87 Rb atoms was initially prepared in
the crossed-dipole-harmonic trap. A 2D optical lattice along the x-y direction was then adiabaticallyramped on to create an array of quasi-1D tubes. The crossed dipole-trap was then turned off, andthe expansion dynamics in the presence of the 2D periodic potential were investigated by direct
imaging. The line of sight was at 45 degrees from the lattice directions (45 degrees from x and y).
The 2D atom density distribution was recorded after various expansion times.

65
4.3.1

Mean field: Coupled GPEs

For the 2D lattice configuration, the total external potential was given by Vext (x, y) =Vlat (x, y) + Vdip (~x) with Vlat (x, y) = V0 [sin2 (πx/a) + sin2 (πy/a)] the lattice potential and Vdip (~x) =
1
1
2
2
2
2 2
2 M ω⊥ (x + y ) + 2 M ω|| z the confinement introduced by the crossed-dipole-trap. The dipole con-finement was turned off during the expansion, and therefore it was only relevant for determining
the initial conditions. A small antitrapping potential remained because of the Gaussian profile ofthe laser beams that generated the lattice. The parameters of the crossed dipole trap used to set
up the initial state in the theoretical model were chosen so that the atom distribution matched
the initial spatial extent and energy of the atoms in the experiment. The experimental parameters
cannot simply be used because atoms in the experiment do not expand by as much as mean-field
theory predicts when the blue-detuned lattice is turned on. The explanation for this may be related
to the anomalous self-trapping that is the subject of this paper, but a full understanding of that
effect will require future work. The suppressed atom expansion during the turn-on of the latticemight also mean that the initial atom correlations are not the same as in the mean-field model.
Assuming that only the lowest band of the 2D optical lattices was populated, a condition
that was satisfied in the experiment, we write Ψ in terms of the lowest-band Wannier orbitals as,Ψ=X

W (x − an)W (y − am)Φnm (z, t),

(4.2)

m,n

where Φnm is the order parameter describing the tube centered at lattice site {n, m}, and Nnm (t) =RPdz|Φnm (z, t)|2 is its corresponding atom number. The normalization condition n,m Nnm = N
with N the total number of atoms, is satisfied.
By assuming Φnm (z, t) = ψnm (t)φnm [z, Nnm (t)] withRdz |φnm [z, Nnm (t)]|2 = 1 and Nnm =|ψnm (t)|2 and using a similar variational method to the one described in Appendix D, we can obtainthe following effective Hamiltonian#
2 2/3 2/3
9N
= N −4J + Uef f
25
πr2


9
3
= N −4J + U ρa ,
25
"

Hef f

(4.3)

66
where r is the width of the wave packet in the lattice directions (in the lattice units). Here we
1/3

2 a2 U 2
9M ω||g R
4
4have defined the 2D interaction strength as U = a dxdy W (x)W (y). Uef f =
32
h 2 2
i
R
~ ∂
plays a similar role as Uα in the Appendix D, and J = dx W (x − a) 2M−VW (x).
lat
2
∂x


1/3
2/3 2M ω||2 a2
is the 3D density at the center of the cloud. Based on the variationalρ = 4a33 r2N
2π
3U
method, we predict that the MQST in the 2D optical lattice should take place when U ρa3 /J ≥ 100
9 .
Although the variational method roughly describes the global dynamics and gives the threshold value ρ of the MQST, it only explicitly accounts for the expansion dynamics across the lattice
but not along the tube’s direction. We first tried to adapt the variational methods to capture theexpansion in both directions. For that purpose, we assumed the density profile along each tube
φnm (z, t) to have either a Thomas-Fermi (TF) or a Gaussian shape with a width Rnm (t) and aphase δnm (t). In the lattice direction, on the other hand, the evolution of the amplitudes ψnm (t)at each lattice site was numerically computed. Unfortunately, we found that it was quantitatively
accurate only at very short times. The reason is that (in the experiment) atoms in the tubes rapidly
expanded vertically with a number-dependent expansion rate. Since the atoms tunneled between
neighboring tubes at the same time, these two processes destroyed the assumed Gaussian or TF
shape in each tube and the ansatz broke down.
Since the variational method failed, we instead numerically solve for the axial dynamics. We
neglect any temporal dependence of the Wannier functions, which is justified because the interwell
number/phase dynamics are much faster than the time associated with the change in shape of suchfunctions. Generically, the wave field Ψ(~x, t) can still be written in terms of Eqn. (4.2), so wekeep Φnm (z, t) as an arbitrary function. After integrating out the lattice directions, we obtain thefollowing nonlinear Schrödinger equations and solve them numerically.

i~Φ̇nm (z, t) =
~2 ∂ 2
2−+ Vdip (an, am, z) + U d|Φnm (z, t)| Φnm (z, t)
2M ∂z 2−J [Φn+1,m (z, t) + Φn−1,m (z, t) + Φn,m+1 (z, t) + Φn,m−1 (z, t)] .

(4.4)

To characterize the MQST, we numerically compute the root mean square (RMS) transverse
width, xrms , of a vertical slice of the cloud centered at z = 0 and with a width corresponding to67±10% of the total axial width (characterized by the Thomas-Fermi radius, Rz ). The reason for

0.25
0.20
0.15
0.10
0.05
0.00

Acc HΜmms2L

Acc HΜmms2L

selecting a slice instead of integrating over the whole cloud is to consider a sample of approximately
pP
P2
constant density. Theoretically, xrms is directly calculated as: xrms = am ρm m /m ρm ,
R Rz/10P
where ρm = −Rz/10 dz n |Φnm (z, t)|2 .

0

0.25
0.20
0.15
0.10
0.05
0.00

5 10 15 20 25 30t HmsL

0

0.25
0.20
0.15
0.10
0.05
0.00

(b)

Acc HΜmms2L

Acc HΜmms2L(a)

0

5 10 15 20 25 30t HmsL(c)

5 10 15 20 25 30t HmsL

0.25
0.20
0.15
0.10
0.05
0.00
0

10

20 30
t HmsL

40(d)

Figure 4.2: The acceleration of xrms (ẍrms ) obtained from the mean-field calculations for all
lattice depths used in the experiment: (a) 7.25 ER , (b) 9.25 ER , (c) 11 ER and (d) 13 ER . Fromthe sign of ẍrms , the evolution can be separated into three different regimes indicated by a solid
green line and a dashed blue line. The three regimes are discussed in detail in the text. MQST issignaled in the xrms as a negative acceleration. It starts at the solid green line and stops at ∼ tc
(indicated by the dashed blue line) when the acceleration curves cross 0 from below. The values ofU ρ(tc )/Ja3 for the four depths under consideration at tc are {0.4, 3.0, 3.0, 3.0}. Those correspond−1to ρ(tc )/J = {60, 350, 320, 300} µm−3 ER, respectively.

The expansion along the tubes makes MQST a time-dependent phenomenon, since as theatoms freely expand, the density decreases with time until the system becomes too dilute to sustain
self-trapping. To determine the critical value of ρ/J at which MQST stops, we numerically calculateẍrms . We construct an approximation function that interpolates the numerical values of xrms (ti )at recorded time ti . As shown in Fig. 4.2, the curves of ẍrms are not completely smooth, especiallyEint, Ek x, Ekz

0.25
0.20
0.15
0.10
0.05
0.00

68

0

5 10 15 20 25 30t HmsL

Figure 4.3: Time evolutions of: the kinetic energy along the lattice direction (Ekx , dashed orangeline), the kinetic energy along z (Ekz , dotted black line), and the interaction energy (Eint , solid
red line). Here energy is in units of ER . This plot is computed using the mean-field model for a
7.25 ER lattice. The vertical lines are at the same positions as in Fig. 4.2 (a).

at long times, because of accumulated numerical errors. We define the critical time tc to be thetime when the acceleration first crosses zero from below, i.e., it is negative for times shorter thanthe critical time due to the MQST (See Fig. 4.2). Note that tc as well as ρ(tc )/J are only roughly
estimated.
The variational method gives the criterion, U ρ(tc )/Ja3 ≈ 11. Numerically, we found U ρ(tc )/Ja3= {0.4, 3.0, 3.0, 3.0} for lattice depths 7.25, 9.25, 11 and 13 ER , which correspond to ρ(tc )/J =−1, respectively. Although the numerical values of ρ(tc )/J are smaller{60, 360, 320, 300} µm−3 ER

than the value predicted by the variational method, they are almost constant, except for the one
associated to the 7.25 ER curve, consistent with the general concept of a self-trapping threshold.
There are three distinct regimes clearly observed in ẍrms that we use to characterize thedynamics. We will discuss those different regimes in detail below.

4.3.2

General behavior: Three different regimes

Based on the sign of ẍrms , we classify the dynamical evolution of the system into threedifferent regimes: Initial expansion (Exp), MQST, and Ballistic expansion (BE). Each regime has adistinct behavior that can be further characterized by other observables such as the ratio between

69

60

60

40

40

-20 -10 0

20

density(at/μm)	
  density(at/μm)	
  

20
10 20

-20 -1010 20t	
  (ms)	
  t	
  (ms)	
  (b)ρ	
  (at/μm3)	
  

(a)

x	
  (μm)	
  (c)Figure 4.4: Expansion dynamics predicted by mean-field calculations for the lowest lattice-depth(7.25 ER ): Panel (a) displays the evolution of the density profiles at z = 0, after integrating along
one transverse direction. Panel (b) displays the same density profiles, but after integrating along a
direction that is 45 degrees between the x and y axes (as done in the experiment). Panel (c) showstransverse density profiles at z = 0, during the MQST regime at t = 8.7 ms. The insets in Panel(a) and (b) are intersecting profiles at t = 8.7 ms marked by white lines in the 3D density plots
and have the same units as in the 3D plots.

70
interaction and kinetic-energies (Fig. 4.3) and the shape of the density profile (See Fig. 4.4).
• Initial expansion (Exp):The initial expansion regime describes early times during which the atoms at the centerof the cloud expand transversely. We determine this regime by looking at the time period overwhich the transverse RMS width of the cloud expands with a nonlinear rate (positive acceleration,ẍrms > 0). During this time, the density-profile remains smooth. The interaction energy decreases
very rapidly and is converted into kinetic energy along the axial and transverse directions. In
contrast to a noninteracting system, in which the initial potential energy in the trap is converted
to kinetic energy while the width of the quasi-momentum distribution stays constant, the meanfield interactions cause a broadening of the quasi-momentum distribution whose width grows until
it reaches values of quasi-momenta at which the effective mass becomes negative. At this point,
sharp peaks in the density profile develop.
For the 7.25 ER lattice, the initial-expansion takes place within the first 5ms [See Fig. 4.2 (a)]after the parabolic-confinement is turned off. In Fig. 4.4 (a) and 4.4 (b), we also show the time
evolution of a horizontal slice of the cloud z = 0, after having integrated along one transverse
direction (along x), and a direction 45 degrees from the transverse direction (45 degrees from x andy), which is the line of sight direction in the experiment.
• Macroscopic Quantum Self-Trapping (MQST):
At intermediate times (t ∼ 5 − 12 ms for the lowest lattice depth) MQST is signaled by
the RMS width of the cloud. During the MQST, regime the interaction energy decreases slowly
and remains comparable to the kinetic energy of the atoms in the lattice. The interplay between
interatomic interactions and the atomic density gradient at the cloud edges prevents atoms fromtunneling outward. In this regime, the atoms start to pile up at the edges, and a hole forms at thecenter of the cloud [see Fig. 4.4 (c)]. The radial expansion slows-down, or even stops, as illustrated
by a negative xrms acceleration. The slowing down can also be linked to the population of a large
number of quasi-momentum states with negative effective mass.
Global observables, such as the RMS width are easily measured and characterized. They

71
encode some important signatures of MQST. However, MQST is an inherently local effect, andimportant information can be hidden in the global probes. For example, although the formationof steep edges and a hole is shared by the MQST phenomena in both 1D and 2D lattices and isobservable in the RMS width, there are more features in the 2D system.
In the MQST regime, the atom distribution in the lattice is no longer smooth, and its profile
evolves from having a circular shape to a square one. Atoms at 45o to the lattice axes are only
self-trapped at larger radii, since the corresponding density gradients are smaller than in the latticedirections. This gives rise to the square fort-like barrier around the edges, as seen in Fig. 4.5. The
formation of steep edges also happens in the 1D system, but in 2D there is always a direction along
which atoms can tunnel. Since they are not fully frozen, the self-trapped edges are not stationary
as in the 1D case, but instead they evolve with time.
The shape of the integrated-density distribution strongly depends on the direction along
which the integration is performed [See Fig. 4.4(a) and Fig. 4.4 (b)]. For example, while the peaks
at the edges, signaling the hole at the center of the cloud, are clearly visible in Fig. 4.4 (a), when
the imaging direction is along the lattices, they are barely visible when the imaging direction isalong the diagonals as in Fig. 4.4 (b).
• Ballistic Expansion Regime (BE):Because atoms expand axially along the tubes, the interatomic interactions decrease as time
evolves, and at some point, they are no longer strong enough to enforce the self-trapping mechanism.
This determines the interruption of the MQST and the onset of the ballistic regime. In this regime
mean-field interactions have decreased so significantly that an analysis in terms of single-particle
eigenstates and eigenmodes can be carried out.

4.3.3

Comparison with experiment

The experimentally measured RMS widths are slightly different from the xrms defined above.
In the experiment, the density-square distribution, instead of the density, was directly measured.
To mimic the RMS widths measured in the experiment, we define another RMS width, denoted as

y"(μm)"ρ"(at/μm3)"72x"(μm)"(a)(b)

Figure 4.5: Panel (a) is a schematic picture of the tunneling in 2D optical lattices. In 2D, nosite is self-trapped in all directions, so atoms are not fully frozen. The arrows show the directions
that atoms may tunnel, while the crosses mean that tunneling along those directions is forbidden
by MQST. Panel (b) is a density (ρ at z = 0) contour plot of a square fort-like barrier developed
during the MQST regime.

σn2 . To get σn2 , we square the 1D distribution ρm (ρm is obtained by integrating the 3D density
distribution along a line 45 degrees from a lattice direction over a slice centered at z = 0) andconvolve it with a Gaussian of RMS width 1.7 µm (the imaging system resolution) to get the
convolved-density-square distribution, n2m or ρ̃m . The convolution procedure mimics the finitesize resolution of the imaging system. Finally, σn2 is directly calculated by the definition: σn2 =pPP2am ρ̃m m /m ρ̃m . This process removes fine features from the theoretical curves that are notresolvable in the experiment.
Another observable measured in the experiment was the Thomas-Fermi radius Rz of theatomic density along z. It was obtained by fitting the atomic density distribution, after integrating
along the transverse directions, to 4R3 z max[0, 1 − (z/Rz )2 ]. Because the initial width of the cloud issensitive to the loading procedure and because the mean field treatment tends to overestimate theinitial width of the cloud compared to the one measured in the experiment, we use different initial
sizes to compare with the experiment, subject to the constraint of matching axial expansion rates
at long times. The latter just ensures that we are using the correct total energy, which is conserved

73
during the evolution. Comparisons between the numerical simulations and the experimental data
are shown in Figs. 4.6 and 4.7. In general, the theory overestimates the expansion rate of the RMS
widths. Only for 7.25 ER does the mean-field theory roughly capture the experimental observationsat short times. For deeper lattice depths, such as 13 ER , the mean-field theory fails to reproduce
the observed dynamics.
The long-time dynamics of σn2 predicted by the mean-field method is sensitive to the initialwidth of the cloud. When the initial RMS width is larger, the MQST boundary is further awayfrom the center, and thus it takes more time for the RMS width to reflect the presence of MQSTat the edges. In contrast, the MQST is more visible for smaller initial RMS widths, as shown in
Figs. 4.6 and 4.7.
The plots in Fig. 4.4 share the same initial condition as the dashed blue line in Fig. 4.6. With
this initial condition, the MQST effect is weak and the self-trapping signatures, such as peaks at
the edges [Fig. 4.4 (a)] and a hole in the middle [Fig. 4.4 (c)], are not very pronounced. These
self-trapping signatures were not seen in the experiment even when imaged from above. These
characteristic signatures might be visible in future experiments with smaller lattice depths, smaller
initial RMS widths, and larger values of Rz .

4.4

Beyond mean-field model
As shown in Sec. 4.3, the mean-field treatment does not correctly describe the experiment.

A fundamental limitation of this treatment is that it assumes all the atoms are in the condensate
and neglects condensate depletion due to quantum correlations or thermal fluctuations (comingfrom an initial thermal component or nonadiabatic effects as the lattice is turned on). Condensate
depletion modifies the dynamics predicted by the single mode approximation.
One way to quantify the role of quantum fluctuations is to view the system as an arrayof approximately one-dimensional Bose gases. 1D Bose gases are typically characterized by thedimensionless parameter γ [91], the ratio of the mean-field interaction energy per particle calculated
without correlations to the kinetic energy calculated with maximal correlations. In a cigar-shaped

74
200

14

æ
æ

æ

10
8
6

Rz HΜmLΣn2 HΜmL

12
æ

æ

æ

æ

æ æ

æ

æ

æ
æ
æ

100

æ
æ
æ

æ æ
ææ
ææ
ææ
ææ

0

150

50

æ

æ
æ
æ
æ
ææ
æ
æææ

0 5 10 15 20 25 30 35t HmsL

5 10 15 20 25 30 35t HmsL(a)(b)

Figure 4.6: The evolution of (a) σn2 and (b) Rz for 7.25ER from mean field calculations andthe experiment. Each black point is the average of 10 experimental measurements. In Panel(a), the error bars represent random uncertainty, but an overall systematic uncertainty of 0.5µm
associated with the imaging resolution is not included. To illustrate the dependence of the meanfield dynamics on the initial width, we show the evolution of σn2 for different initial conditions forthe 7.25ER lattice. Each line corresponds to a particular initial condition, subject to the constraintof matching axial expansion rates at long times. In Panel (b), error bars are not shown, since
we are mainly interested in the long-time dynamics when Rz is very large and the error bars are
negligible. Even with different initial conditions, the expansion rates along z are almost the same,
therefore in Panel (b) all lines are almost on top of each other.

7.5
æ

6.5
æ

6.0
5.5

æ
ææææ

5.0
4.5

0

æ

æ

æ
æ

ææ ææ
æ
æ ææ æ
æ
æ
æ

æ

10

20 30
t HmsL(a)40

Rz HΜmLΣn2 HΜmL

7.0

æ æ
350
æ
300
æ
æ
250
æ
ææ
200
ææ
ææ
æ
150
æ
æ
æ
æ
100
æ
æ
æ
50
æ
æ
ææ
0
0
10 20 30 40t HmsL(b)

Figure 4.7: Time evolution of (a) σn2 and (b) Rz for 13 ER from mean-field calculations and theexperiment. As in Fig. 4.6, the solid and dashed lines are simulation results, and each corresponds
to a particular initial condition. The long-time dynamics of σn2 are sensitive to the initial widthof the cloud also at 13 ER . The black points are experimental data.

75

Γ
2.
1.
0.5
0.3
0.2
-10 -5 0

510

x HΜmL

Figure 4.8: Dimensionless coupling strength at t = 0 as a function of tube location along the yaxis, if the tubes are approximated to be independent. Here V0 = 13 ER . The approximation isnever fully warranted in the regime of this work, where the tubes are manifestly not independent.
Still, we suspect that this γ gives a sense of the initial axial correlations in the relatively deep
lattices.

trap (ω⊥  ωz ), γ = n1D 2|a1D | is inversely proportional to the one-dimensional gas parametera2

n1D |a1D |, where n1D is the one-dimensional number density of particles, and a1D ≈ − 2a⊥s is theq
one-dimensional scattering length. Here a⊥ = µω~⊥ is the size of the ground state in the transverse
direction, and µ = M
2 is the reduced mass of two colliding particles. A weakly interacting 1D gas(γ  1) is well described by mean-field theory. In contrast, in the strongly interacting regime, γ 
1, atoms avoid each other to reduce mean-field energy, and the bosons behave like noninteracting
fermions [89].
In this approximate treatment of each tube as an independent 1D system, for 7.25 ER , 9.25
ER , 11 ER and 13 ER , the corresponding γ values at time t = 0 in the central tubes (around z = 0)
are 0.088, 0.115, 0.132 and 0.144 and increase monotonically from the centers to the edges. To be
more specific, in Fig. 4.8, we plot γ as a function of tube position for this idealized scenario ofindependent tubes and find that γ can become larger than one at the edges.
The values of γ are not small even at time t = 0, particularly at the edges where MQST is
expected to take place. Moreover, they are expected to increase at short times because of the rapid
decrease in density along each tube caused by the axial expansion. This observation is consistent76
with the failure of the mean-field GP equation to describe the dynamics and the need to accountfor beyond mean-field corrections. The latter is especially true for deeper lattices, for which thepicture of an array of 1D tubes is more appropriate. For lattices shallower than the ones used in
the current experiment, tunneling between tubes is expected to become significant. In this casethe system must be treated as a three-dimensional object. In 3D, quantum correlations are lessrelevant and a pure mean-field analysis should correctly describe the dynamics.
At longer times, the system has expanded so much that the probability of having two particles
close to each other becomes negligibly small, and a pure single-particle analysis of the dynamics
becomes accurate. Note, however, that even in this regime a quantitative analysis is still not trivial,
since the state of the gas in the ballistic regime depends on the conditions at the time at whichcorrelations become irrelevant.
In an attempt to include quantum fluctuations, we use an approximate TWA [66, 92, 93]. The
TWA incorporates leading-order corrections to the dynamics expanded around the classical (GPE)limit. The whole idea of the TWA is that the expectation value of a quantity at time t can be
determined by solving the classical equations of motion from 0 to t and sampling over the Wigner
distribution at time t = 0. The TWA is guaranteed to be accurate at short times. For longer times,
there are also higher-order corrections to the classical trajectories. In a 3D lattice (0D condensate
in each well) there are two limiting regimes where the Wigner function can be easily found and has
been shown to capture the exact dynamics well [66]. One is the very weakly interacting regime inwhich the system is almost an eigenstate of the noninteracting Hamiltonian, a product of coherent
states. The other case is the strongly interacting regime with commensurate filling, in which thesystem is a Mott insulator. In this case the many-body wave function is mostly a product FockQ
QL √state L
j=1 |nij with n the filling factor. The semiclassical wave function is
j=1 n exp (iφj ) andthe corresponding Wigner function is characterized by independent random variables φj uniformly
distributed between 0 and 2π.
The case that we are dealing with now has the complication that not only is it in the intermediate regime where neither of the two limiting cases holds, but also instead of a 0D condensate

77
within a lattice site, we have a quasi-1D gas. Finding the Wigner function therefore is not a trivial
task. Prior work done for two-coupled 1D tubes [94] has already shown the importance of taking
into account quantum correlations for the proper characterization of the dynamics. Consequently,
we expect that beyond-mean-field corrections will play an important role in our system.
With this goal in mind and subject to the limitation that the pure mean-field dynamics
involve the propagation of thousands of coupled GPEs, we implement the TWA in an ad-hoc waythat we refer to as aTWA. We account for initial phase fluctuations within each tube by adding
random phase factors: φnm (t = 0) → ηθnm , where θnm is a uniformly distributed random variable
between 0 and 2π, and η parameterizes the strength of the phase fluctuations. When η = 0, theinitial conditions correspond to a fully coherent array of 1D gases, and when η = 1, we have an
initially fully incoherent array. 0 < η < 1 corresponds to a partially coherent array. Since η is
not determined in the theory, we compare the results for different η to experiment, and choose themost similar one. We average the results over at least 20 initial configurations.
Besides phase fluctuations, we also study the effect of including number fluctuations in theinitial conditions. The results obtained using the aTWA are described in detail in the following
subsections.
The aTWA considerably improves the agreement between theory and experiment. The valueof η that optimizes the agreement with the experiment increases with increasing lattice depth.
This is expected, since as tunneling between tubes becomes weaker, phase coherence is suppressed.
The noise introduced by those random phases quickly generates local self-trapping everywhere,
suppressing the transverse dynamics until, because of the axial expansion, interactions drop so
significantly that the system goes to the ballistic expansion regime. Noise also tends to suppress
the hole formation. See Fig. 4.9.

4.4.1

Comparison between theory and experiment: phase fluctuations

Phase fluctuation among tubes substantially suppresses the expansion dynamics. To illustrate
this effect, consider the toy model of a double well with initial relative phase difference φ(0). Theρ	
  (at/μm3)	
  

78

x	
  (μm)	
  Figure 4.9: Transverse density profiles at z = 0 during the MQST regime at t = 8.7 ms for the 13
ER lattice obtained from the aTWA results.

79
critical value of Λ (Λ is proportional to the ratio of the on-site interaction energy and the tunnelingmatrix element between two wells, see the Appendix D for details) that determines the MQST-todiffusive transition depends on φ as [80],Λc = 2!
p
1 − z(0)2 cos[φ(0)] + 1
,
z(0)2

(4.5)

NL −NR. In this case,where z is the fractional population difference between the two wells, and z = N
L +NR

one can see that even when the initial population imbalance is small, z(0)  1, the system canbecome self-trapped if φ(0) → π (Λc → 1).
Our simulation agrees with this expectation. Fig. 4.10–4.11 show the transverse RMS σn2 andthe vertical Thomas-Fermi radius Rz . The aTWA treatment agrees much better with experimental
observations, especially for the deepest lattices. In Fig. 4.10, we can see that the greater the η, theslower is the expansion rate of σn2 . When η = 1, i.e., when the phase fluctuations are maximal,
the expansion rate of σn2 reaches its smallest value.
Although the modeled curves with the appropriate η are never too far from the data, forV0 ≥ 11 ER , the experiment shows a more delayed and more sudden onset of ballistic expansion
than the model. For V0 = 7.25 ER , η ≈ 0.4 best fits the early evolution, when there is MQST.
For larger V0 , η = 1 clearly fits best for the short evolution time. For V0 = 7.25 ER , η ≈ 0.4,
the three regimes predicted by the pure mean-field model are still visible, although ẍrms becomes
barely negative in the MQST regime. On the other hand, for the deeper lattices, the random
initial phases cause a fast development of site-to-site density fluctuations and the formation of
localized domains randomly distributed throughout the tube array. In contrast to the pure meanfield model, localization occurs without the formation of sharp edges, there is no reflection fromthem, and ẍrms never becomes negative. Instead, ẍrms asymptotically approaches zero from above
as shown in Fig. 4.12. MQST always disappears when the density drops below a critical value that
depends on V0 . The model predicts the following values for the transition to a ballistic regime,−1
ρ(tc )/J = {840, 1900, 1970} µm−3 ER
for V0 /ER = {9.25, 11, 13}. We can compare these values

14

9

12

8

æ

10

æ

8

æ

æ

æ æ

æ

æ

æ

Σn2 HΜmL

Σn2 HΜmL

80

æ æ
ææ
ææ
ææ
ææ

6

0

æ
æ
æ æ

7

æ
æ

æ
æ æ
ææ

6
æ æ

5æææ
0

5 10 15 20 25 30 35t HmsL

5 10 15 20 25 30t HmsL

(a)

(b)

7.5

8

æ
æ

æ

æ
æ æ
æææ
æ

6
æ

5

æ æ æ ææ æ

Σn2 HΜmL

Σn2 HΜmL

7.0

7

æ

6.5
æ
æ
æ
æ ææ
æ
æ æ æ
æ ææ æ

6.0
5.5

æ
ææææ

0

5 10 15 20 25 30t HmsL
(c)

5.0

0

æ

æ
æ

æ

10

20 30
t HmsL

40(d)

Figure 4.10: Mean-field, aTWA, and experimental results for the evolution of σn2 for all lattice
depths used in the experiment. The lattice depths are (a) 7.25 ER , (b) 9.25 ER , (c) 11 ER and (d)13 ER . At low lattice depth, a small η seems to account for the observed behavior at short times.
The dashed blue line is the result of a mean-field calculation of the dynamics with no random phasebetween the tubes (η = 0). The dotted purple line is for η = 0.2, the dashed-dotted green lines
are for (a) η = 0.4 and (b) η = 0.5, and the solid orange line is for η = 1. The dashed-vertical red
lines indicate the transition from MQST to BE predicted by the aTWA. The yellow shadow regions
indicate the tc inferred from the experimental data.

81

Figure 4.11:
Evolution of Rz for all lattice depths used in the experiment. The dots show
experimental data, and the solid and dashed lines show the predictions obtained from the aTWA
for different η. The lattice depths are (a) 7.25 ER , (b) 9.25 ER , (c) 11 ER and (d) 13 ER . Rz isobtained by fitting the density along z to a Thomas-Fermi profile. The theoretical curves use thesame parameters as the ones shown in Fig. 4.10.

250

200
æ
æ

150

Rz HΜmLRz HΜmL

200
æ
æ
æ

100
50

æ
æ
æ
ææ
ææ
æ
æææ

æ

æ

æ

0

ææ
ææ
æ
ææ
æ

æ

æ

200
æ
æ
ææ
æ
æ
ææ

100
æ
æ
æ
æ
æ
ææ

0 5 10 15 20 25 30t HmsL(c)

Rz HΜmL

æ

150

æ æ

0 5 10 15 20 25 30t HmsL(b)

250

Rz HΜmL

100

æ

(a)

0

æ

æ

æ

50

0 5 10 15 20 25 30 35t HmsL

50

æ

150

æ æ
350
æ
300
æ
æ
250
æ
ææ
200
ææ
ææ
æ
150
æ
æ
æ
æ
100
æ
æ
æ
50 æ æ
æ0æ
0 10 20 30 40t HmsL

(d)

82

0.15

Acc HΜmms2L

Acc HΜmms2L

Figure 4.12: ẍrms calculated using the aTWA for different lattice depths and for the corresponding
optimal η, (a) 7.25 ER , η = 0.4 and (b) 13 ER , η = 1. The boundaries between MQST and BE
are indicated by the dashed blue lines. In panel (a), the boundary between EXP and MQST isindicated by the solid-vertical green line. In panel (b), because ẍrms is almost 0 at all times, it canbe interpreted that the self-trapping starts at t = 0 ms.

0.1
0.05
0
0

0.15
0.1
0.05

5 10 15 20 25 30t HmsL

0
0

10

(a)

20 30
t HmsL40(b)

r2Hx,z=0LH1êmmLFigure 4.13: Transverse density profiles: The solid and dashed lines are obtained from the aTWA
results for the V = 7.25 ER lattice and for η = 0.4. The profiles are taken at z = 0 and viewed
transversely at an angle 45o from the lattice axis and normalized to 1. The solid blue lines indicate
the theoretical profiles before convolution, and the dashed black lines after convolution. The redpoints are averaged experimental data of 10 measurements.
0.08
t=0.1ms

0.06

t=4ms

0.04
0.02
0.00
-40

-20

0

20

40 -40

-20

0

20

40

r2Hx,z=0LH1êmmLX (µ m)0.08
t=20mst=10ms

0.06
0.04
0.02
0.00
-40

-20

0

20

40 -40

-20

0

20

40

r2Hx,z=0LH1êmmLX (µ m)0.08
t=26ms0.06

t=32ms

0.04
0.02
0.00
-40

-20

0

20

40 -40X (µ m)

-20

0

20

40

83

r2Hx,z=0LH1êmmLFigure 4.14: Transverse density profiles: The solid and dashed lines are obtained from the aTWA
results for the V = 13 ER lattice and for η = 1. The profiles are taken at z = 0 and viewed
transversely at an angle 45o from the lattice axis and normalized to 1. The solid blue lines indicate
the theoretical profiles before convolution and the dashed black lines after convolution. The redpoints are averaged experimental data of 10 measurements.
0.08
t=4ms

t=0.1ms

0.06
0.04
0.02
0.00

-20 -100

10-20 -10

20

0

10

20

r2Hx,z=0LH1êmmLX (µ m)0.08
t=20mst=10ms

0.06
0.04
0.02
0.00

-20 -100

1020-20 -100

1020

r2Hx,z=0LH1êmmLX (µ m)0.08
t=26ms0.06

t=32ms

0.04
0.02
0.00

-20 -100

1020-20 -10X (µ m)

0

10

20

84
with the ones measured in the experiment that occur when the ratio ρ(tc )/J = {1190 ± 310, 860 ±−1190, 910 ± 210} µm−3 ER
for V0 /ER = {9.25, 11, 13}, respectively. We can see that the ratio of

ρ(tc )/J in the experiment is roughly constant. However, while the values of ρ(tc )/J predicted by
the aTWA are almost the same for 11 ER and 13 ER , ρ(tc )/J for 9.25 ER is only half of that value.
On the other hand, the actual ρ(tc )/J value predicted by the aTWA at 9.25 ER is closer to theobserved value.
Figures 4.13–4.14 show the density profiles measured experimentally for the shallowest anddeepest lattices and the corresponding theoretically computed profiles. For the 7.25 ER lattice, one
can see a good agreement after the theory curves are convolved to account for the imaging system
resolution. However, the theoretical distributions are not well fitted by a Gaussian at longer times.
The experimental distributions, on the other hand, always retain a Gaussian profile. In general,phase fluctuations prevent the hole formation and help to keep the density distributions closer toa Gaussian profile.

4.4.2

Number fluctuations

So far we have included only phase fluctuations. However, we know that in the superfluid
regime, relative atom number fluctuations provide the dominant beyond-mean-field corrections.
To investigate the role of number fluctuations in the expansion dynamics, we introduce number√
fluctuations in the initial conditions. We implement this by allowing ≈ 1/ Nnm fluctuations in
the atom number in each tube, with Nnm (0) being the total number of atoms in a tube centered at{n, m}. We focus on the case V0 = 7.25 ER . The results are shown in Fig. 4.15, where we can seethat number fluctuations do suppress the expansion rate of the cloud, but by only a small amount.

4.4.3

Phase and number fluctuations

In the intermediate regime—between the superfluid and Mott phases—both quantum andnumber fluctuations are relevant, but constrained by the corresponding Heisenberg uncertainty

85

Σn2 HΜmL

14
12
10
8
ææ

6 ææææææ
0 5

æ æ

æ

æ æ
æ æ

æ

æ

æ

æ

10 15 20 25 30t HmsL

Figure 4.15: Time evolution of σn2 computed for a 7.25 ER lattice, and for η = 0. The dashed
red line shows the result when number fluctuations
are neglected, and the solid blue line shows the√result when number fluctuations (≈ 1/ Nmn ) are included. The plot clearly shows that number
fluctuations suppress the expansion, but only by a small amount.

Σn2 HΜmL

12
10

æ

8
ææ

6 ææææææ
0 5

æ æ

æ

æ

æ
æ æ

æ

æ

æ

10 15 20 25 30t HmsL

Figure 4.16: aTWA results for the time evolution of σn2 for a 7.25 ER lattice, and for η = 0.35.
The dashed red line shows the result when only phase fluctuations
are included, and the solid√
blue line shows the result when both number fluctuations (≈ 1/ Nmn ) and phase fluctuations areincluded. The plot shows that number fluctuations slightly suppress the initial expansion rate.

relation. To get an upper bound of the amount of localization predicted by the aTWA, we add√
both number fluctuations (≈ 1/ Nnm ) and phase fluctuations in the initial conditions. The results
are shown in Fig. 4.16. Indeed the addition of both phase and number fluctuations helps to suppress
the expansion rate of the cloud for moderate η (as in the 7.25 ER case) but for large η, the additionof number fluctuations is barely noticeable. In other words, the curves with number and phasefluctuations are quite similar to the curves with phase fluctuations only.
In summary, so far we have accounted for the quantum correlations developed in the tubes

86by adding an overall random phase in a phenomenological way. Although this procedure pushes
the theory closer to the experimental observations, a fundamental problem with the aTWA is thatit misses the correlations within each tube. Those might be the key ingredient responsible forthe observed localization at short times because of the additional energy cost they impose whenatoms tunnel from one tube to the next. Those correlations are expected to play a major role atshort times when interactions are dominant. A significantly better treatment could be obtained bycomputing the full Wigner function of the interacting systems. The latter could, in principle, be
achieved by starting with the Wigner function of the noninteracting system (which is known) andthen by adiabatically increasing the interactions to the desired value. The large number of degrees
of freedom in consideration substantially complicates this implementation.

4.5

Beyond the aTWA—two coupled tubes
To test the validity of the various approximations used so far, we now turn to a simplified

model system that is amenable to a numerically exact treatment using the t-DMRG. In the following, we will consider a configuration of two-coupled tubes (see Fig. 4.17) and compare thedynamics of this minimal system to the one obtained with the pure GPE and the aTWA methods.
To be more specific, we replace in our toy two-tube model the axial parabolic-trapping potential
along the tubes by a box potential that is easier to treat numerically, and we create the population
imbalance at the beginning of the dynamics by choosing different chemical potentials on the two
tubes. We then set the chemical potentials to zero instantaneously to mimic the time evolution
in the experiment. Note that when removing the box potential in this setup, the dynamics atthe edges of the cloud can dominate the tunneling between the two tubes. Therefore, we find it
more convenient to keep the box turned on throughout the time evolution. Moreover, to avoid thecomplications of dealing with a continuum system, we assume that along each tube there is a weaklattice (with L lattice sites) and use the tight-binding approximation to describe the dynamics.

87

Figure 4.17: Schematic configuration of the two-tube system.

This gives rise to a Bose-Hubbard model,
H = −Jk

 XU
 X
X X  †X †âl,i âl,j + h.c +n̂l,i (n̂l,i − 1) − J⊥â1,i â2,i + h.c −µl n̂l,i , (4.6)
2

l=1,2 <i,j>

l,i

i

l,i

where â†l,i and âl,i are bosonic creation and annihilation operators, and n̂l,i = â†l,i âl,i gives the number of particles at site i on tube l. Jk is the tunneling matrix element between nearest neighboringsites within each tube, J⊥ is the tunneling matrix element between tubes, U is the on-site interaction strength and µl is the chemical potential on tube l (l = 1, 2). The difference in chemical
potentials µ1 − µ2 equals the bias between the two tubes that is adjusted to obtain the desired
population on each of the tubes. Note that we assume tunneling between the tubes takes place atnearest neighbouring sites only, so that we obtain a ladder geometry.

4.5.1

The DMRG method

DMRG is a numerical method that is capable of obtaining ground-state properties of (quasi)onedimensional systems with very high efficiency and accuracy for lattices with up to several thousandsites. This is achieved by working in a truncated basis of eigenstates of reduced density matrices
obtained for different bipartitions of the lattice. The so-called discarded weight, which is the sumof the weights of the density-matrix eigenstates that are neglected and which should be as small
as possible [67] quantifies the error of the method. Also, its time-dependent extensions (t-DMRG)
can treat the real-time evolution of strongly correlated quantum many-body systems substantiallylarger than the ones amenable to exactly diagonalizing the Hamiltonian and with an accuracy that

88case 1case 2case 3case 4case 5

U Ntotal
30
30
30
60
Ntotal
30

Jk
1
1
1
1
Jk
1

J⊥
0.025
0.05
0.1
0.05
J⊥
0.05

L ( lattice sites)
40
40
40
40
L ( lattice sites)
40

Table 4.1: Set of parameters and different initial conditions used to study the dynamics. Ntotal is thetotal particle number in the ladder. U is the on-site interaction strength, Jk is the nearest-neighbortunneling matrix element within each tube, and J⊥ is the tunneling matrix element between tubes.

can be, at short and intermediate times, similar to what is achieved in ground-state calculations.
Here, we use this accuracy to obtain high-precision numerical results. During the evolution, we aim
for a discarded weight of < 10−9 and keep up to 500 density-matrix eigenstates for systems with
L = 40 sites in each tube. We apply a time step of ∆t = 0.01/Jk and find a discarded weight of. 10−6 at the end of the time evolution displayed in the plots.
In the 1D Bose-Hubbard model, γ ≈ UJ . Therefore in our ladder configuration, we willuse γ = JU to quantify the interaction strength of the system. This is expected to be a goodkcharacterization when the intertube tunneling is much smaller than the intratube tunneling (J⊥ Jk ). We set Jk = 1 in all cases and vary both γ (by changing U ) and the density. Table 4.1
summarizes the set of parameters that we use to study the dynamics. Ntotal is the total particle
number in the system.

4.5.2

Short-time dynamics

Even in this simpler two-tube system, the many-body dynamics after a quench can be quite
complicated. To compare the differences between results obtained using the mean-field theory,
the aTWA, and the t-DMRG, as well as to check the validity of the aTWA, we first focus on the
short-time dynamics, for which we can obtain analytical expressions and numerical results.
A quantity that gives relevant information about the system’s dynamics is the total particleP
number in each tube, i.e., Nl = i nl,i with nl,i = hn̂l,i i. Heisenberg’s equations of motion for the89
atom number in the tubes then read (~ ≡ 1),X †dNl = −iJ⊥hâl0 ,i âl,i i − hâ†l,i âl0 ,i i ,dtioJ⊥ U X n
d2
†
†
2
0
0
N
=
−
h
n̂
−
n̂
,
â
â
+
â
â
l
l,i
l ,i l0 ,i l,il,i l ,i i − 2J⊥ (Nl − Nl0 ),
dt2
2

(4.7)
(4.8)

i

where l0 6= l denotes the position in either tube 1 or 2.
Generally, for the ground state, Im[ha†l,i al0 ,j i] = 0 (there are no currents in the system), andd
Nl t=0 = 0.
therefore the first derivative of Nl is zero at t = 0, i.e., dt

The second derivative of Nl is generally nonzero when there is a population imbalance between
the two tubes, because of the second term on the right hand side of Eqn. (4.8). The first term
accounts for the role of interactions in the dynamics and involves four-particle correlation functions,which can play a fundamental role in the dynamics. For strongly repulsive bosons (hard-corebosons), on-site occupancies > 1 are energetically forbidden, so the four-particle correlations vanish.
Therefore only the term proportional to the population imbalance contributes to the dynamics.
Hence, the dynamics of hard-core bosons are governed by,d22(Nl − Nl0 ).
Nl = −2J⊥
dt2

(4.9)

Analytical expressions for the short-time dynamics can also be obtained using the mean-field (GPE)approximations and the aTWA.
At the mean-field (GPE) level, the operators âl,i are replaced with c-numbers ψl,i . At t = 0,
the ground state corresponds to ψl,i =√

nl,i eiφ , and there is a global phase for all sites. Without

loss of generality, ψl,i (0) can be chosen as√

nl,i . Under the aTWA, ψl,i is chosen asa global phase factor for each tube, i.e., ψl,i (0) →√√nl,i times

nl,i eiθl ·η . As we explained in section 4.4, η

is a tunable parameter thatcharacterizes the strength of phase fluctuations, and θl is a uniformly
distributed random variable between 0 and 2π. The average can be computed analytically, and in
the aTWA we obtain that the dynamics is governed byd2
sin2 (πη)
Nl − Nl0 p
2
0
0
N
N
N
.
=
−2J
(N
−
N
)
−
2J
U
l
l
⊥
l l
l
⊥
dt2
L
π2η2
t=0

(4.10)

90
U=1
U=2
U=3
U=5

η for case 1
0.31
0.46
0.56
0.69

η for case 2
0.30
0.45
0.56
0.70

η for case 3
0.30
0.46
0.57
0.71

η for case 4
0.27
0.39
0.48
0.61

η for case 5
0.30
0.39
0.45
0.54

Table 4.2: Fitted values of η for different set of parameters used to investigate the dynamics: at
short times the t-DMRG dynamics are well captured by the aTWA when η is used as a fitting
parameter.

The second term in the above equation equals −2J⊥ U(Nl −Nl0 ) √Nl Nl0 when η = 0 and goes to 0Ldwhen η = 1, agreeing with the GPE and hard-core boson limits. Since dt
Nl t=0 = 0, to lowest

 2
d
2
order, Nl (t) = 21 dt
2 Nl t=0 t + Nl (0) at short times. Numerically, we can either directly calculated2Nfrom the initial wave-function at t = 0 or extract it from Nl (t) by fitting it to a quadratic
dt2 l t=02dpolynomial. Here we use the second (fitting) method to get dt
2 Nl t=0 .

At t = 0, Eqn. (4.10) becomes,ẍ|t=0 = f (x) + g(x)sin2 (πη),π2η22 (1 − 2x), and g(x) ≡ 2J U Ntotal (1 − 2x)
1, f (x) ≡ 2J⊥
where x = NN
⊥
L
total

(4.11)px(1 − x). When x = 0, thepopulation imbalance between the two tubes is maximum, and all particles occupy the same tube.
When x = 1/2, atoms are equally distributed between the tubes.
The fitted values of η that best reproduce the t-DMRG calculations obtained for different
parameter regimes are shown in Table 4.2. Figure 4.18 shows comparisons between the aTWA
results and the t-DMRG results for the optimal value of η found for one of the parameters displayedin Table II. In the GPE limit (η = 0), only the product U Ntotal enters in the equations of motion,
meaning that, as long as U Ntotal is kept constant, the dynamics of the system are independent of
U . This feature can be clearly seen in Eqn. (4.10). This is certainly not the behavior predictedby the t-DMRG solutions that show a slow down of the dynamics with increasing U . The TWA
is expected to fully capture the role of quantum correlations at short times. Interestingly, we find2that ẍ|t=0 computed by t-DMRG is extremely well described by f (x) + g(x) sinπ2(πη), where η isη2
the only fitting parameter. The excellent agreement between t-DMRG and aTWA allows us toÐ
N 1 HtL t=0
Ntotal

91

0.020
0.015
0.010
0.005

à

à

à

à

à

à
à
à
ò ò
à
ò ò
à
ò
à ò
à
ò
àòò
ì
ì
à
ì
ò
àòò ìì ì
ì ì
à
ì
ò
ò
ì
ì
òì
ò
ìì æ æ æ
ì
à
ò
æ æ æ
æ
æ
æ
ì
æ
æ
æ
æ æ
ì ò
à
æ
ì ò
æ
æ
ì
æ
à

0.1

0.2

0.3

0.4

à
ò
ì
æ

N1 H0L
0.5 Ntotal

1
Figure 4.18: NN̈|t=0 (in units of Jk2 ) as a function of NN1 (t=0), J = 0.05, U Ntotal = 30. The
total
total
solid black line shows the GPE limit, and the dash-dotted gray line the hard-core boson limit [seeEqn. (4.11)]. The red squares, magenta upper-triangles, blue diamonds, and orange disks show the1
initial curvature NN̈
|t=0 extracted from the t-DMRG results at U = 1, U = 2, U = 3 and U = 5,
total
respectively. The dashed black lines show the aTWA solutions using η as a fitting parameter. The
values of η are listed in Table 4.2 . This plot is computed using the “case 2” parameters shown in
Table 4.2.

conclude that the aTWA captures well the quantum correlations parameterized by η at short times
in this simple two-tube model. We also see, that when U Ntotal is fixed in the parameter regime
under consideration, the fitted values of η that best reproduce the t-DMRG calculations are almost
independent of the initial population imbalance, i.e., 1 − 2x, and the tunneling between tubes, i.e.,J⊥ , when J⊥  Jk .

4.5.3

Long-time dynamics

In Fig. 4.19 we compare the GPE, the aTWA and the t-DMRG dynamics. We see significant
deviations between the approximate methods and the t-DMRG at longer times. This shows that the
dynamics of the Bose-Hubbard model on this ladder geometry is dominated by correlation effects.
Although at short times the t-DMRG results are well described by the aTWA results with anappropriate value of η, at longer times the aTWA fails to reproduce the many-body dynamics. Thisbehavior confirms the relevance of quantum correlations neglected by the aTWA. Those deviations

92
are consistent with the deviations we saw when trying to model the experimental data with theaTWA.
Figure 4.19 shows that the weakly interacting regime is the regime where the GPE and theaTWA (with small η) predictions are closer to the t-DMRG results for longer times. The aTWA
gives slightly better agreement at short times, indicating that for weak interactions, it better
captures correlation effects. In the strongly interacting regime, tunneling is greatly suppressed,
especially in the presence of a large initial population imbalance. While in the GPE and theaTWA, N1 (t) oscillates with more or less constant frequency and amplitude, in the t-DMRG, N1 (t)oscillates with a lower frequency and gradually “relaxes” towards an equal distribution of atoms in
the two tubes.
To further quantify the quantum correlations (which are not captured in the GPE and are
only partially captured in the aTWA) in the ladder system, we calculate the local pair correlation
function, g (2) , which is proportional to the probability of observing two particles in the same latticesite [39, 1]. For a 1D gas, when γ  1, g (2) → 1, and when γ  1, g (2) → 0. In our quasi 1D(2)

system, we define glin the tube l as,
(2)

gl
(2)

Figure 4.20 shows that gl(2)

γ goes to infinity. gl

=

L X 2
hn̂l,i i − hn̂l,i i.
Nl2 i(4.12)

decreases with increasing U (larger γ) and approaches zero whenis a slowly varying function when compared to the evolution of the total

atomic population in each tube (Nl ).
By using t-DMRG as a means to benchmark the GPE and the aTWA, we clearly see thatwhile the aTWA is accurate at short times, both approaches break down at longer times. We canalso clearly observe the complexity of the nonequilibrium behavior even in this simple coupled-tube
model.
It is important to keep in mind that the two-tube model used in this section is much simpler
than the multitube system of the experiment, where in addition the axial external confinement was
turned off, and no lattice was present along the tubes. Consequently, the conclusions of this analysis

25

25

20

20

UN1

UN1

93

15
10

0

2

4

6
t

8

15
10

10 12

0

2

4

(a)

10 12

3
t

4

5

25
N1

N1

8

(b)

25
20
15

6
t

0

1

2

3
t

4

5

6

20
15

0

1

2

(c)

6

(d)

0.5
0.4
0.3
0.2
0.1
0.0
0

gH2L
2

gH2L
1

Figure 4.19: Evolutions of U N1 (in units of Jk ) and N1 as functions of time (in units of 1/Jk )calculated using the DMRG and GPE methods. In panel (a) and (b), the values of U are 1 and
5, respectively, while U Ntotal is fixed to 30. In panel (c) and (d), the values of U are 2 and 5,
respectively, while Ntotal is fixed to 30. In each panel, the solid line is the t-DMRG solution, andthe dashed black line is the mean-field calculation with no random phase (GPE). The dotted grayline is for (a) η = 0.3, (b) η = 0.7, (c) η = 0.4 and (d) η = 0.53.

2

4
(a)

(2)

6
t

8

10 12

0.6
0.5
0.4
0.3
0.2
0.1
0.0
0

2

4

6
t

8

10 12

(b)

Figure 4.20: gl as a function of time in different regimes. The correlations are extracted from theP(2)t-DMRG results for the parameters describing case 2. Panel (a) shows g1 , NL2 i hn̂21,i i − n1,i , and1P(2)
Panel (b) shows g2 , NL2 i hn̂22,i i − n2,i . In each panel, from top to bottom, the lines correspond2
to U = 1, U = 2, U = 3, and U = 5, respectively.

94
cannot be simply compared to the experimentally observed dynamics. However, this simple twotube model clearly shows that correlations can play a dominant role in suppressing the tunnelingdynamics, without the need for density gradients. 1

4.5.4

Discussion

Although we see decent agreement between the aTWA and the DMRG results at short times inthe two-tube model, in the real experiment, we find that the aTWA deviates from the measurements
almost as soon as the expansion starts. We will now discuss the various possible explanations forthe experiment-theory discrepancy.
First, although the aTWA adds phase fluctuations that arise from quantum correlations, they
are added in an ad hoc way, and not as the result of the system lowering its energy. In a complete
theory (and presumably in the experiment) correlations build up in order to decrease the meaninteraction energy. Since the wave function of a tunneling atom will not in general be appropriately
correlated with the atoms in an adjacent tube, it would have to pay much of the interaction energy
cost that the correlations avoid. If that energy cost is larger than the tunneling energy, then hopping
between tubes will be suppressed. Although γ is not initially large near the central tubes, the fastdensity decrease due to the axial expansion could dynamically increase the quantum correlations,so that this mechanism could play an important role in suppressing transverse expansion.
Second, it is possible that the initial phase variations along the tubes, which are left out of
the aTWA, could be important. Note, however, that phase variations along tubes are generated
dynamically in the aTWA when atoms tunnel between tubes. Still, the resulting phase fluctuations
may not be sufficient to capture the quantum dynamics.
Finally, we have seen recent experimental evidence that the initial atom distribution among
the central tubes is flatter than a Thomas-Fermi profile. Such a distribution is probably a conse1

In the strong coupling limit, correlation-suppressed tunneling is similar to Pauli Blocking. It does not requiredensity gradients, which is why in the experiment the expansion of the central part of the tube bundle is fully
suppressed until the overall density, and hence the mean-field cost of tunneling, gets sufficiently small. We speculatethat such a process might tend to suppress mass transport whenever there is intermediate or strong coupling, regardlessof the dimension [48, 95].

95
quence of the known nonadiabaticity of the 2D lattice turn-on, which leads to much-higher initial
densities than adiabatic turn-on predicts, in a way that is quite insensitive to turn-on time. We
expect that a more flat-top initial distribution would lead to less initial expansion, but it is notclear that it would be decreased enough to explain the experimental results. The limited spatial
resolution of the experiment prevents us from directly checking the distribution of atoms among
tubes and performing a calculation based on the result.
In summary, there remain several logical possibilities for the discrepancy between experimentand theory here. Further investigation is needed.

4.6

Conclusion
We have presented various mean-field and beyond-mean-field models to describe our recent

experimental observations of interaction-induced localization effects during the expansion of anarray of coupled 1D tubes. In contrast to previous 1D-lattice experiments, where a pure mean-fieldmodel was able to capture the observed MQST behavior, in our case, we find important corrections
induced by quantum fluctuations. Thermal fluctuations and nonadibatic loading conditions, which
are not accounted for in our analysis, may give rise to similar effects, but we have not explored
those in this work. 2
The addition of phase fluctuations in the initial conditions (aTWA) coarsely captured themain self-trapping features seen in the experiment, but underestimated the observed localization.
The fact that the aTWA does not properly account for the full quantum correlations could be
responsible for the failure of the theory to reproduce the experiment. Comparisons with t-DMRG
calculations performed in a simple two-tube model find very good agreement between the aTWA
and the t-DMRG at short times, but a break down of the aTWA at longer times. The two-tube
calculations support the fact that quantum correlations, not captured by the mean-field model, canplay an important role in suppressing tunneling. Further comparisons among the various methods
2

It is hard to adiabatically transform a 3D gas into a 1D gas in theory [96], and there is some irreversibility inturning on the lattice.

96
used for modeling the dynamics and experiments are needed to help shed light on the complex
many-body dynamics of coupled-correlated systems.

Bibliography[1] I. Bloch, J. Dalibard, and W. Zwerger. Many-body physics with ultracold gases. Rev. Mod.
Phys., 80:885, 2008.
[2] S. Peil et al. Patterned loading of a bose-einstein condensate into an optical lattice. Phys.
Rev. A, 67:051603, 2003.
[3] P. Blair Blakie and Charles W. Clark. Wannier states and bose-hubbard parameters for 2d
optical lattices. J. Phys. B: At. Mol., Opt. Phys., 37:1391, 2004.
[4] C. Orzel et al. Squeezed states in a bose-einstein condensate. Science, 291:2386, 2001.
[5] T. Kinoshita, T. Wenger, and D. S. Weiss. Observation of a one-dimensional tonks-girardeaugas. Science, 305:1125, 2004.
[6] B. Paredes et al. Tonks-girardeau gas of ultracold atoms in an optical lattice. Nature, 429:277,
2004.
[7] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch, and I. Bloch. Quantum phase transition
from a superfluid to a mott insulator in a gas of ultracold atoms. Nature, 415:39, 2002.
[8] M. Lewenstein et al. Ultracold atomic gases in optical lattices: mimicking condensed matterphysics and beyond. Advances in Physics, 56:243, 2007.
[9] J. D. Jackson. Classical Electrodynamics. John Wiley & Sons, Inc., New York, 1999.
[10] R. Grimm, M. Weidemüller, and Y. B. Ovchinnikov. Optical dipole traps for neutral atoms.
arXiv: physics/9902072v1, 1999.
[11] A. M. Rey. Ultracold bosonic atoms in optical lattices. PhD thesis, University of Maryland,
2004.
[12] N. W. Ashcroft and N. D. Mermin. Solid State Physics. Brooks/Cole, Cengage Learning,
1976.
[13] J. M. Ziman. Principles of the Theory of Solids. Cambridge University Press, 1964.
[14] M. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher. Boson localization and thesuperfluid-insulator transition. Phys. Rev. B, 40:546, 1989.
[15] J. Hubbard. Electron correlations in narrow energy bands. Proceedings of the Royal Societyof London Series a-Mathematical and Physical Sciences, 276:238, 1963.

98[16] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller. Cold bosonic atoms in opticallattices. Phys. Rev. Lett., 81:3108, 1998.
[17] T. Stöferle et al. Transition from a strongly interacting 1d superfluid to a mott insulator.
Phys. Rev. Lett., 92:130403, 2004.
[18] I. B. Spielman, W. D. Phillips, and J. V. Porto. Mott-insulator transition in a two-dimensional
atomic bose gas. Phys. Rev. Lett., 98:080404, 2007.
[19] I. B. Spielman, W. D. Phillips, and J. V. Porto. Condensate fraction in a 2d bose gas measured
across the mott-insulator transition. Phys. Rev. Lett., 100:120402, 2008.
[20] T. Esslinger. Fermi-hubbard physics with atoms in an optical lattice. arXiv:1007.0012, 2010.
[21] E. Altman, E. Demler, and M. D. Lukin. Probing many-body states of ultracold atoms via
noise correlations. Phys. Rev. A, 70:013603, 2004.
[22] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell. Observationof bose-einstein condensation in a dilute atomic vapor. science, 269:198, 1995.
[23] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and
W. Ketterle. Bose-einstein condensation in a gas of sodium atoms. Phys. Rev. Lett., 75:3969,
1995.
[24] C. J. Pethick and H. Smith. Bose-Einstein Condensation in Dilute Gases. Cambridge University
Press, Cambridge, 2002.
[25] M. R. Andrews, C. G. Townsend, H.-J. Miesner, D. S. Durfee, D. M. Kurn, and W. Ketterle.
Observation of interference between two bose condensates. Science, 275:637, 1997.
[26] R. Hanbury Brown and R. Q. Twiss. Correlation between photons in two coherent beams oflight. Nature, 177:27, 1956.
[27] N. Gemelke, X. Zhang, C-L. Hung, and C. Chin. In situ observation of incompressible mottinsulating domains in ultracold atomic gases. Nature, 460:995, 2008.
[28] W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Fölling, L. Pollet, and
M. Greiner. Probing the superfluid-to-mott insulator transition at the single-atom level.
Science, 329:547, 2010.
[29] J. F. Sherson, C. Weitenberg, M. Endres, M. Cheneau, I. Bloch, and S. Kuhr. Single-atomresolved fluorescence imaging of an atomic mott insulator. Nature, 467:68, 2010.
[30] E. Toth, A. M. Rey, and P. B. Blakie. Theory of correlations between ultracold bosons released
from an optical lattice. Phys. Rev. A, 78:013627, 2008.
[31] M. Köhl, H. Moritz, T. Stöferle, K. Günter, and T. Esslinger. Fermionic atoms in a three
dimensional optical lattice: Observing fermi surfaces, dynamics, and interactions. Phys. Rev.
Lett., 94:080403, 2005.
[32] R. Jördens, N. Strohmaier, K. Gunther, H. Moritz, and T. Esslinger. A mott insulator offermionic atoms in an optical lattice. Nature, 455:204, 2008.

99[33] U. Schneider et al. Metallic and insulating phases of repulsively interacting fermions in a 3d
optical lattice. Science, 322:1520, 2008.
[34] S. Folling et al. Spatial quantum noise interferometry in expanding ultracold atom clouds.
Nature, 434:481, 2005.
[35] T. Rom et al. Free fermion antibunching in a degenerate atomic fermi gas released from an
optical lattice. Nature, 444:733, 2006.
[36] E. H. Lieb and Liniger. Exact analysis of an interacting bose gas. i. the general solution andthe ground state. Phys. Rev, 130:1605, 1963.
[37] M. Olshanii. Atomic scattering in the presence of an external confinement and a gas of impenetrable bosons. Phys. Rev. Lett., 81:938, 1998.
[38] D. M. Gangardt and G. V. Shlyapnikov. Stability and phase coherence of trapped 1d bosegases. Phys. Rev. Lett., 90:010401, 2003.
[39] T. Kinoshita, T. Wenger, and D. S. Weiss. Local pair correlations in one-dimensional bosegases. Phys. Rev. Lett., 95:190406, 2005.
[40] D. R. Hofstadter. Energy levels and wave functions of bloch electrons in rational and irrational
magnetic fields. Phys. Rev. B, 14:2239, 1976.
[41] M. Y. Azbel and M. Rubinstein. Resistance and phase diagram of quasi periodic systems.
Phys. Rev. B, 27:6530, 1983.
[42] A. N. Kolmogorov. O dinamicheskikh sistemakh s integralnym invariantom na tore. Doklady
Akad. Nauk SSSR, 93:763, 1953.
[43] I. C. Percival. Variational principle for invariant tori of fixed frequency. Journal of Physics A:
Mathematical and General, 12:L57, 1979.
[44] K. Drese and M. Holthaus. Exploring a metal-insulator transition with ultracold atoms instanding light waves? Phys. Rev. Lett., 78:2932, 1997.
[45] S. Aubry and G. Andre. Proceedings of the Israel Physical Society, volume 3. Hilger, Bristol,
1979.
[46] J. B. Sokoloff. Unusual band structure, wave functions and electrical conductance in crystals
with incommensurate periodic potentials. Physics Reports, 126:189, 1985.
[47] B. Laburthe-Tolra et al. Observation of reduced three-body recombination in a correlated 1ddegenerate bose gas. Phys. Rev. Lett., 92:190401, 2004.
[48] C. D. Fertig et al. Strongly inhibited transport of a degenerate 1d bose gas in a lattice. Phys.
Rev. Lett., 94:120403, 2005.
[49] L. Fallani, J. E. Lye, V. Guarrera, C. Fort, and M. Inguscio. Ultracold atoms in a disordered
crystal of light: Towards a bose glass. Phys. Rev. Lett., 98:130404, 2007.
[50] V. Guarrera et al. Noise correlation spectroscopy of the broken order of a mott insulatingphase. Phys. Rev. Lett., 100:250403, 2008.

100[51] D.J. Thouless. J. Phys., C5:77, 1972.
[52] S. Ostlund and R. Pandit. Renormalization-group analysis of the discrete quasiperiodic
schrodinger equation. Phys. Rev. B, 29:1394, 1984.
[53] S. Ostlund, R. Pandit, D. Rand, H. J. Schellnhuber, and E. D. Siggia. One-dimensional
schrodinger equation with an almost periodic potential. Phys. Rev. Lett., 50:1873, 1983.
[54] D. Rand, S. Ostlund, J. Sethna, and E. Siggia. Phys Rev. Lett., 49:132, 1983.
[55] S.J. Shenker and L.P. Kadanoff. J. Stat. Phys., 27:631, 1982.
[56] E.I. Dinaburg and Ya. G. Sinai. Func. Anal. Appl., 9:279, 1976.
[57] E. Ott. Chaos in Dynamical Systems. Cambridge University Press, 1993.
[58] A. M. Rey, I. I. Satija, and C. W. Clark. Hanbury brown-twiss interferometry for fractional
and integer mott phases. New Journal of Physics, 8:155, 2006.
[59] A. M. Rey, I. I. Satija, and C. W. Clark. Quantum coherence of hard-core bosons: Extended,
glassy, and mott phases. Phys. Rev. A, 73:063610, 2006.
[60] A. M. Rey, I. I. Satija, and C. W. Clark. Noise correlations of fermions and hard core bosons
in a quasi-periodic potential. Laser Physics, 17:205, 2007.
[61] T. Roscilde. Bosons in one-dimensional incommensurate superlattices.
77:063605, 2008.

Phys. Rev. A,[62] Jukka A.Ketoja and Indubala I. Satija. Self-similarity and localization. Phys. Rev. Lett.,
75:2762, 1995.
[63] A. Óttl, S. Ritter, M. Kóhl, and T. Esslinger. Correlations and counting statistics of an atom
laser. Phys. Rev. Lett, 95:090404, 2005.
[64] C. L. Hung, X. Zhang, N. Gemelke, and C. Chin. Slow mass transport and statistical evolutionof an atomic gas across the superfluidcmott-insulator transition. Phys. Rev. Lett, 104:160403,
2010.
[65] Aaron Reinhard, Jean-Felix Riou, Laura A. Zundel, David S. Weiss, Shuming Li, Ana M. Rey,and Rafael Hipolito. Self-trapping in an array of coupled 1d bose gases. Phys. Rev. Lett.,
110:033001, 2013.
[66] Anatoli Polkovnikov. Phase space representation of quantum dynamics. Annals of Physics,
325:1790–1852, 2010.
[67] U. Schollwöck. The density-matrix renormalization group. Rev. Mod. Phys, 77:259, 2005.
[68] A. J. Daley, C. Kollath, U. Schollwöck, and G. Vidal. J. Stat. Mech, page 04005, 2004.
[69] T. Kinoshita, T. Wenger, and D. S. Weiss. All-optical bose-einstein condensation using acompressible crossed dipole trap. Phys. Rev. A, 71:011602, 2005.
[70] L. D. Landau. Z. Phys. (N.Y.), 13:905, 1933.

101[71] A. L. Shluger and A. M. Stoneham. Small polarons in real crystals: concepts and problems.
J. Phys.: Condens. Matter, 5:3049, 1993.
[72] E. K. H. Salje, A. S. Alexandrov, and W. Y. Liang. Polarons and Bipolarons in High
Temperature Superconductors and Related Materials. Cambridge University Press, Cambridge, England, 1995.
[73] A. J. Millis, P. B. Littlewood, and B. I. Shraiman. Phys. Rev. Lett., 74:5144, 1995.
[74] G. Wellein and H. Fehske. Self-trapping problem of electrons or excitons in one dimension.
Phys. Rev. B, 58:6208, 1998.
[75] Z. H. Musslimani and J. Yang. J. Opt. Soc. Am. B, 21:973, 2004.
[76] M. Hartmann, F. Brandao, and M. Plenio. Laser Photon. Rev., 2:527, 2008.
[77] P. Martinoli and C. Leemann. J. Low. Temp. Phys., 118:699, 2000.
[78] G. J. Milburn, J. Corney, E. M. Wright, and D. F. Walls. Quantum dynamics of an atomic
bose-einstein condensate in a double-well potential. Phys. Rev. A, 55:4318, 1997.
[79] A. Smerzi, S. Fantoni, S. Giovanazzi, and S. R. Shenoy. Quantum coherent atomic tunnelingbetween two trapped bose-einstein condensates. Phys. Rev. Lett., 79:4950, 1997.
[80] S. Raghavan, A. Smerzi, S. Fantoni, and S. R. Shenoy. Coherent oscillations between two
weakly coupled bose-einstein condensates: Josephson effects, pi oscillations, and macroscopic
quantum self-trapping. Phys. Rev. A, 59:620, 1999.
[81] A. Trombettoni and A. Smerzi. Discrete solitons and breathers with dilute bose-einsteincondensates. Phys. Rev. Lett., 86:2353, 2001.
[82] O. Morsch, M. Cristiani, J. H. Muller, D. Ciampini, and E. Arimondo. Free expansion of abose-einstein condensate in a one-dimensional optical lattice. Phys. Rev. A, 66:021601(R),2002.
[83] A. Smerzi and A. Trombettoni. Nonlinear tight-binding approximation for bose-einstein condensates in a lattice. Phys. Rev. A, 68:023613, 2003.
[84] T. J. Alexander, E. A. Ostrovskaya, and Y. S. Kivshar. Phys. Rev. Lett., 96:040401, 2006.
[85] Ju-Kui Xue, Ai-Xia Zhang, and Jie Liu. Self-trapping of bose-einstein condensates in opticallattices:the effect of the lattice dimension. Phys. Rev. A, 77:013602, 2008.
[86] S. Wüster, B. J. Dabrowska-Wüster, and M. J. Davis. Phys. Rev. Lett, 109:080401, 2012.
[87] M. Albiez et al. Direct observation of tunneling and nonlinear self-trapping in a single bosonicjosephson junction. Phys. Rev. Lett., 95:010402, 2005.
[88] T. Anker et al. Nonlinear self-trapping of matter waves in periodic potentials. Phys. Rev.
Lett., 94:020403, 2005.
[89] M. Girardeau. Relationship between systems of impenetrable bosons and fermions in one
dimension. Journal of Mathematical Physics, 1:516, 1960.

102[90] B. Paredes, A. Widera, V. Murg, O. Mandel, S. Fölling, I. Cirac, G. V. Shlyapnikov, T. W.
Hänsch, and I. Bloch. Tonks-girardeau gas of ultracold atoms in an optical lattice. Nature,
429:277, 2004.
[91] V. Dunjko, V. Lorent, and M. Olshanii. Bosons in cigar-shaped traps: Thomas-fermi regime,
tonks-girardeau regime, and in between. Phys. Rev. Lett., 86:5413, 2001.
[92] P. B. Blakie, A. S. Bradley, M. J. Davis, R. J. Ballagh, and C. W. Gardiner. Dynamics andstatistical mechanics of ultra-cold bose gases using c-field techniques. Advances in Physcs,
57:5:363–455, 2008.
[93] A.Sinatra, C. Lobo, and Y.Castin. Rigorous results on valence-bond ground-states in antiferromagnets. J. Phys. B, 35:3599, 2002.
[94] Rafael Hipolito and Anatoli Polkovnikov. Breakdown of macroscopic quantum self-trappingin coupled mesoscopic one-dimensional bose gases. Phys. Rev. A, 81:013621, 2010.
[95] J. Mun et al. Phase diagram for a bose-einstein condensate moving in an optical lattice. Phys.
Rev. Lett., 99:150604, 2007.
[96] A. Polkovnikov and V. Gritsev. Breakdown of the adiabatic limit in low-dimensional gaplesssystems. Nat. Phys., 4:477, 2008.

Appendix A

Time-of-flight image in a single-particle expansion

Considering one atom initially trapped in the optical lattice at lattice site i, the wave functionof the atom is w(x − xi ). Assuming at time t = 0, the lattice potential is suddenly turned off, andthe atom freely expands. The wave function of the atom after some time t is,1
ψ(x, t) = √
2π

Z

Z

2

t
−i ~q
2M

1
2π

Z

2

t
iq(x−x0 )−i ~q2Mdq w(x0 − xi ),where w̃(q) is the Fourier transform of w(x), i.e., w̃(q) = √12π

e−iqx w(x)dx. Therefore,ew̃(q)e

dq =

dx

0



e

iq(x−xi )RZ~q 2 t~q 02 t10n(x, t) = |ψ(x, t)| =eiq(x−xi ) w̃(q)e−i 2M e−iq (x−xi ) w̃∗ (q 0 )ei 2M dq 0 ,
2π
Z
2
M
M 2
M
=
e−i ~t (x−xi )x1 +i 2~t x1 w(x1 )dx1 .
2π~t
2

(A.1)

(A.2)
(A.3)

Assuming the width of the Wannier function w(x − xi ) is x0 , when the expansion time is longM 2enough such that M~t xi x0  1, and ~t x0  1, we haven(x, t) ≈

M
2π~t

Z

Mx

e−i ~t x1 w(x1 )dx1

2=

M|w̃(Q)|2 ,
~t

(A.4)

with Q = M~tx . This expression represents an important approximation in the time-of-flight technique, i.e., the density distributions in the time-of-flight images actually show the momentum
distribution before the lattice is turned off.

Appendix B

Mapping from position-to-momentum space

As we have mentioned in chapter 3, ηk is connected to zn by the relation that zn = ηk e−inφwith k = Mod[nFM −1 , FM ], while θ is set to 0. Combining the above relation with Eqn. (3.10), we
can shown that the experimentally relevant variables, ηk , satisfy the following equation:
 2− ηMod[k+FM −1 ,FM ] e−iφ + ηMod[k−FM −1 ,FM ] eiφ + cos(2πk/FM )ηk = − ηk
λ
λ

(B.1)

Here we provide an example (φ = 3π/4) of the mapping between position coordinates j and quasimomentum coordinates k, represented by K(j) that we use to link the quasi-momentum-position
observables. In the plot, we highlight the Fibonacci sites. A small system size is used to make thevisualization clearer.

80

KHjL

60
40
20
0

0

20

40
j

60

80

Figure B.1: Mapping from j to k , denoted by K(j) used to determine momentum-space observables from position-space observables and vice versa. Larger points are Fibonacci sites. Here wechoose FM = 87, and K(j) is the relation that maps j to k.

Appendix C

Perturbation theory and dimerized states

We begin our analysis with the Harper equation (3.2). For λ → ∞, the single-particle wave(m)functions are localized at individual lattice sites, ψj= δj,Lm and (m) = cos(2πσLm + φ), where

Lm is defined by cos(2πσLm−1 + φ) ≤ cos(2πσLm + φ). When λ  1, we can get the single-particle
wave function through exact diagonalization, or we can use the perturbation theory to obtain theapproximate eigenvalues and eigenfunctions.
We focus first on the site Lm . Assuming that 2λ| cos(2πσLm +φ)−cos(2πσ(Lm ±1)+φ)|  1,non-degenerate perturbation theory can be used, and the results are:
(m,0)=

2λ cos(2πσLm + φ),(C.1)

(m,1)

=

0,

(C.2)

(m,2)=

1
1+,
(m,0) − 2λ cos[2πσ(Lm + 1) + φ] (m,0) − cos[2πσ(Lm − 1) + φ]

(C.3)

(m,0)

= δLm ,j , andδLm +1,jδLm −1,j+ (m,0),= − (m,0)− 2λ cos[2πσ(Lm + 1) + φ] 
− 2λ cos[2πσ(Lm − 1) + φ]

ψj

(m,1)

ψj

(C.4)
(C.5)

where (m,0) , (m,1) , (m,2) are the zeroth, first, and second order terms of the eigenvalue (m) , with
h  i
2
(m,0)
(m,1)
(m)
a similar meaning for ψjand ψj. Also, |ψj |2 = δLm ,j + O λ1
. So, up to first order of
(m) 2
1 P
= 1.
j ψj
λ,
For other sites Lm that satisfy 1 ≈ 2λ| cos(2πσLm +φ)−cos[2πσ(Lm +1)+φ]|, nondegenerate
perturbation theory doesn’t work and degenerate perturbation theory will be used. In that case,

106
the zero-order energies  and eigenfunctions ψ satisfy:





(m,0)
(m,0)
−1

 2λ cos(2πσLm + φ)
  ψLm

 ψLm


 = 
.
(m,0)
(m,0)
−1
2λ cos[2πσ(Lm + 1) + φ]
ψLm +1
ψLm +1(C.6)

The above equations result in a pair of energies that we denote as u and d :
d = λ{cos(2πσLm + φ) + cos[2πσ(Lm + 1) + φ]} −p
1 + Λ2m , and(C.7)

u = λ{cos(2πσLm + φ) + cos[2πσ(Lm + 1) + φ]} +p
1 + Λ2m ,(C.8)

where Λm = λ{cos(2πσLm + φ) − cos[2πσ(Lm + 1) + φ]}. The corresponding orthogonal eigenfunctions can be written as:




 A 
ψd = 
Br√where A =

Λ2m +1−Λm√ 2

2

of A and B: |α|2 =



 B 
and ψu = 
,
−A

(C.9)

Λ2m +1+Λm√ 2. We can choose α (β) to be the greater (lesser)
Λm +1√
1+Λ2m +|Λm |1+Λ2m −|Λm |
2 =
√
√
,
|β|
. From these results one can see that |α|2
2
2

Λm +1

√

r√



, and B =

2

1+Λm

2

2

1+Λm

increases as |Λm | increases. For example, when |Λm | = 0, |α|2 = 0.5 and when |Λm | = ∞, |α|2 = 1.
When |Λm |  1, |β|2 is comparable to |α|2 , and the states are localized at the same twoneighboring sites. Hence, they will be referred as a pair of dimerized states [ The numerical results
are shown in Fig.C.1 (a) ]. The dimerized states are found by noticing that:cos(2πσLm + φ) − cos[2πσ(Lm + 1) + φ] = 2 sin [2πσ(Lm + 1/2) + φ] sin(πσ).

(C.10)


So for the points that satisfy | sin 2πσ(Lm + 12 ) + φ |  λ1 , dimerized sites exist at Lm ,and Lm + 1. Since σ =

FM −1
FM ,

it is clear that we can construct the relation Mod[2πσLm +φ, 2π] ≈ 2π Mod [Bm /FM , 1], where m ∈ (1, 2, . . . , FM ) are indexes of increasing energy level,
Lm ∈ (1, 2, . . . , FM ) are indexes of sites’ positions, and Bm = Mod[Lm × FM −1 + l, FM ] ∈(0, 1, 2, . . . , FM − 1) are indexes that we introduce for the convenience of discussion (l is the closest
φinteger to 2π
FM ). Using this relation, it can be shown that m ≈ 2|Bm − FM /2| (we use ≈ here,

since the exact expression for even and odd FM are slightly different). The relations of m, Lm and
Bm are plotted for the case of FM = 610 in Fig. C.1 (b).

107

LmΛ
1.5
1.0
0.5
-0.5
-1.0
-1.5

100 200 300 400 500 600

m

(a)

Lm
600

Bm
600

400

400

200

200

0

200

400

m
600

0

200

400

m
600

(b)

Figure C.1: Panel (a) shows numerically that at the two red points, m1 = FM −1 and m2 = FM −2 ,Λm vanishes, so that degenerate perturbation theory is required. Panel (b) shows the relations of
Lm , Bm and m (see text) for φ = 3π
4 . The relations between Lm and m are very sensitive to thevalue of φ while the relations between Bm and m are not sensitive at all.

As we discussed, the dimerized states Lm satisfy:sin(2πσ(Lm + 1/2) + φ) ≈ 0.
Writing the above equation in terms of Bm , we get:

sin

2πBm
+ πσ
FM


≈ 0.

(C.11)

108
The solutions are

2Bm +FM −1≈ 1 or 2. Therefore,FM

Bm1 ≈ (FM − FM −1 )/2.

(C.12)

Bm2 ≈ (2FM − FM −1 )/2.

(C.13)

m1 ≈ |2(Bm1 − FM/2 )| ≈ FM −1 .

(C.14)

m2 ≈ |2(Bm2 − FM /2)| ≈ FM −2 .

(C.15)

which correspond to:

A detailed analysis shows that the band opens exactly at m1 = FM −1 and m2 = FM −2 ,regardless of FM being even or odd. This demonstrates, at perturbative level, the special behaviorof the return map at the irrational filling number FM −2 or FM −1 (See Fig.3.2).
In our calculations, we find out those Lm with the lowest values of | sin (2πσ(Lm + 1/2) + φ) |(avoiding double counting of Lm in different pairs) and then apply to those points the perturbation
theory we discussed above.
Because of the latter considerations, it is possible to demonstrate that around m1 = FM −1
or m2 = FM −2 , there exists a sequence of paired states, exhibiting the following relationships as
λ → ∞,
ψdi,τ

=

p
p
Ai δj,Lmτ −i + 1 − Ai δj,Lmτ −i +1 ,

ψui,τ

=

p
p
1 − Ai δj,Lmτ −i − Ai δj,Lmτ −i +1 ,with i  FM , and τ = 1, 2. These states determine the properties of the system when the Fermi
energy is close to the major gaps, i.e., ν = σ, or ν = 1 − σ, and thus the main properties of theband insulator phases.

Appendix D

Variational method for 1D optical lattice

In a deep 1D optical lattice, the external potential Vext is composed of a 1D lattice potential
Vlat = V0 sin2 (πx/d) and an external magnetic or dipole confinement Vdip = 12 M (ωx2 x2 + ωy2 y 2 +ωz2 z 2 ), where V0 is the lattice depth, and d is the lattice spacing. In the tight-binding approximation,
the macroscopic wave function can be written as,Ψ=Xψj (t)φj [~x, Nj (t)].

(D.1)

j

φj [~x, Nj (t)], which is centered at the minimum of the jth well, is a time-dependent function that
depends on the particle number on site j at a given instant and on the trapping frequencies alongRthe two transverse directions. φj [~x, Nj (t)] is normalized to 1, i.e., d~xφj [~x, Nj (t)]2 =1.
By using the tight-binding approximation, Eqn. (D.1), and after integrating out the spatial
degrees of freedom in the GPE, one gets the so-called discrete nonlinear Schrödinger equation(DNLS)[83],i~∂ψj= j ψj − J(ψj+1 + ψj−1 ) + µj ψj .
∂t

(D.2)

The on-site energies j arise from the external potential (along the lattice direction) superimposed on the optical lattice, i.e.,
Z
j

=

1
1
d~x φj (~x)2 M ωx2 x2 = M ωx2 d2 j 2 .
2
2

(D.3)

J is the tunneling matrix between adjacent sites,ZJ=~2 ∂ 2d~x φj (~x)
− Vlat φj±1 (~x),
2M ∂x2


(D.4)

110
and µloc
j is the “local” chemical potential,potµloc= µkin+ µintjj
j + µjZ
ZZ~2  ~ 2
1=d~x∇φj + g |ψj (t)|2 d~x |φj |4 + d~x φ2j M ωy2 y 2 + ωz2 z 2 .
2M
2

(D.5)

Depending on the relative values of the dipole-trapping frequencies and the on-site chemical
potential, the condensates in each well of the lattice can be regarded as 0D, 1D and 2D. The
effective dimensionality of the condensate gives a different scaling of the local potential with thenumber of atoms [83],αµloc
j = Uα |ψj | ,

(D.6)

4
, D = 0, 1, 2 is the dimensionality of the condensate, and Uα is a constant that does
where α = 2+D

not depend on the particle number or site index. The effective Hamiltonian of the system is [83]
Hef f

=

X

2

j |ψj |− J(ψj∗ ψj+1 + c.c) +

j


2
α+2
Uα |ψj |
.
2+α

(D.7)

To proceed further, a variational method can be used. For this method, we assume that theamplitudes ψj can be, to a good approximation, parameterized by a Gaussian function, ψj (t) =
i
h
√ q 2
2
δ(t)
2 . From this assumption, we can determine theN 4 πr(t)2 exp − (j−ξ)
+
ip(j
−
ξ)
+
i
(j
−
ξ)
2
2
r(t)
equations of motion of qi (t) = {r, δ, p, ξ}. Here N is the particle number, ξ(t) and r(t) are, respectively, the center and the width of the atom distribution, and p(t) and δ(t) are their associatedconjugate momenta (in the lattice units). For the case j = 0, the effective Hamiltonian under thisassumption becomes"Hef f = 2JN − cos pe
2 2

−σ


+

2
2+α

#
3  α
2
2 4 Λ
,
π
rα/2

(D.8)

α/2

where σ = 2r12 + r 8δ , Λ = N 2J Uα . The quasi-momentum dependence of the effective mass,2∗
m∗ −1 ≡ ∂∂pH2 = cos pe−σ , and group velocity, νg ≡ ∂H
∂p = tan p/m , allows a rich variety of dynamical

regimes. A diverging effective mass m∗ → ∞ is the signature of wave-packet self-trapping.
The variational method has been shown to be a great success in characterizing various regimesof the dynamics. It also gives a critical value of MQST-to-diffusive transition (when r0  1, where111
r0 is the initial value of r) as
α/2
Λc = r0

2+α2

3/2  
π α/4
.
2

(D.9)

Based on the variational method, when Λ > Λc (MQST regime), m∗ (t → ∞) → ∞. The
2r0= 1 − ΛΛc α .
ratio between r0 and the asymptotic width rmax (t → ∞) is given by the relation rmaxcWhen Λ ≤ Λc (the diffusive regime), m∗ (t → ∞) → ΛcΛ−Λ, r(t → ∞) → ∞. The same variational

method has also been extended into 2D and 3D optical lattices with 0D condensates [85].