Quantum Dynamics of Condensates, Atomtronic
Systems, and Photon Fluids
byBrian Thomas Seaman
B.S., Rowan University, 2003

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
2008This thesis entitled:
Quantum Dynamics of Condensates, Atomtronic Systems, and Photon Fluids
written by Brian Thomas Seaman
has been approved for the Department of Physics

Prof. Murray Holland

Prof. Victor Gurarie

Date

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

iii
Seaman, Brian Thomas (Ph.D.)
Quantum Dynamics of Condensates, Atomtronic Systems, and Photon Fluids
Thesis directed by Prof. Murray Holland

In the first part of this thesis, the dynamics of interacting ultracold bosonic atomic
gases trapped in an optical lattice are examined from the perspective of nonlinear band
theory. The mean-field Gross-Pitaevskii equation is used to model the Bloch waves forweakly and strongly interacting gases with a Kronig-Penney potential, i.e. a lattice ofdelta functions. The appearance of looped swallowtail structures in the energy bands isa highly nonlinear effect. These swallowtails are then related to period-doubled Blochstates by examining a two-color lattice. A stability analysis shows that the effective
mass of the atoms is a main feature describing the stability properties of the system.
In the second part of this thesis, the dynamics of atomtronic systems, ultracold atom analogs of semiconductor devices, is examined. Atomtronic systems are firstpresented from the perspective of the fundamental building blocks needed to create acircuit. An atomtronic battery creates a chemical potential difference. An optical lattice can play the role of a wire in electronics. The combination of N-type and P-typesemiconductors lead to an atomtronic diode. The physics behind this device are due toa transition between superfluid and insulating states for the atomic device, contrary to
the workings of an electronic diode. An atomtronic transistor is also presented.
In the third part of this thesis, the method of evaporative cooling is applied to a
photon fluid confined to a nonlinear Fabry-Perot cavity. A photon fluid is a collectionof interacting photons that exhibit fluidic hydrodynamic properties. The effects of photon recombination due to atom-mediated interactions and evaporation with an energy
dependent reflectivity of the cavity lead to Bose amplified stimulated emission into thelowest energy mode.

Dedication

To my Rachel, my wife, companion and best friend, without whom this would
never have been possible.

v

Acknowledgements

I would like to take this time to acknowledge all of the people who have helped
me to get to where I am. First, my advisor Murray Holland has given me both theindependence and the guidance necessary to stimulate my many interests in physics.
Many other people have had a hand in contributing to this thesis including Jinx Cooper,
Dana Anderson, Hong Ling, Lincoln Carr, Simon Gardiner, Meret Krämer, Dominic
Meiser, Rajiv Bhat, Satyan Bhongale, Rob Chiaramonte, Brandon Peden, Josh Milstein,Kang-Kuen Ni, Ron Pepino, Dave Tieri, and Jochen Wacter. Several others have always
been available for an engaging conversation that has helped to shape my view of theworld including Nrupen Baxi, Amar Doshi, John Ferro, Tom Joad, Glenn Koslowsky,
Anthony Leone, and John Lombardini.

viContents

Chapter
1 An Introduction to Atom-Photon Interactions

1

1.1

Atom-Photon Interactions . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.2

Bose-Einstein Condensates . . . . . . . . . . . . . . . . . . . . . . . . . .

3

1.3

Optical Lattices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.4

Photon Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.5

Overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

102 Interacting Atoms in an Optical Lattice
2.1

12Single Atom Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.1.1

AC Stark Shift . . . . . . . . . . . . . . . . . . . . . . . . . . . .

13

2.1.2

Dissipation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1.3

Optical Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1.4

Band Structure . . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

2.2

Many-Atom Second-Quantized Hamiltonian . . . . . . . . . . . . . . . .

18

2.3

Gross-Pitaevskii Equation . . . . . . . . . . . . . . . . . . . . . . . . . .

19

2.4

Bose-Hubbard Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

3 Nonlinear Band Structure

24

3.1

Constant Potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28

3.2

Potential Step . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

31

vii
3.2.1

General Solution . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

3.2.2

Particular Examples . . . . . . . . . . . . . . . . . . . . . . . . .

34

Point-like Impurity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

3.3.1

General Solution . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

3.3.2

Particular Examples . . . . . . . . . . . . . . . . . . . . . . . . .

37

Linear Limits of Transmitted and Reflected Waves . . . . . . . . . . . .

41

3.4.1

Transmitted Waves . . . . . . . . . . . . . . . . . . . . . . . . . .

41

3.4.2

Evanescent Waves . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3.5

Single Boundary Conclusions . . . . . . . . . . . . . . . . . . . . . . . .

45

3.6

A Kronig-Penney Potential and Bloch Waves . . . . . . . . . . . . . . .

49

3.7

Nonlinear Band Structure . . . . . . . . . . . . . . . . . . . . . . . . . .

52

3.8

Density Profiles of Bloch Waves . . . . . . . . . . . . . . . . . . . . . . .

57

3.9

Stability of Bloch Waves . . . . . . . . . . . . . . . . . . . . . . . . . . .

59

3.9.1

Attractive Atomic Interactions . . . . . . . . . . . . . . . . . . .

62

3.9.2

Repulsive Atomic Interactions . . . . . . . . . . . . . . . . . . . .

63

3.10 Analytical Methods in Nonlinear Band Theory . . . . . . . . . . . . . .

66

3.10.1 Solution by Cancellation . . . . . . . . . . . . . . . . . . . . . . .

67

3.10.2 Solution by Three Mode Approximation . . . . . . . . . . . . . .

70

3.10.3 Solution via a Piece-wise Constant Lattice . . . . . . . . . . . . .

71

3.11 Periodic Potential Conclusions . . . . . . . . . . . . . . . . . . . . . . .

73

3.12 Two Color Lattice . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

75

3.13 The Two Color Lattice and Formation of Swallowtails . . . . . . . . . .

77

3.14 Stability Properties and Comparison with Experiments . . . . . . . . . .

83

3.15 Two-Color Lattice Conclusion . . . . . . . . . . . . . . . . . . . . . . . .

87

3.3

3.4

4 Atomtronics
4.1

Quantum Phase Transition in the Bose-Hubbard Formalism . . . . . . .

89
93

viii
4.2

Atomtronic Battery. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

96

4.3

Atomtronic Conductors . . . . . . . . . . . . . . . . . . . . . . . . . . .

98

4.3.1

Doped Materials . . . . . . . . . . . . . . . . . . . . . . . . . . .

99

4.3.2

Material Current Properties . . . . . . . . . . . . . . . . . . . . . 101

4.4

Atomtronic Diode. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

4.4.1

Diode PN-junction configuration . . . . . . . . . . . . . . . . . . 104

4.4.2

Diode current-voltage characteristics . . . . . . . . . . . . . . . . 106

4.5

Atomtronic Transistor . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

4.6

Atomtronics Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . 112

5 Photons Interacting in a Cavity

116

5.1

Photon Modes of the Cavity . . . . . . . . . . . . . . . . . . . . . . . . . 117

5.2

Quantum Hamiltonian . . . . . . . . . . . . . . . . . . . . . . . . . . . . 117

5.3

5.2.1

Polarization of Atoms . . . . . . . . . . . . . . . . . . . . . . . . 118

5.2.2

Fourth Order Hamiltonian . . . . . . . . . . . . . . . . . . . . . . 121

Photon Fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 123

6 Evaporative Cooling of a Photon Fluid to Quantum Degeneracy

127

6.1

System Dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128

6.2

Photon Number Distribution . . . . . . . . . . . . . . . . . . . . . . . . 132

6.3

Photon Spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 135

6.4

Summary and Outlook . . . . . . . . . . . . . . . . . . . . . . . . . . . . 138

Appendix
A Jacobian Elliptic Functions

140

B Completeness of Constant Potential Solution Set

142

ix
C Atomtronics Calculational Details

144

C.1 Small System Considerations . . . . . . . . . . . . . . . . . . . . . . . . 144
C.2 Battery Contacts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 146
C.3 Maximum current solution . . . . . . . . . . . . . . . . . . . . . . . . . . 147
C.4 Current Response Calculation . . . . . . . . . . . . . . . . . . . . . . . . 147
C.5 Atomtronic Wires and Diodes . . . . . . . . . . . . . . . . . . . . . . . . 148
C.6 Transistor . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 149

Bibliography

151

x

FiguresFigure1.1

Velocity distribution showing formation of a Bose-Einstein condensate. .

4

1.2

Velocity distributions showing the Mott-insulator to superfluid transition.

8

2.1

Band structure of atoms in a one dimensional lattice. . . . . . . . . . . .

18

3.1

Transmitted states to the nonlinear Schrödinger equation with a potential
step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3.2

35

Evanescent states to the nonlinear Schrödinger equation with a potential
step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

36

3.3

Localized states of a repulsive condensate in the presence of an impurity.

40

3.4

Localized states of an attractive condensate in the presence of an impurity. 41

3.5

Nonlocalized condensates near an impurity . . . . . . . . . . . . . . . .

42

3.6

Nonsymmetric condensates in the presence of an impurity. . . . . . . . .

43

3.7

Band structure of a weakly repulsive condensate. . . . . . . . . . . . . .

53

3.8

Band structure of a strongly repulsive condensate. . . . . . . . . . . . .

54

3.9

Band structure of for systems with different interaction strengths to describe the presence of swallowtails. . . . . . . . . . . . . . . . . . . . . .

55

3.10 Band structure of a weakly attractive condensate. . . . . . . . . . . . . .

56

3.11 Band structure of a strongly attractive condensate. . . . . . . . . . . . .

57

3.12 Density profiles for an attractive condensate. . . . . . . . . . . . . . . .

59

xi
3.13 Density profiles for a repulsive condensate. . . . . . . . . . . . . . . . . .

60

3.14 Instability spectrum of the condensate. . . . . . . . . . . . . . . . . . . .

64

3.15 Growth of instability fo the condensate. . . . . . . . . . . . . . . . . . .

65

3.16 Instability time of the condensate. . . . . . . . . . . . . . . . . . . . . .

66

3.17 Band structure calculated from different analytic methods.

. . . . . . .

69

3.18 Sketch of two color lattices. . . . . . . . . . . . . . . . . . . . . . . . . .

77

3.19 Band structure for a weakly repulsive condensate in a two-color lattice.

79

3.20 Band structure of a strongly repulsive condensate in a two-color lattice.

81

3.21 Condensate density profiles. . . . . . . . . . . . . . . . . . . . . . . . . .

82

3.22 Instability time in a one-color lattice. . . . . . . . . . . . . . . . . . . . .

84

3.23 Quasimomentum at which instability begins. . . . . . . . . . . . . . . .

86

4.1

Atomtronic energy band structure. . . . . . . . . . . . . . . . . . . . . .

92

4.2

Zero-temperature phase diagram of a condensate in an optical lattice. .

94

4.3

Schematic of atoms in a lattice connected to an atomtronic battery.

. .

98

4.4

Schematics of analogies of doped semiconductors. . . . . . . . . . . . . .

99

4.5

Doped semiconductors in the phase diagram. . . . . . . . . . . . . . . . 101

4.6

Current curves for various materials. . . . . . . . . . . . . . . . . . . . . 104

4.7

Schematics of PN-junction diodes. . . . . . . . . . . . . . . . . . . . . . 105

4.8

Phase diagram of PN-junction configuration. . . . . . . . . . . . . . . . 107

4.9

Current-voltage curve of an atomtronic diode. . . . . . . . . . . . . . . . 109

4.10 Schematic of an atomtronic bipolar junction transistor. . . . . . . . . . . 110
4.11 Phase diagram of an atomtronic transistor. . . . . . . . . . . . . . . . . 111
5.1

Bogoliubov dispersion of a photon fluid. . . . . . . . . . . . . . . . . . . 125

6.1

Photon fluid processes in a Fabry-Perot cavity. . . . . . . . . . . . . . . 129

6.2

Photon number distribution of lowest energy mode. . . . . . . . . . . . . 134

xii
6.3

Spectrum with vanishing collision rates. . . . . . . . . . . . . . . . . . . 137

6.4

Spectrum with different collision rates. . . . . . . . . . . . . . . . . . . . 139

Chapter 1

An Introduction to Atom-Photon Interactions

With the creation of a Bose-Einstein condensate with a trapped atomic gas in1995 [6, 39, 19], a new and diverse field of physics reliant on the ultracold properties ofbosons was born. The quantum statistics of bosons lead to the occupation of a single
quantum state at ultracold temperatures. Many properties, such as superfluidity andphase coherence, become experimentally accessible with these new systems. In this thesis, we examine three areas of this field of ultracold bosonic systems: the nonlinear band
structure of interacting atoms confined to an optical lattice potential, the creation of
atomtronic systems analogous to electronic semiconductor devices, and the evaporative
cooling of a photon fluid to quantum degeneracy. These three works are connected byseveral over-arching themes. The fundamental particles under examination are ultracold bosons. Each system is inherently dynamical, whether considering the motion of
atoms or of photons. And finally, the ability to control the interactions between atoms
and light make all three systems possible.
This introduction covers several of the important topics necessary to analyze
the three areas of examination. The use of atom and photon interactions to control
each of these bosonic species is discussed. The properties and creation of a BoseEinstein condensate is then examined, followed by an analysis of a trapping opticallattice potential. The system of a photon fluid is then introduced.

2

1.1

Atom-Photon Interactions
The interactions between atoms and light are fundamental to nearly all fields ofphysics. It is common to view these interactions through one of two lenses. Either
light alters the properties of the atoms or atoms alter the properties of the light. There
are several ways of examining the atoms and photons: classically, semiclassically, or
quantum mechanically. The system under consideration usually dictates the appropriate
method of examination.
For many physical systems, the fields of light do not need to be considered in fullquantum mechanical detail and can be instead accurately described as classical fields.
Classical fields are characterized solely by Maxwell’s equations with the appropriate
boundary conditions and this description is generally a valid approximation [83]. The
presence of atoms then allows for the light to behave in new ways due to the atoms’
polarization or magnetization. For instance, a nonlinear medium instigates the process
of four-wave mixing for light driven by strong pump fields [110]. Two pump fields of thesame frequency but opposite direction drive a nonlinear medium and with the aid of a
weak probe field create a sideband component. The laser is another excellent exampleof how atoms can be used to enhance the properties of photons [148]. An active gain
medium inside a mirrored cavity enhances the effects of stimulated emission and createsa high intensity coherent beam of light.
A fully quantized view of light is necessary to understand many phenomena. For
instance, the details of spontaneous and stimulated emission depend crucially on thequantum nature of light [110]. Emission is the process of creation of a photon by lossof energy by an atom. Stimulated emission is the enhancement of photon creation ofa particular mode, due to the presence of additional photons in the same mode, withthe same phase. The bosonic nature of photons play a large role in this effect. The
quantum nature also frequently appears in the noise profile of the light in phenomenon

3like resonance fluorescence [91, 90] and squeezing [159].
Atom-field interactions can also lead to mechanical effects on atoms such as the
Doppler cooling of atoms [100]. In laser cooling, the moving atoms see a Doppler shifted
field of light. For a red-detuned laser, atoms moving against the light’s direction will
have an enhanced interaction strength. Laser cooling is one of the traditional methods
used to cool atomic gases into quantum degeneracy.
The creation of an optical lattice is another example of how light can be used tocontrol the motion of atoms [62]. Counter-propagating lasers create a sinusoidal standing wave intensity pattern and due to the AC Stark shift have the effect of establishingan external periodic potential for the atoms. The use of an optical lattice has led to theexperimental verification of the Mott-insulator superfluid transition [85, 64].
This thesis examines the nonlinear interactions of light and atoms from both
perspectives. In the first two systems, atom-atom interactions are enhanced by thepresence of an optical lattice. The effects of increased interactions on the energy band
structure is examined and new features are found. The interactions are also used tocreate atomic analogs of electronic devices replicating the motion of electrons in a solid.
In the final system, a nonlinear atomic medium inserted in a Fabry-Perot cavity createsan effective interaction between photons. This system can reach quantum degeneracy
via a process analogous to evaporative cooling of atomic gases.

1.2

Bose-Einstein Condensates
The observation of Bose-Einstein condensation (BEC) of alkali-metal atoms in1995 [6, 39, 19] confirmed predictions that Bose and Einstein first developed for noninteracting bosons in 1924 [16, 47]. A Bose-Einstein condensate is created when the de
Broglie wavelength of the bosons is larger than the average interparticle spacing, leading
to a macroscopic occupation of a single quantum state. This is a low temperature phasetransition that is only accessible to bosons and not fermions, due to the differences in4

Figure 1.1: The velocity distribution of ultracold rubidium atoms after an expansion
performed at JILA [6]. From left to right, the temperature of the system is decreased
from 400nK to 50nK, and a condensate - a macroscopic population of the ground state- appears. The left most frame has a negligible condensate fraction since it is just above
the transition temperature. In the right most expansion, nearly all of the atoms arecondensed, corresponding to the sharp peak.

their quantum statistics.
In atomic systems, Bose-Einstein condensation was first observed for rubidium,
sodium, and lithium atoms but has later gone on to include many other atomic species.
The atoms are first collected in a magneto-optical trap and then cooled to micro-Kelvintemperatures [100]. Evaporative cooling in a magnetic trap leads to temperatures in
the nano-Kelvin range. Condensed atoms can number well into the millions.
After cooling, the velocity distribution of the atoms can be examined by performing a time-of-flight analysis. Figure 1.1 presents the velocity distributions for systems
with three different temperatures, just above the transition temperature at 400nK and
two below the transition temperature at 200nK and 50nK, as presented by the team at
JILA [6]. At the lowest temperature, nearly all of the atoms have been condensed andthe strong peak corresponds approximately to the velocity distribution of the harmonic
oscillator ground state.

5
Before the examination of atomic gases, superfluid liquid Helium was the firstsystem that was thought to exhibit properties of a Bose-Einstein condensate [87, 3, 103].
As a liquid, the effects of interactions in the superfluid Helium are considerable. This
is in contrast to the new atomic gas systems, in which interaction effects can be quite
small. This property of weak interactions has allowed for the control of many of thesystem parameters and examination of interesting new phenomena. Since the atoms areat very low temperatures, a single parameter, the s-wave scattering length, can usually
be used to characterize inter-atomic interactions. In several current BEC systems, theinteractions can be made repulsion and attraction using magnetically tunable Feshbachresonances [81]. A Feshbach resonance occurs when the energy of a bound state of theinteratomic potential equals the kinetic energy of two incoming atoms. The relative
energy of the bound state can be controlled by an external magnetic field and has theeffect of dramatically increasing the s-wave scattering length. The sign of the scattering
length can also be altered.
Bose-Einstein condensates exhibit several interesting properties. For instance,
the gases are superfluid and sustain quantized vortices [65]. A superfluid has zeroviscosity and sustains frictionless flow. The formation of quantized vortices is anothermanifestation of superfluidity. Below a critical rotation rate, the condensate will remain
stationary. However, a faster rotation causes the superfluid to rotate around vortex
cores. Each vortex core can be described by discrete amounts of angular momentum.
Current work is being performed trying to create Bose-Einstein condensation in
new systems. Molecular BECs have been created by tuning a Feshbach resonance to create ultracold molecules from condensed atoms [44]. Quantum degenerate fermions havebeen coaxed into pairing in ways similar to Cooper pairing [131]. Since the fermions
are at very low temperatures, the Cooper pairs are also created in a quantum degenerate state. The pairs, however, are effectively bosons and so they create a “fermioniccondensate”.

6
A Bose-Einstein condensate, confined to an optical lattice as described in thefollowing section, will be the primary system of interest for the first two parts of thisthesis. The final part describes how methods traditionally associated with atomic condensates, such as evaporative cooling, can be used to create a quantum degenerate setof photons, a condensed and superfluid state of light.

1.3

Optical Lattices
The use of lasers has led to the creation of many new tools to control, manipulate

and examine atoms. One such useful tool is an optical lattice that has been developed
to examine and alter the properties of ultracold atomic gases [82]. An optical latticeis created by aligning counter-propagating lasers to form a sinusoidal intensity profile.
This intensity pattern creates a sinusoidal AC Stark shift that acts as an external
periodic potential for the atoms. Since the parameters of the laser, such as intensityand polarization, are easily tuned, the optical lattices that are formed can be easily
manipulated [62]. Atoms are generally loaded into an optical lattice by first forminga Bose-Einstein condensate using traditional methods. The lattice is then created byadiabatically turning on and increasing the laser intensity.
Bose-Einstein condensates of atomic gases are generally weakly interacting since
they are created at very low densities. Therefore, the effects of interactions are usually
not very prominent. An optical lattice can be used to convert the weakly interacting
condensate into a strongly correlated gas by simply increasing the laser intensity. This
causes a greater localization of the atoms into the lattice wells. The greater localization
is able to increase the effects of interactions compared with the kinetic energy associatedwith moving from one lattice site to the next.
Ultracold atoms in an optical lattice exhibit many interesting properties. They
are particles in a periodic potential and so exhibit a gapped energy band structure. Theenergy structure of a many-particle system can also exhibit a gapped band structure as it

7
is the condensate wave function that is in the presence of an effective potential due to thelattice and interactions [105]. For large interactions, new features are shown to appear,such as the occurrence of looped swallowtail structures in the energy band [153, 154].
Ultracold atoms in optical lattices allow for the examination of both attractive andrepulsive interactions due to the alterability of the sign of the s-wave scattering lengthwith a Feshbach resonance.
The effects of strong interactions and strong correlations present themselves inother ways [85]. For instance, ultracold atoms exhibit a quantum phase transition as thelattice depth is increased. For a shallow lattice, the effects of the potential are negligible
and the atoms remain a condensate with a well defined system-wide phase. This is thesuperfluid phase of the gas. As the depth is increased, the atoms become localized inseparate wells. For deep enough lattices, the phase in each well is independent of theothers. In addition, the hopping of atoms from one site to the next becomes severely
restricted and the Mott-insulator phase is formed. When the atoms are in a Mott
insulator state, each well contains on average an integer number of atoms. The number
variance on each site decreases quickly as the well depth increases.
The Mott insulator to superfluid transition was clearly demonstrated in an experiment by Greiner, et. al. [64], in which the interference pattern of atoms trapped in an
optical lattice shows a transition from a phase coherent state to a state with no global
phase coherence, as seen in Fig. 1.2. The figure presents eight different lattice depths,
increasing from 0ER to 20ER , where ER is the recoil energy of the photons in the lattice.
The interference pattern shows an initially coherent system lose its phase coherence as
the superfluid transitions into a Mott insulator. Another feature of the phase transitionis that, for a system with a small hopping energy, there is a jump in chemical potentialof a system with slightly less than integer filling to a state with slightly higher thaninteger filling. The gap is on the order of the interaction energy and leads to regions
where the system is unaffected by changes in the chemical potential. This transition is

8

a

b

c

d
1

e

f

g

h
0

Figure 1.2: The velocity distribution of ultracold rubidium atoms trapped in an opticallattice after an expansion performed by Greiner, et. al. [64]. From (a) to (h), the lattice
depth is increased from 0ER to 20ER , where ER is the atomic recoil energy due to thelattice lasers. Initially, the atoms all have a similar momentum and are a collectively
a superfluid state. As the lattice depth is increased, additional momentum peaks form
due to the interference patterns of partially localized atoms. A single phase still defines
the entire system. Further increases in the lattice depth eventually lead to a systemwith no phase coherence.

an integral component to the workings of the atomtronic devices that will be discussed.

1.4

Photon Fluid
A photon fluid is a way to describe photons with atom-mediated interactions [33]and is in general a system of photons that exhibit hydrodynamic properties. Atom
mediated photon-photon interactions can be obtained for light trapped between twoparallel mirrors placed close together in the presence of a nonlinear medium. The
nonlinear polarizability creates the higher order interaction effect. A linear polarizability
creates only an energy shift of the bare modes. The nonlinear medium can be any ofa number of substances. For instance, a nonlinear crystal deposited onto the mirrors
or a Rydberg gas could be used. The sign of the interactions depends on whether thecavity mode is red or blue shifted relative to the closest atomic transition. The mirrors

9forming the cavity cannot be perfectly reflecting so the system must be continuously
pumped. The pumping rate and loss rates set the number of photons in the system.
A photon fluid, or weakly interacting photon gas in the nonlinear Fabry-Perotcavity, is an effectively two dimensional system if excitations to higher longitudinalmodes are minimized. The reduced dimensionality and small transverse momenta leadto an effective mass, m ≡ h̄ω/c2 , determined by the frequency of the lowest longitudinalmode. This effective mass is due to the tight confinement in one direction and therelatively free movement in the other directions and will be presented in more detailin Chapter 5. A photon fluid can have interesting properties, especially at low temperatures. Although the system is composed of photons, the chemical potential does
not vanish due to the interactions. The energy required to put a photon in the cavity
depends on the number of photons already in the cavity. This is in stark contrast to the
traditional Planck blackbody photon system in which the chemical potential vanishes.
For low temperatures, since the system is composed of weakly interacting bosons,the photons should populate the lowest energy mode with a macroscopic occupation.
The number of photons in this mode will be determined by the intensity of the pumping.
The excitation spectrum for the ground state of the fluid is given by the Bogoliubov
dispersion relation [100]. For small momenta, the dispersion relation is linear. Thisimplies that elementary excitations are phonon modes, i.e. that sound-like waves canpropagate through the photon fluid. The presence of phonon modes means that different long wavelength excitations can propagate together without distortion. For larger
momenta, the dispersion relation becomes quadratic. These excitations behave like freenonrelativistic particles with a finite mass. For a completely non-interacting system,
the dispersion relation does not have the linear component and is instead only quadratic
in nature. The interactions are critical to the linear aspect of the dispersion relation.
The Bogoliubov form of the dispersion relation implies that the system is superfluid andshould support dissipationless flow and quantized vortices [38].

10
In this thesis, it is shown how a quantum degenerate photon fluid can be createdusing the method of evaporative cooling in which the highest energy photons leave thesystem.

1.5

Overview
Chapter 2 presents the theory of an optical lattice. It is shown how the AC Stark

shift due to the interference pattern of counter-propagating electromagnetic fields cancreate a periodic lattice potential. The energy band structure of non-interacting atoms is
described. The many-atom second-quantized Hamiltonian is introduced and two models
or approaches, the Bose-Hubbard Hamiltonian and the Gross-Pitaevskii equation, arepresented.
Chapter 3 examines the properties of interacting atoms confined to a periodic
potential. This study uses the mean-field Gross-Pitaevskii equation with a KronigPenney potential to model the system. The complete set of stationary solutions tothe single boundary problem, step potential or localized delta function potential arepresented. These states are the nonlinear equivalents to transmitted and reflected waves.
The Bloch states calculated from the Kronig-Penney potential are compared to othermethods for calculating stationary states in lattices. The presence of swallowtails isthen connected to period-doubled states of a two-color lattice.
Chapter 4 presents the concept of atomtronics, in which ultracold atom analogs
play the role of traditional semiconductor devices. The tight binding Bose-Hubbard
model is used to model the dynamics of the atoms as they flow through the lattice configurations. The quantum phase transition between the Mott insulator and superfluid
phase is a key component to the workings of the atomtronic devices. Key electronic
analogs are presented, such as the connection between electronic charge and presenceof atoms. Atomtronic semiconductors are connected to form an atomtronic diode andtransistor.

11
Chapter 5 presents a Hamiltonian that describes the atom-mediated photonphoton interactions that occur in a Fabry-Perot cavity filled with a nonlinear medium.
The presence of the photons polarizes the nonlinear medium, which then has a backeffect on the photons. The interactions lead to the concept of a photon fluid in which
the presence of Bogoliubov excitations implies the existence of a superfluid state of light.
Chapter 6 presents a method for evaporatively cooling a photon fluid to quantum
degeneracy. The five key processes that occur in the system are described in detail.
The cavity is constantly pumped by an incoherent broad spectrum source. Photons
trapped in the cavity slowly decay through the mirrors. Photons self-interact to createa mean field energy shift. Photons recombine to form new photon modes. Finally, highenergy photons quickly leave the cavity either out the cavity sides or due to a frequency
dependent reflectivity of the mirrors. The photons quickly populate the ground statemode and possess a nearly Poissonian number distribution. The spectral width of thesystem is analyzed to show a narrow spectrum and long coherence time.

Chapter 2

Interacting Atoms in an Optical Lattice

In current experiments with Bose-Einstein condensates, the densities, ρ, are generally low compared to interaction scattering lengths, as , ρa3s " 1, so that the system
is weakly interacting. Many interesting properties of condensates, however, occur forlarge interactions and strong correlations. There are several ways in which the effects
of interactions can be enhanced. The use of a Feshbach resonance, in which an external
magnetic field can be used to tune the two-body interactions, has provided many new
avenues of exploration [81]. The confinement of atoms in an optical lattice, formedby the AC Stark shift due to counter-propagating lasers, has also allowed for the effects of interactions to be enhanced [85]. The quantum phase transition between the
Mott insulator and superfluid phases observed in optical lattices has been an impressive
demonstration of the effects of strong interactions [52, 64]. Atoms are usually loaded
into optical lattices by first condensing the atomic system into a condensate and thenadiabatically increasing the intensity of the lasers that form the optical lattice. Opticallattices have also been used to trap and cool atoms and ions [75, 66]. The use of ionsin an optical lattice present an interesting opportunity to create a quantum registry fora quantum computer.
Atoms confined in optical lattices open many possibilities for exploration. The
similarity of cold atoms in an optical lattice to electrons in a solid can be leveraged to
explore phenomena associated with condensed matter systems. Optical lattices have13
several advantages over electronic systems. The structure and geometry of an opticallattice can be very well controlled and are largely free of unintended defects. These
properties are utilized in the following two chapters to describe the nonlinear band
structure of atoms in a lattice and to reinterpret electronic systems in the context of
ultracold atoms.
This chapter first presents the Hamiltonian of a single atom in the presence of a
monochromatic field and introduces the origin of an optical lattice. The energy structureof the atom is presented. This Hamiltonian is then extended to a many-body system anda second-quantized Hamiltonian is presented. It is difficult to directly perform calculations with the second-quantized Hamiltonian, and so several approaches are commonly
taken. This chapter then presents two models that are used to describe the properties,
both stationary and dynamic, of a dilute gas of atoms confined in an optical lattice. The
Gross-Pitaevskii equation, a mean field theory that well describes the condensate wave
function for large particle numbers, is presented, as well as the Bose-Hubbard model
which describes atoms in well separated lattices sites for large lattice depths.

2.1

Single Atom Hamiltonian
An optical lattice is created by interactions between atoms and light. A gasof neutral atoms, traditionally alkali-metal due to ease of manipulation, loaded into
sets of counter-propagating laser beams experience a sinusoidal potential. The primarymechanism controlling the potential is the AC Stark shift between a single atom andthe field of light. This section explores the Hamiltonian of a single atom interacting via
the optical dipole potential with a monochromatic field of a large number of photons.

2.1.1

AC Stark Shift

The atom is considered a two-level system with a dipole moment between theinternal states and is confined to a fixed volume V . The internal atomic states are labeled

14|g# and |e# for the ground and excited states, with energies h̄ωg and h̄ωe , respectively.
The electric field will be considered classically, and given by a monochromatic planewave,# r , t) = #εE(#r ) cos(ωt) ,
E(#

(2.1)

where #ε is the polarization of the field, E(#r ) contains the spatial dependence of the fieldand ω is the frequency of the field. The Hamiltonian that describes the internal statesof the atom, with a position of #r, and the electric field is given byĤ0 = h̄ωg |g#$g| + h̄ωe |e#$e| − dE(#r ) cos(ωt) ,(2.2)

where d ≡ $e|er̂|g# · #ε is the dipole moment of the atom in the direction of the electricfield. Note that the dipole approximation has been employed which is valid as long as
the wavelength of the field is much larger than the size of the atom. The frequency
difference between the internal states is defined as ω0 ≡ ωe − ωg . The detuning of thephoton frequency with the atom energy difference is given by ∆ ≡ ω − ω0 .
When the photons are near resonant with the atom energy difference, |∆| " ω0 ,the rotating wave approximation can be used to eliminate off-resonant interactions. The
Hamiltonian in the interaction picture becomesĤI ≈h̄Ω(#r)h̄Ω∗ (#r )|e#$g|ei∆t +|g#$e|e−i∆t ,
2
2

(2.3)

where the Rabi frequency is given by,Ω(#r ) =dE(#r )
.
h̄

(2.4)

The lack of a permanent electric dipole moment means that there is no first orderenergy shift. However, the field can induce a dipole, which in term interacts with thelight to create the energy shift. If the detuning is large compared to the Rabi frequency,|Ω| " |∆|, then second order perturbation theory can be used to determine the AC
Stark shift,∆E = ±|Ω(x)|2|$e|HI |g#|2
= ±h̄
.
h̄∆
4∆

(2.5)

15
The plus sign of the energy shift corresponds to the energy change of the ground state,while the minus refers to the excited state energy shift. It is the spatial dependence of
the Rabi frequency, and hence the spatial intensity of the light, that allows for a spatiallydependent external potential. This is utilized to create an optical lattice. Note thatfor blue detuning, ∆ > 0, the ground state atoms will prefer to occupy a minimum inthe laser intensity, while for red detuning, ∆ < 0, the atoms will prefer the maximum
intensity.

2.1.2

Dissipation

In addition to the conservative processes that lead to the AC Stark shift, additional dissipative processes can also occur. Dissipation can be caused by spontaneousemission from the excited state. For a lattice with many atoms, interactions can also
cause decay to the ground state. The effective spontaneous emission rate, due to partial
occupation of the excited state, is given by
Γef f =|Ω(x)|2 Γse
,
8∆2

(2.6)

where Γse is the linewidth of atoms in the excited state. The effects of dissipation are
therefore mitigated when the population of the excited state is kept to a minimum,and so a blue shifted lattice is preferred since this configuration places atoms in the
low intensity portions of the lattice. Since the dissipation scales as ∆−2 and the latticepotential scales as ∆−1 , an off-resonant, but high intensity, laser optimizes the effectsof the lattice while minimizing the effects of dissipation.

2.1.3

Optical Lattice

An optical lattice is created by aligning sets of counter propagating laser beams.
The beams create a standing wave interference pattern with sinusoidal intensity and16Rabi frequency. In one dimension, the Rabi frequency becomes,|Ω(x)| = 2Ω0 sin(kx) ,(2.7)

where k = 2π/λ is the wavenumber of the light and λ is the wavelength of the light.
This leads to a corresponding external potential,
V (x) =h̄Ω20
sin2 (kx) .
∆

(2.8)

Two and three dimensional lattices can be created by using multiple sets of counter
propagating lasers. In addition to the intensities, the relative polarizations of the beams
can also influence the structure of the optical lattice. It is possible to create triangular,
trapezoidal and rectangular lattice configurations.
A single set of counter-propagating beams can be used to create a series of flat
pancake-like structures in which atoms may easily move in two dimensions but are
restricted in the third. The use of two sets of beams can create a set of interacting
tube-like structures in which only one direction allows easy motion. If a third set of
weak intensity beams is added, an effectively one-dimensional system with a periodic
potential can be created. If the beam intensity is increased, a true three dimensionallattice can be created.

2.1.4

Band Structure

One of the many important features of particles in a periodic potential is thatthe single-particle energy spectrum takes the form of a band structure where bands ofallowable energy are separated by finite gaps. This is a feature common to all systems inwhich particles move in a periodic potential. There is a well-developed theory describing
these effects [9].
In the case of a one dimensional system, it is possible to solve for the stationary
states analytically. The Hamiltonian of an atom in an optical lattice is given by
Ĥ =

p̂2
+ V0 sin2 (kx) ,
2m

(2.9)

17
where the kinetic energy of the atom has been specifically included and the periodicityof the lattice is a ≡ π/k = λ/2. According to Bloch’s theorem [9], the eigenstates ofsuch systems can be written as the product of a plane wave and a function periodic in
the lattice spacing, as follows:iqx (n)φ(n)uq (x) ,q (x) = e(2.10)(n)u(n)
q (x) = uq (x + a) .

(2.11)

where

The stationary Schrödinger equation, Ĥφ = Eφ, then yields!

"(p̂ + h̄q)2(n)+ V0 sin2 (kx) u(n)q (x) = Euq (x) .
2m(2.12)

This leads to the interpretation of h̄q as the quasi-momentum of the atom since it isassociated with the discrete translational symmetry of the lattice just as momentum isassociated with translational symmetry of a free particle. The Schrödinger equation is
simply a form of the Mathieu equation and the Bloch waves are, therefore, given by
Mathieu functions [1].
The eigen-energy structure exhibits periodicity with respect to the quasi-momentum
given by the reciprocal of the lattice spacing, K = 2π/a. According to Bloch’s theorem,
(n)

(n)

(n)

states given by φq (x) and φq+K (x) have the same energy Eq . Therefore, it is only
necessary to examine the energy diagram in the restricted quasi-momentum domain of−π/a < q < π/a. All other quasi-momenta map back into this restricted region.
Figure 2.1 presents the band structure for various lattice depths of a one dimensional lattice and tight confinement in the remaining two directions. For a vanishingly
small lattice, the atoms possess the traditional quadratic free particle spectrum. As thelattice depth is increased, a gap appears between bands. Initially the gap at the edgesof the Brillouin zone of the first band is given by V0 /2. The gap increases as the lattice
depth increases and the energy band width shrinks. Interesting new features will appear
for interacting systems.

18

V=1E

V=0E

R

E/ER

15

R

15

15

10

10

10

5

5

5

0

0

0

−1

0

1

qa/π

−1

0

qa/π

V=10ER

−1

1

0

qa/π

1

Figure 2.1: Band structure for an atom confined to a one dimensional optical lattice.
For small lattice depths the energy structure has a nearly free-particle quadratic spectrum. Increased lattice depth creates gaps between bands, which themselves become
increasingly flat. Note that ER ≡ h̄2 π 2 /2ma2 .

2.2

Many-Atom Second-Quantized Hamiltonian
The AC Stark shift is an effect between a single atom and the field of light.

There are also the interactions between the atoms themselves. The formalism of second
quantization is useful to describe these features. The second-quantized many-body
Hamiltonian in grand canonical ensemble is given by
Ĥ =

#

h̄2 2
dxΦ̂ (x) −∇ + VL (x) Φ̂(x)2m(2.13)+#dxΦ̂† (x)(VE (x) − µ)Φ̂(x)(2.14)#(2.15)

+

†

1
2

!

"

dxdx# Φ̂† (x# )Φ̂† (x)VI (x − x# )Φ̂(x)Φ̂(x# ) ,where Φ̂(x) is the boson field annihilation operator for a boson at position x, VL (x) is thelattice potential, VE (x) is an arbitrary external potential, µ is the chemical potential,
and VI (x − x# ) is the two-body interaction potential of atoms at positions x and x# . The

19
two-body interaction can be accurately described by an effective contact potential [38],VI (x − x# ) ≈4πah̄2δ(x − x# ) ,m(2.16)

where a is the s-wave scattering length and m is the mass of the atoms. The assumptionof a contact interaction characterized solely by the s-wave scattering length is valid
as long as the temperature of the atoms is low enough to ignore p-wave and higher
interaction terms. In this regime, the specifics of the interatomic potential are onlya small perturbation to the s-wave scattering length. The interaction can be either
repulsive or attractive depending on the sign of the scattering length. As described
previously, the interaction strength and sign can be altered by the use of a Feshbachresonance.

2.3

Gross-Pitaevskii Equation
Exact calculations with the full second-quantized Hamiltonian can be impracti-

cable for more than a small number of atoms. In cases where a large number of atoms
are present, a mean-field theory can be used to describe the system using physicallymeaningful quantities as parameters [38].
The atomic field annihilation operator can be expanded exactly asΦ̂(x) =$

Φi (x)âi ,(2.17)iwhere the set {Φi (x)} is a complete basis of single particle wave function, and {âi }are the corresponding bosonic annihilation operators. The annihilation operators aredefined such thatâi |Ni # =%Ni |Ni − 1# ,(2.18)

where Ni are the eigenvalues of the â†i âi operator and give the number of atoms in thei state. The annihilation operators obey the usual Bosonic commutation relations. In a
Bose-Einstein condensate, most of the atoms are in a single state, so that N0 ( Ni for20
all i )= 0. Therefore, N0 ≈ N0 ± 1 and the approximation â0 ≈ â†0 ≈√N0 can be made.

For a uniform gas confined to a volume V , the wave function can then be expressed asΦ̂(x) =&

N0 /V + δΦ̂(x) .

(2.19)

The correction δΦ̂(x) corresponds to excitations of the system out of the condensate
mode. With a nonuniform system, the field operator can be expanded in a more generalform asΦ̂(x) = φ(x) + δΦ̂(x)(2.20)

where φ(x) ≡ $Φ̂(x)# is a complex number called the condensate wave function, ρ(x) =|φ(x)|2 is the density of condensed atoms at position x, and $δΦ̂(x)# = 0. The condensate
function also has a well-defined phase. Note that for a finite system, the wave function
φ(x) corresponds to the eigenstate with the largest eigenvalue of the single-particle
density matrix, $Φ̂† (x# )Φ̂(x)#.
The equations that determine the dynamic motion of the condensate wave function are found by calculating the Heisenberg equation of motion for the Φ̂ operatorih̄∂ Φ̂(x)∂t= [Φ̂, Ĥ]=!

−(2.21)"h̄2 24πah̄2 †∇ + VL (x) + VE (x) +Φ̂ (x)Φ̂(x) Φ̂(x) .
2mm(2.22)

To zeroth order, one can replace Φ̂(x) → φ(x) to give
∂φ(x)
=
ih̄
∂t

!

h̄2 2
4πah̄2−∇ + VL (x) + VE (x) +|φ(x)|2 φ(x) ,2mm"(2.23)

which is known as the Gross-Pitaevskii (GP) equation [65, 127], or the nonlinear Schrödinger(NLS) equation. The GP equation is valid when the interparticle spacing is much greaterthan the s-wave scattering length, the number of condensed atoms is much greater thanone, and the temperature is low so that the number of condensed atoms is much greaterthan the number of atoms in any of the excited states. However, the theory is only

21
strictly valid for δΦ̂(x) = 0 which corresponds to zero temperature. Although the system is dilute, ρ|a|3 " 1, interactions can have large effects on the system when theinteraction energy is comparable in size to the kinetic energy of the atoms.
The Gross-Pitaevskii equation will be used to describe the nonlinear band structure of atoms in a periodic potential.

2.4

Bose-Hubbard Model
Another method used to describe the motion of atoms in a periodic potential isgiven by the Bose-Hubbard Hamiltonian. This model focuses on the relative influenceof the interaction energy of two-body collisions and the kinetic energy associated withmotion between lattice sites. The use of the Bose-Hubbard Hamiltonian is very common in condensed matter and solid state systems in which the motion of electrons is
examined. The presence of a quantum phase transition between Mott insulator andsuperfluid states has been shown theoretically and experimentally [52, 64]. The detailsof this quantum phase transition will be described in a later chapter.
As with the examination of the GP equation, the Bose-Hubbard Hamiltonian isderived from the second-quantized many-body Hamiltonian,
Ĥ =

!

"

#

h̄2 2
dxΦ̂† (x) −∇ + VL (x) Φ̂(x)
2m

+

#

+

1
2

(2.24)dxΦ̂† (x)(VE (x) − µ)Φ̂(x)#dxdx# Φ̂† (x# )Φ̂† (x)VI (x − x# )Φ̂(x)Φ̂(x# ) .

The difference lies in the way in the approximations used to decompose the bosonicannihilation operator. The operator will be decomposed into a set of site localized
Wannier states that are related to the Fourier transform of the Bloch states. This
decomposition is valid as long as the interactions are smaller than the band gap and thelattice is strong enough to prevent next-to-nearest neighbor hopping between sites.

22
The creation operator can be expanded in any complete set of basis states. Although the Bloch states may seem the most direct basis states to use, there is a set of
localized states that often proves more convenient. The Wannier states are defined as
1 $ −iqxi (n)eψq (x) ,N q

wn (x − xi ) =(2.25)

where N is the total number of lattice sites, the sum is over all states in the first Brillouin(n)zone, xi is positioned in the center of site i, and ψq (x) is the Bloch state solution of asingle particle in the lattice with quasi-momentum q in the nth band. With the use of
the Wannier states, the annihilation operator can be given by,Φ̂(x) =$
in(n)

where âi(n)

âi wn (x − xi ) ,(2.26)

is the annihilation operator of an atom at site i in the nth band. If thechemical potential is less than the energy required for a single particle excitation to
the second band, only the states from the lowest band need to be considered. The
annihilation operator then simplifies toΦ̂(x) =$
iâi w0 (x − xi ) ,(2.27)

which will be used in describing the Bose-Hubbard Hamiltonian. These are not the only
way of defining Wannier states. For instance, a phase gradient can be applied to each
Wannier state to better describe the effects of rotation [13].
The expansion of annihilation operator in terms of Wannier states can be inserted
into the Hamiltonian to find
H = −J$
1 $ † †ân ân ân ân +(Vn − µ)â†n ân ,â†n âm + U
2
n
n
$n,m%

J

= −

#

U

=

$

(2.28)

where

Vn =

"

h̄2 ∂ 2dxw0∗ (x − xn ) −+ V (x) w0 (x − xn+1 ) ,2m ∂x24πah̄2m#!

#

dx|w0 (x)|4 ,dx|w0 (x)|2 V (x) ,(2.29)(2.30)(2.31)

23
and $n, m# indicates a sum over nearest neighbors. It should be noted that only terms
in which the overlap integral of Wannier states that are at most one lattice site inseparation are considered since the remaining terms are much smaller. The higher
order terms are often approximately two orders of magnitude smaller than the nearest
neighbor terms. This is called the tight-binding approximation and assumes that theatoms are primarily described by a set of atoms in separate lattice sites with a small
correction due to hopping between lattice sites. The tight-binding regime occurs when
the lattice depth becomes approximately V0 > 5ER , where ER ≡ h̄2 k2 /2m is the atomic
recoil energy.
The interaction strength can be characterized by the parameter U while the kinetic energy can be characterized by the hopping parameter J. In an optical lattice
it is possible to alter the ratio J/U by adjusting the laser intensity. The hopping rate
decreases exponentially as the intensity is increased while the interaction strength increases according to a power law. The alteration of the ratio J/U changes the relative
effects of the kinetic energy and the interaction energy. Therefore, it is possible tohave a dilute system, ρa3s " 1, with large effects due to interactions, J/U ≈ 1. The
assumptions of the Bose-Hubbard model are valid as long as the interaction energy U N
is much smaller than the gap between the first and second bands.
A mean field approximation can be applied to the Bose-Hubbard model in which
the bosonic operators are replaced with a c-number, as was done in the GP equation.
If interactions are neglected, U = 0, the Hamiltonian can be used in what is called thediscrete Schrödinger equation or tight-binding Schrödinger equation. If interactions, arenot excluded, the discrete nonlinear Schrödinger equation can be used.
The Bose-Hubbard Hamiltonian with a site dependent external potential, will be
the basis for the model of an atomtronic system.

Chapter 3

Nonlinear Band Structure

The nonlinear Schrödinger equation (NLS), as introduced in the previous chapter,models many kinds of wave phenomena. The NLS appears in diverse fields such as
nonlinear optics [167], gravity waves on deep water [73], magneto-static spin waves[151], solitons in liquid crystals [36], and magneto-sonic solitons in the atmosphericmagneto-pause boundary layer [160]. It also describes the dynamics of the mean field ofa weakly interacting atomic or molecular Bose-Einstein condensate (BEC) [6, 39, 19, 18],
where it is known as the Gross-Pitaevskii equation [65, 127].
Periodic potentials are ubiquitous in physics, appearing in electron transport inmetals [4], Josephson junction arrays [171], nonlinear photonic crystals and waveguidearrays [35], and Bose-Einstein condensates (BEC’s) [63]. With the realization of BEC’sof alkali atoms in a sinusoidal optical lattice, there has been an explosion in studies of
BEC’s in these potentials, both experimentally and theoretically [64, 49, 46, 48, 45, 50].
BEC’s in optical lattices have been used to study phase coherence of atom lasers [5, 70]
and matter-wave diffraction [121]. Therefore, in the context of the BEC, the study of
periodic potentials provides an excellent connection between condensed matter physics
and atomic physics. In contrast to other physical contexts, the lattice geometry andstrength, as well as the interatomic interactions [136, 81], are all tunable parameters for
the BEC.
Experimentally, in order to create a BEC in a lattice, alkali atoms are first cooled

25
to a quantum degenerate regime by laser cooling and evaporation in harmonic electromagnetic traps [6, 39]. A sinusoidal optical lattice is then made by the interference of
two counter-propagating laser beams, which creates an effective sinusoidal potential proportional to the intensity of the beams [64]. The potential is due to the AC Stark shift
induced by the dipole interaction with the electromagnetic field on the atoms’ center ofmass [110]. For large detuning of the optical field from the atomic transition, dissipative
processes, such as spontaneous emission, can be minimized and the potential becomes
conservative. Nonzero quasi-momentum can be examined by slightly detuning the twolasers by a frequency δν [50]. The resulting interference pattern is then a traveling wave
moving at the velocity v = (λ/2)δν, where λ is the wavelength of the first beam. This
produces a system with quasi-momentum q = mv/h̄, where m is the atomic mass. Notethat this is not necessarily the velocity and momentum of the atoms in the lattice since
dragging may occur. The actual velocity is derived from the dispersion relation of theband structure. After a given evolution time, the traps are switched off. The BEC isallowed to expand and the density is then imaged.
With the experimental demonstration of Feshbach resonances in BEC’s of diluteatomic gases [136, 81], it is possible to alter the s-wave scattering length and, hence,
the nonlinearity of the NLS. Near a Feshbach resonance, the scattering length becomesa function of a uniform background magnetic field. By altering the magnetic field, theeffects of the nonlinearity can be experimentally controlled. It is therefore important
to be able to characterize the complete set of solutions as a function of the interaction
strength.
The complete set of stationary solutions to the NLS with a constant potential
on the infinite line were discovered by Zakharov and Shabat [177, 178]. The stationary solutions of the NLS under periodic and box boundary conditions has also beensolved analytically [26, 27], as well as the finite well [30]. The parabolic potential has
been solved numerically [93]. The potential step has been examined theoretically and26experimentally for the linear Schrödinger equation with a constant [172, 118] and oscillating [51] step. Klein examined similar problems experimentally by deflecting neutron
beams with a vibrating crystal [96]. In addition, symmetric steady-state solutions with apoint-like impurity potential have been studied by Hakim [71] and Taras-Semchuk [168]in certain limiting cases. Point-like impurity potentials have been studied extensively,
such as with helium impurities in a BEC [25, 72], BEC formation initiated by point
like impurities [32] and impurity scattering in a BEC of sodium [34]. The superfluid
transmission of matter waves across various potentials has been studied [99, 122]. In
addition, bound solutions to the NLS are the one-dimensional analog to the pinned vortex solutions which occur when a discontinuity is present in a two dimensional system,
such as in a two dimensional high-Tc superconducting system.
The case of a BEC trapped in a sinusoidal potential has been studied theoreticallyin great detail by a number of researchers [101, 174, 175, 23, 24, 28, 164, 165, 43, 117,
173, 40, 107, 108, 105, 134, 169, 120, 133, 129]. In addition, a subset of the solutions
to the Kronig-Penney potential has been found [170, 55, 168, 2, 102]. Recently, forexample, Li and Smerzi [102] investigated generalized Bloch states for constant phaseand zero current. In contrast, we present the full set of Bloch wave stationary solutions
for both repulsive and attractive BEC’s. An important point is that the sinusoidal
optical lattice potential is composed of a single Fourier component. If more counterpropagating laser beams of different frequencies are added as is experimentally possible,
more Fourier components are introduced, and the potential becomes a lattice of well
separated peaks. In the limit that the width of the these peaks becomes much smaller
than the healing length of the BEC, the potential effectively becomes a Kronig-Penneylattice. However, we show that the Kronig-Penney potential serves as a good modeleven for experiments with a single Fourier component.
Here, we consider the steady state response of the mean field of a BEC to a
constant potential, a potential step, a delta function potential, and a Kronig-Penney

27
potential as modeled by the NLS. These potentials may be easily produced in present
experiments on the BEC. The step potential may be realized by a detuned laser beamshined over a razor edge to make a sharp barrier, where the diffraction-limited fall-offof the laser intensity is smaller than the healing length of the condensate, so that thepotential is effectively a step function. The delta function potential models the responseof the condensate to an impurity of a length scale smaller than the healing length, which
could be realized by a tightly focused laser beam, by another spin state of the sameatom, or by any other object, as for instance another alkali atom, confined in an opticaltrap.
The delta-function lattice is an approximation to the easily created optical latticepotential. Specifically, we consider the steady state response of the mean field of a
BEC to a Kronig-Penney potential using the Bloch ansatz. We first obtain the fullset of stationary solutions for a single delta function analytically [152, 14, 10, 130].
Using these results, we are able to rigorously describe the band structure. We find thatabove a critical nonlinearity, swallowtails, or loops, form in the bands [173, 107, 153].
Stability properties are studied numerically by time evolution of perturbed stationary
state solutions to the NLS. It is found that stable, as well as unstable, regimes exist
for both repulsive and attractive BEC’s. Finally, the appearance of swallowtails isconnected to period doubled states by examining the band structure of a two-color
lattice [154].
We show that the mean field Bloch states of a BEC in a Kronig-Penney potential,i.e., a lattice of delta functions, exhibits the same band structure and stability properties
as the experimental case of a sinusoidal potential. Unlike in the case of the sinusoidal
potential, Bloch state solutions to the Kronig-Penney potential can be described bystraightforward analytic expressions. The Kronig-Penney potential, therefore, has distinct advantages as a model of BEC’s trapped in periodic potentials.
In this chapter, we first examine a condensate subject to a constant potential,

28
essentially the free-particle states. These results can be used to analyze the effects on acondensate due to a step potential or delta function potential, both of which are single
boundary potentials. A lattice of delta functions, a Kronig-Penney potential, models
an optical lattice and is then used to create a full nonlinear band structure analytically.
Finally, the appearance of swallowtails, or loops in the band structure, are connected
to period-doubled states by examining a two-color lattice.

3.1

Constant Potential
We examine the BEC in the quasi-one-dimensional regime and obtain the fullset of stationary states in closed analytic form. When the transverse dimensions of
the BEC are on the order of its healing length, and its longitudinal dimension is much
longer than its transverse ones, the 1D limit of the 3D NLS is appropriate to describe
the system [29]. This corresponds to tight confinement with an optical lattice in two
directions creating a 2D set of one dimensional tubes. A third weaker lattice is createdin the free dimension. The 1D NLS, with an external potential, V (x), may be written
1i∂t Ψ = − ∂xx Ψ + g|Ψ|2 Ψ + V (x)Ψ ,
2

(3.1)

where a harmonic oscillator confinement in the transverse directions with frequency ω
has been assumed [119] for atoms of mass M, the length has been rescaled according
to units of the oscillator length, lho = (h̄/M ω)1/2 , and energy rescaled according tounits of the oscillator energy, h̄ω. The renormalized 1D coupling, g ≡ 2as , where
as is the s-wave scattering length, characterizes the short-ranged pairwise interactions
between atoms. The wave function or order parameter, Ψ(x, t), can be written asΨ(x, t) =%ρ(x, t) exp[iφ(x, t)], where ρ(x, t) is the longitudinal line density and thelongitudinal superfluid velocity is given by v(x, t) = ∂φ(x, t)/∂x. Both attractive andrepulsive atomic interactions, i.e., g < 0 and g > 0, shall be considered.
In the case where the harmonic oscillator length approaches the s-wave scattering

29
length, lho , as , the 1D NLS no longer models the system and a one-dimensional field
theory with the appropriate effective coupling constant must be considered instead [119].
Since as is on the order of angstroms for typical BEC’s, this regime is not relevant to
the present study. Thus, it should be noted that this study does not examine theTonks-Girardeau regime [57], where quantum fluctuations become important and theGross-Pitaevskii equation no longer models the system.
The ability to analytically describe the full nonlinear band structure of a condensate derives from the stationary solutions of a free condensate [177, 178]. Stationary
states are given by the form
Ψ(x, t) = ψ(x)e−iµt ,(3.2)

where µ is the eigenvalue of the stationary solution and converts the partial differential
GP equation, Eq. (3.1), into the ordinary differential equation1µψ = − ∂xx ψ + g|ψ|2 ψ + V (x)ψ .
2

(3.3)

The wave function can be further broken up to describe the density and phase asψ(x, t) = R(x)eiφ(x) ,(3.4)

where |R(x)|2 is the condensate density at x and φ(x) is the condensate phase at x.
There are several excellent references which consider stationary solutions in a constantpotential [80, 95, 29, 79, 176]. Assuming a constant external potential V , the separationof Eq. (3.1) into real and imaginary parts gives
R ∂xx φ + 2(∂x R) (∂x φ) = 0 ,1
1∂xx R + (∂x φ)2 R + g R3 + V R = µ R .
2
2

(3.5)
(3.6)

Equation (3.5) can be integrated to give∂x φ =

α
,
ρ

(3.7)

30
where α is an undetermined constant of integration, and ρ(x) = |R(x)|2 is the singleparticle density. Substituting Eq. (3.7) into Eq. (3.6) and integrating yields&dρ= 2 g ρ3 − 2(µ − V )ρ2 + C ρ − α2 .
dx

(3.8)

Integrating a second time gives#1%gρ3 − 2(µ − V )ρ2 + Cρ − α2dρ = 2x + x0 ,(3.9)

where C and x0 are undetermined constants of integration. It is shown in App. B
that the only solutions to this integral equation are given by the Jacobian elliptic functions [17, 1]. The Jacobi elliptic functions can be considered generalizations of thetrigonometric and hyperbolic functions and interpolate between the two. In App. A
these special functions are briefly reviewed. The most general form of the solution isthen given byρ = A sn2 (b x + x0 , k) + B ,(3.10)

where sn is one of the Jacobian elliptic functions, A is a density prefactor, b is atranslational scaling, x0 is a translational offset, k is the elliptic parameter, and B isa vertical density offset. The period of the density is given by 2 K(k)/b, where K(k)is the complete elliptic integral of the first kind. Since the density is positive definite,
the variables are restricted such that B ≥ 0 and B + A ≥ 0. In Sec. 3.5 the relation
between these variables and the mean number density, energy density, and momentum
density are calculated and discussed for the nonlinear wave of Eq. (3.10).
It should be noted that since many of the solutions are unbounded, the normof the wave function remains unconstrained. This is in contrast to bound or localized
solutions when the normalization,# +∞−∞(ρ − ρ) dx = 1 ,(3.11)

may be used [95], where ρ is the average density. This allows for a normalization ofsystems with a localized disturbance to a constant background field. It is possible to31
use a nonlinear scaling of the position and density, x → ax and ρ → a−2 ρ, to scale themaximum density on one side of a boundary to unity. However, for considerations of
single boundary systems the wave function remains unnormalized. Normalization for aperiodic potential will be such that the norm of each lattice site is set to unity.
One may then determine the variables µ, α, and k, in terms of A, B, b, and g by
substituting Eq. (3.10) into Eq. (3.8) and equating powers of the elliptic sn function to
give
µ =

1 2
(b + (A + 3B)g) + V ,2

α2 = B(A + B)(b2 + Bg) ,
k2 =

A
g.
b2

(3.12)(3.13)(3.14)

This leaves the eigenvalue, µ, the constant of integration of the phase, α, and theelliptic parameter, k, determined up to A, B, b, and the interaction strength, g. Notethat the fact that α enters into the equations only as α2 implies that all nontrivial phasesolutions, i.e., those for which α )= 0, are doubly degenerate, as ±α lead to the samevalue of the eigenvalue, µ, without otherwise changing the form of the density or phase.
This corresponds to the energetic degeneracy of left moving and right moving systems.
We shall use the term trivial phase to refer to solutions for which the phase is spatiallyconstant.
In the following two sections these results are applied to piecewise constant potentials. In particular, the potential step and the delta function potential are examined.

3.2

Potential Step
In this section the complete set of solutions to the NLS with an external step

potential of height V0 ,V (x) = V0 θ(x) ,(3.15)

32
are presented analytically, given the solution parameters on the lower side of the step.
In the following two subsections, the general solution to the NLS with a potential stepand some particular examples are discussed.

3.2.1

General Solution

Applying the two boundary conditions of continuity of the wave function andcontinuity of the derivative of the wave function gives the following five conditions,ρ(0+ ) = ρ(0− ) ,(3.16)∂x ρ(0+ ) = ∂x ρ(0− ) ,(3.17)

φ(0+ ) = φ(0− ) + 2πn ,(3.18)

α(0+ ) = α(0− ) ,(3.19)µ(0+ ) = µ(0− ) ,(3.20)

where n is an integer. The first two conditions represent the continuity of the density,
Eq. (3.16), and the derivative of the density, Eq. (3.17). The next two conditions represent the continuity of the phase, Eq. (3.18), and the derivative of the phase, Eq. (3.19).
Equation (3.20) demands that the eigenvalue be the same on either side of the boundary.
Note that in Eq. (3.12) the eigenvalue on the right hand side is offset by V0 as comparedto the left hand side. Since n only enters into the phase and does not effect such quantities as the eigenvalue and density, only the n = 0 state is considered and, consequently,
the phase, φ, is continuous across the boundary. However, it is important to note thatall solutions are of denumerably infinite degeneracy, according to Eq. (3.18). Since φ isgiven byφ(x) =# x0αdx + const. ,ρ(x)(3.21)

continuity in the phase is easily achieved by setting the constant phase shift equal oneither side of the boundary and, therefore, Eq. (3.18) is satisfied.

33
In the following derivation, it is assumed that the wave function parameters on
the left side of the step are known completely. Therefore the density prefactor, AL , thevertical density offset, BL , the translational scaling, bL , and the horizontal offset, x0L
are all known, where the L subscript refers to variables on the left side; an R subscript
will refer to variables on the right side. In addition, the experimental parameters of
the interaction strength, g, and the potential step height, V0 , are both known. Fromthe variables on the left and Eq. (3.10), the density at the boundary, ρL ≡ ρ(0− ), andits derivative at the boundary, ∂x ρL ≡ ∂x ρ(0− ), can be determined. The eigenvalue, µ,and the phase constant, α, can be determined from Eqs. (3.12) and (3.13).
From Eqs. (3.10) and (3.16) the square of the Jacobian elliptic sn function canbe solved for,
sn

2

!

x0R ,

'

AR g
b2R

"

=

ρL − BR.
AR(3.22)

From Eqs. (3.12) and (3.20), the horizontal scaling, bR , is given by
b2R = 2(µ − V0 ) − (AR + 3BR )g .

(3.23)

These variables are substituted into Eqs. (3.19) and (3.17), using Eq. (3.13), to give
α2L = BR (AR + BR )× [2(µ − V0 ) − (AR + 2BR )g] ,(3.24)(∂x ρL )2 = −4(BR − ρL )(AR + BR − ρL )× [2(µ − V0 ) − (AR + 2BR + ρL )g] .

(3.25)

Equations (3.24) and (3.25) are quadratic in AR and cubic in BR and can be solved
analytically for AR and BR to give six solutions. The remaining variables on the rightside can then be found by substituting the values of AR and BR into Eq. (3.23), to findbR , and by taking the inverse Jacobi sn function of Eq. (3.22) to givex0R = sn−1!

ρ L − BR
,
AR

'

AR g
b2R

"

.

(3.26)

34
The full solution is then completely known. It is therefore possible to completely describe
the system analytically knowing only the parameters on one side of the step. This not
only introduces computational ease in evaluating solutions, it also provides all possible
solutions, most of which cannot be determined using purely numerical methods. In thefollowing section, specific examples of a BEC in the presence of a step potential will be
examined.

3.2.2

Particular Examples

The solutions to the potential step problem can be divided into two categories.
The eigenvalue, µ, can be large enough that particles are free to move across the boundary. It is then possible to examine the nonlinear analog to a transmitted wave for the
linear Schrödinger equation. However, if the eigenvalue is too small, then the wave
function must decay under the step.
When µ is larger than the effective potential,
Vef f (x) = V0 + gρ(x) ,(3.27)

the wave can be transmitted across the boundary. Note that in the case of an attractive
interaction, g < 0, the eigenvalue can be less than the step height, V0 , and for a
repulsive interaction, g > 0, the eigenvalue must be strictly greater than the step height.
In Figs. 3.1(a) and (b), the density and phase of a nonlinear state with a repulsive
interaction strength is shown. A step of height V0 = 1, positioned at x = 0, and acondensate with an interaction strength of g = 0.2 and eigenvalue of µ = 2.404 wereused. The increased interaction strength and nonlinearity has caused the peaks of the
wave function to become much broader than in the linear case, which would be simply
sine waves. Figure 3.1(c) and (d) present the density and phase of a similar solution
but with an attractive interaction strength. This potential is again given by a stepwith height of V0 = 1, positioned at x = 0. An interaction strength of g = −0.2 and

(a)

0
-10

0

10

Position (units of lho)

(b)

0
-5

20 -10
-10
1.5

(c)

Phase / 2 π

Density (units of lho-1)

3.5

40

35

7.5
Phase / 2 π

Density (units of lho-1)

7

20
0
-10 0

30

Position (units of lho)

60

0

0

10

Position (units of lho)

20

(d)

-2
-5
-10 0

30

Position (units of lho)

60

Figure 3.1: Stationary solutions to the NLS with a potential step of form V (x) = V0 θ(x).
These solutions, which are the nonlinear analogs of the continuum of linear stationary
scattering states, exhibit a large deviation from the traditional linear solutions. Shownare particular examples of (a) the density and (b) the phase for a repulsive interaction
strength and (c) the density and (d) the phase for attractive interaction strength.

eigenvalue of µ = 0.98 was used. In this case, the peaks have instead narrowed due tothe attractive interaction.
When the eigenvalue is smaller than the effective potential, µ < Vef f (x), the wave
function must decay under the step. In Fig. 3.2(a) the density of a nonlinear solution
with a repulsive interaction strength that decays as it crosses the boundary of the stepis shown. A step with height of V0 = 1, positioned at x = 0, and a condensate withan interaction strength of g = 0.12 and eigenvalue of µ = 0.5 were used. Figure 3.2(b)shows a similar solution but with an attractive interaction strength. This potential isagain given by a step with height of V0 = 1, positioned at x = 0. An interaction strengthof g = −10 and eigenvalue of µ = −49 were used. For both wave functions, the phaseis necessarily trivial, since all wave function that approach zero at infinity must have
A = −B, and hence from Eqs. (3.7) and (3.13), the phase is constant. In order to stress

2.5

-10

0

Position(units of lho)

10Density (units of lho-1)

2

(b)

5
0-2

-1

0 0.5

Position (units of lho)

2

(c)

1

0
-10

Density (units of lho-1)

(a)

0
-25

36

10

0

10

Position (units of lho)

20

Density (units of lho-1)

Density (units of lho-1)

5

(d)

1

0
-10

0

10

Position (units of lho)

20

Figure 3.2: Stationary solutions to the NLS with a potential step for which the wave
function decay, which are the nonlinear analogs to evanescent waves. Shown are particular examples of the densities of nonlinear waves with (a) repulsive interaction strength,(b) attractive interaction strength, and (c) and (d) with an attractive interaction
strength that decay on the lower side of the step.

the importance of the effective potential, and not just the step potential, Figs. 3.2(c)and (d) show the density of two nonlinear wave function with attractive interactions,g = −1, that decay occurs on the lower side of the potential. Both wave function havean eigenvalue of µ = −0.5 and a potential of height of V0 = 0.01 and V0 = 1 were used
for Fig. (3.2(c) and Fig. 3.2(d), respectively.
Thus the NLS with a potential step has solutions to the step potential that provide
a connection between the linear solutions and a wide variety of exotic nonlinear wave
function, as shall be discussed in Sec. 3.4.

3.3

Point-like Impurity
In this section, the case of a potential of form
V (x) = V0 δ(x)(3.28)

37
is considered. Such a potential models an impurity which deforms the constant background potential on a length scale much less than that of the healing length. Given
the state of the system on the negative x, or left side of the impurity, the variableson the right side of the impurity are determined. A negative value of V0 represents
an attractive impurity, such as due to defects in hydrogen-bonded chains [94, 92], anda positive value of V0 represents a repulsive impurity, such as with helium atoms in a
BEC [25, 72].

3.3.1

General Solution

The boundary conditions for an impurity are similar to those for the potential
step, except that the derivative of the wave function experiences a discontinuity at theboundary. Therefore, Eqs. (3.16), (3.20), (3.19) and (3.18) must still be satisfied, alongwith∂x ρ(0+ ) − ∂x ρ(0− ) = −4ρ(0)V0 .

(3.29)

It is again assumed that all variables on the left side of the impurity are known aswell as the experimental parameters of interaction strength, g, and impurity strength,
V0 . Using a treatment that is exactly analogous to that for the step function, all of theparameters on the right side, given those on the left, are analytically determined. The
only difference is that in Eqs. (3.23), (3.24), (3.25), and (3.26), the quantity ∂x ρL must
be replaced with (∂x ρL − 4V0 ρL ) and (µ − V0 ) must be replaced with µ. It is therefore
possible to completely describe the system analytically knowing only the parameters on
one side of the impurity. In the following section, examples of wave functions subject
to an impurity are examined.

3.3.2

Particular Examples

For the delta function potential, both symmetric and nonsymmetric wave function
are possible. Of particular interest are the symmetric wave function in the k = 1

38limit of the Jacobian elliptic functions. In this case all solutions become hyperbolic
trigonometric functions with a localized change in the density around the impurity andno oscillations at ±∞. Solutions of this type we term localized [30] and are of four
different forms:%ρ = α2 + (1 − α2 ) tanh2 ( 1 − α2 |x| + x0 ) ,%(3.30)

ρ = α2 + (1 − α2 ) coth2 ( 1 − α2 |x| + x0 ) ,(3.31)

ρ = b2 sech2 (b|x| + x0 ) ,(3.32)

ρ = b2 csch2 (b|x| + x0 ) ,(3.33)

where the translational offset, x0 , is determined by the impurity strength, V0 , and thedensity has been normalized according to Eq. (3.11). Equations (3.30) and (3.31) are
valid for repulsive interactions, while Eqs. (3.32) and (3.33) are valid for attractive
interactions.
Figure 3.3 shows the possible wave function for repulsive interactions. An exampleof the solution described by Eq. (3.30) is plotted in Fig. 3.3(a) with a repulsive impurityof strength V0 = 0.5. This may be interpreted as a single dark soliton bound by an
impurity. This is similar to the solution found by Hakim [71] for a soliton that ismoving with an impurity. In addition, a bound state of two dark solitons, Fig. 3.3(c),can be created when the strength of the impurity is attractive and exactly balances therepulsion between the two dark solitons. Figure 3.3(e) shows the hyperbolic cotangent
function solution with an impurity strength of V0 = −0.5; this may be interpreted as adeformation of the ground state constant solution to the NLS with a constant potential.
In all plots an interaction strength of g = 1 and phase constant of α = 0.5 were used.
It should be noted that there exists a bound state of a repulsive condensate withan attractive impurity. This solution is given by,
ρ=

−V0
b2
csch2 bx + coth−1
g
b
(

(

))

,

(3.34)

39
where the interaction strength, g has specifically been included and b must be determinedsuch that the density is normalized to unity. These requirements place a limit on how
repulsive the interaction may become and is given by,gmax = −4V0 ,(3.35)

where gmax is the most repulsive interaction the condensate may have. If the interaction
is increased past this point, the condensate will spill away from the impurity and will
no longer be bound.
The set of symmetric localized solutions for the case of attractive interactions do
not allow for nontrivial phases, in contrast to the case of repulsive interactions. For
attractive interactions, g < 0, the hyperbolic secant function solution, Eq. (3.32), isvalid for both V0 > 0, Fig. 3.4(a), and V0 < 0, Fig. 3.4(b), where potential strengths of
V0 = 0.9 and V0 = −0.9 were used, respectively. These solutions may be interpreted as asingle bright soliton, which is the ground state solution to the 1D-NLS, deformed by an
impurity. The hyperbolic cosecant function solution, Eq. (3.33), is only valid if V0 < 0and is similar in form to the hyperbolic secant solution of Fig. 3.4(b). The two solutions
types are degenerate for V0 < 0, with an eigenvalue of µ = −b2 /2. In Figs. 3.4(a) and(b) an interaction strength of g = −1 and translational scaling of b = 1 were used.
Nonsymmetric wave function which oscillate at infinity are also possible and comein two forms, oscillations on one side of the delta function and oscillations on both sides.
Figure 3.5 shows two possible nonsymmetric wave function subject to a delta function
positioned at x = 0 that oscillate on one side of the potential. The density of a wave
function with a repulsive interaction strength of g = 0.21 and eigenvalue of µ = 2.4,
distorted by a delta function, V0 = 2, is shown in Fig. 3.5(a), where the left side
reproduces the hyperbolic tangent function of Eq. (3.30). In Fig. 3.5(b), an attractively
interacting, g = −50, wave function with an eigenvalue of µ = −50, distorted by a delta
function, V0 = 10, is shown that appears similar to the evanescent wave function of

(a)

1

0

Position (units of lho)

0

5

0

5

0

5

Position (units of lho)

1.5
(c)

1
0-5

0

(d)

Position (lho)05 -1.5-52

Position (units of lho)1(e)Phase / 2 π

Density (units of lho-1)

0

5 -1-5

2

1

0-5

(b)

Phase / 2 π

Density (units of lho-1)

0-5

40

1
Phase / 2 π

Density (units of lho-1)

20

Position (units of lho)(f)0

5 -1-5

Position (units of lho)

Figure 3.3: Localized, symmetric solutions to the NLS with repulsive interaction
strength in the presence of an impurity, V (x) = V0 δ(x). Shown are particular examples of (a) the density and (b) the phase of a dark soliton bound by a repulsive
impurity, (c) the density and (d) the phase of a pair of dark solitons bound by an attractive impurity, and (e) the density and (f) the phase of a supercurrent deformed by
an attractive impurity. Note that (a) and (c) may also be interpreted as deformationsof a supercurrent.

Fig. 3.2 that decay beneath a step. Note that Fig. 3.5(a) has a nontrivial phase while
the phase of Fig. 3.5(b) is trivial.
Figure 3.6 shows two possible nonsymmetric wave function subject to a delta
function positioned at x = 0 with strength V0 = 2 that oscillate on both sides ofthe potential. A repulsive interaction strength produces the characteristic wideningof the pulse peaks, Fig. 3.6(a). The corresponding phase is given in Fig. 3.6(b). An
attractive interaction strength creates a narrowing of the pulse peaks, Fig. 3.6(c). The

(a)

41

1
(b)

0.5

1

0-5

Density (units of lho-1)

Density (units of lho-1)

20

Position (units of lho)

5

0-5

0

Position (units of lho)

5

Figure 3.4: Localized, symmetric solutions to the NLS with attractive interaction
strength in the presence of an impurity. Shown are particular examples of (a) thedensity of a bright soliton, which is the ground state solution to the NLS, deformed bya repulsive impurity and (b) the density of a bright soliton deformed by an attractive
impurity.

corresponding phase is given in Fig. 3.6(d). The condensates in Figs. 3.6(a) and (c) arecharacterized by eigenvalues of µ = 2.4 and µ = −1.3, respectively.

3.4

Linear Limits of Transmitted and Reflected Waves
In this section it is shown how the solutions to the NLS connect to the solutionsof the linear Schrödinger equation. There are two distinct types of possible waves. If
the wave has enough energy it is possible to make the wave propagate through space.
However, if it does not have enough energy, the wave can carry no current and is referred
to as an evanescent wave. These waves must decay. In the next two sections, these two
types of waves are discussed concerning their role as the linear limit of the nonlinear
solutions.

3.4.1

Transmitted Waves

In a linear system, if the energy of the system is greater than the potential energy,
then the wave can be transmitted through space as sine waves. The nonlinear Jacobian
elliptic sn function waves are the nonlinear analogue of the sine waves of a linear system.

Density (units of lho-1)

Density (units of lho-1)

(a)

10

(b)

20
10

5
0-5

42

30

15

0

5

Position (units of lho)

0
-0.5

0

Position (units of lho)

0.5

Figure 3.5: Solutions to the NLS with an impurity that oscillate on one side of the potential. Shown are particular examples of (a) the density of a nonlinear wave with repulsive
interaction strength and (b) the density of a nonlinear wave with attractive interaction
strength. These waves have no analog with the solutions of the linear Schrödinger
equation.

In this section, it is shown how the nonlinear solutions, in the linear limit, recreate thelinear solutions of propagating waves.
For the linear case, the usual representation of the system by incident, transmitted, and reflected waves is given by,*

ΨL (x, t) =ΨR (x, t) = T eikR x e−iµt ,where kL =%2(µ − VL ), kR =+eikL x + Re−ikL x e−iµt ,(3.36)(3.37)%

2(µ − VR ), µ is the eigenvalue of the Schrödinger

equation, R is the reflection coefficient, T is the transmission coefficient, and the incident
wave is assumed to be coming in from the left. The potentials, VL and VR , and thescattering coefficients, R and T , are determined by the type of boundary. For the caseof a potential step, the potentials are given by VL = 0 and VR = V0 . For a delta
function potential, the potentials are given by VL = VR = 0. The wave function can bealternately described by an amplitude and phase as follows:ρL = (1 + r)2 − 4 r sin2 ( 2 µ x − s/2) ,√(r − 1)(r + 1) 2 µ√,∂x φL =
(1 + r)2 − 4 r sin2 ( 2 µ x − s/2)%

(3.38)(3.39)

Phase / 2 π

Density (units of lho-1)

(a)

10
5
0

5

Position (units of lho)

3

0

0

5

10

15

0

4

Position (units of lho)

6
(c)

2
1

0
-10

(b)

15 -3-5

10

Phase / 2 π

Density (units of lho-1)

0-5

43

3

15

-5

(d)

3
0

4 -3
-10

0

Position (units of lho)

-5

Position (units of lho)

Figure 3.6: Nonsymmetric solutions to the NLS with an impurity. Shown are particular
examples of (a) the density and (b) the phase of a nonlinear wave with repulsive interaction strength and (c) the density and (d) the phase of a nonlinear wave with attractive
interaction strength. These solutions are the nonlinear analogs to the continuum of
linear stationary scattering states.

whereΨL (x, t) =&ρL (x) exp[iφL (x)] e−iµt ,R = reis ,(3.40)(3.41)

with r and s real. The nonlinear solution is then connected to the linear solution by
x0 = −s/2 ,(3.42)

b2 = 2µ ,(3.43)

A = −4r ,(3.44)

B = (r + 1)2 .

(3.45)The transmitted wave can easily be given byρR = t 2 ,√2µ∂x φR =,
t2

(3.46)(3.47)

44
whereΨR (x, t) =T&ρR (x) exp[iφR (x)] e−iµt ,= teiw ,(3.48)(3.49)

with t and w real. The nonlinear solution is then connected to the linear solution by
x0 = 0 ,(3.50)

b2 = 2µ ,(3.51)

A = 0,(3.52)

B = t2 .

(3.53)

When the eigenvalue is greater than the effective potential, the nonlinear solutions are
therefore adiabatically connected to the linear transmitted wave solutions. The next
section connects the decaying evanescent waves with nonlinear solutions.

3.4.2

Evanescent Waves

If the eigenvalue, µ, is smaller than the effective potential then the wave function
must decay. In the linear case, the decay is precisely exponential. The Jacobian elliptic
sn function solution, Eq. (3.10), of the NLS can provide the appropriate exponential
decay. In the limits g → 0 and B, x0 → +∞ under the constraints A=−B and A g/b2 =1, Eq. (3.10) gives
lim

B,x0 →+∞

B sech2 (bx + x0 ) = ρ0 e−bx ,(3.54)

where ρ0 is the density at the boundary and is given by,
ρ0 =

4B
.
e2x0

(3.55)

When the limit that g approach zero is not enforced,
ρ(x) ∝

1(N+ ebx + N− e−bx )2,(3.56)

45
where N+ and N− are constants related to the magnitude and sign of the interaction
strength. The decay of the density is, therefore, not strictly an exponential decay. It isinteresting to note that all decaying solutions, whether linear or nonlinear, must possessa trivial phase. This is due to the restriction that A=−B, and, hence from Eqs. (3.13)
and (3.7), the phase must vanish. This is also consistent with the physical interpretationof the phase since a nontrivial phase corresponds to a superfluid velocity and the velocity
must vanish if the wave cannot be transmitted.

3.5Single Boundary Conclusions
The full set of stationary states of the mean field of a Bose-Einstein condensate,modeled by the nonlinear Schrödinger equation in one dimension, in the presence ofa potential step or pointlike impurity were presented in closed analytic form. Nondecaying solutions were divided into two categories: localized soliton-like solutions, andsolutions that oscillate out to infinity. The localized solutions are of a purely nonlinear
character, as they have no linear analog. On the other hand, the oscillating solutions
were shown to be adiabatically connected to the solutions to the linear Schrödinger
equation.
With a delta function potential, the localized solution can be interpreted as asingle bright or dark soliton trapped by the impurity. In addition, it was shown thatan impurity can also bind a soliton pair. If the impurity is attractive, the natural
repulsion between two dark solitons can be exactly canceled by the attraction of the
impurity, while if the impurity is repulsive, it can balance the natural attraction of inphase bright solitons. Since these solutions are conjected to be stable (see below), they
are excellent candidates for the experimental realization of stationary excited states ofa Bose-Einstein condensate. In addition, the maximum repulsive interaction strengthof the condensate with an attractive impurity that allows for a bound state has been
determined.

46
The oscillating solutions of the NLS, despite being adiabatically connected tooscillating solutions of the Schrödinger equation, have very different properties due tothe concept of an effective potential. For the attractive interaction solution for anevanescent wave decaying under a step, as illustrated in Fig. 3.2(b), the eigenvalue islarger than the effective potential in the regions of high density and is less than theeffective potential in regions of low density. Figures 3.2(c) and (d) show more radical
deviations from the linear solutions since the wave function decay on the lower side of
the step.
It is possible to characterize the general solution of the NLS, Eq. (3.10), in termsof physical quantities such as the mean linear number density, n , mean energy density,
E, and mean momentum density, P. The densities are given by
E(k) +,k
k2 K(k)3B 2 g − b2 A 2= nµ ++ (µ − V )(n − B) ,63
*1

n = B+AE−
2

P = α,(3.57)(3.58)(3.59)

where µ, α, and k are given by Eqs. (3.12), (3.13), and (3.14), respectively. In general,the number, energy, and momentum densities must be calculated separately for the left
and right sides of a boundary. These densities can be used to determine the variables
A, B, and b of Eq. (3.10), leaving only the translational offset, x0 , as a free variable
determined by the boundary conditions. It should be noted that the factor multiplying
A in the mean number density, Eq. (3.57), approaches one half when k approaches zeroand approaches one when k approaches unity. This is to be expected since the mean
number density of a linear wave is given by B + A/2 and the mean number density ofan extremely nonlinear wave is given by B + A. The mean energy density can easily becalculated fromg
E = nµ − ρ2 ,2

(3.60)and so the second and third terms on the right side of Eq. (3.58) are due to this47
nonlinear correction. Since the mean momentum density, Eq. (3.59), is equal to α, themean momentum density must be continuous across the boundary due to the boundary
condition on α, Eq. (3.19). In addition, the momentum density as a function of positionis also given by α, and so the momentum density is equal everywhere.
The healing length, ξ, of the NLS in the quasi-one-dimensional regime, where thetransverse dimensions are trapped by a harmonic potential of frequency ω, is given by
ξ2 =

2
lho
,
8πas n(3.61)

where n is given by Eq. (3.57). Since the mean number density can vary across theboundary, it is possible for the condensate to have a different healing length on eitherside of a boundary. Since the speed of sound in the condensate is inversely proportional
to the healing length, the speed with which phonon-like excitations can travel vary as
they cross the boundary.
While finding the complete set of solutions to the Gross-Pitaevskii equation,
Eq. (3.1), with an impurity or step potential provides much information about thesystem, only stable solutions are experimentally observable. Previous works have examined the stability of stationary states for a constant external potential (see, for example, [80, 79, 60, 61, 74, 95, 166, 176]), as well as for periodic and harmonic potentials(see, for example, [24, 28]). Most studies are ultimately numerical: linear stability canbe solved in a few special cases, while nonlinear stability is analytically intractable. For
a constant external potential, single bright and dark solitons and dark soliton trains arestable. A finite number of bright solitons may form bound states, as for example ordern solitons (n > 1). Bright soliton trains are always unstable. However, they may be
experimentally stable since their lifetime may be longer than experiments, which typically require stability timescales of from milliseconds to seconds. Bright soliton trains
which have a phase difference ∆φ between adjacent peaks such that −π/2 < ∆φ < 3π/2
exhibit this experimental stability, with the lifetime being longer the closer ∆φ is to π.

48
Bright soliton trains with −π/2 < ∆φ < π/2 are unstable but become quasi-periodic intime in a finite system.
Based on these known results from the case of a constant potential, as well as thefrom the stability analysis with an impurity performed by Bogdan et al. [14], the stabilityof an attractive condensate with an impurity is as follows. According to Bogdan, thebound state of two bright solitons, as in Fig. 3.4(b), are stable, since, so long as they
are strongly overlapping, they will be in phase (∆φ = 0) and remain bound to theimpurity. This is also the ground state of the system. However, the kind of solution
shown in Fig. 3.4(a) is unstable [14]. We now conject on the stability of a repulsive
condensate that is not bound and whose density approaches a nonzero constant at
infinity. Localized solutions in the case of repulsive nonlinearity are obviously stablein the cases of Fig. 3.3(a)-(b) and 3.3(e)-(f), since they are the ground states. The
bound pair of dark solitons, illustrated in Fig. 3.3(c)-(d), should be likewise stable,
so long as the impurity is sufficiently strong. All of these solution types, except for
attractive solitons bound by a repulsive impurity, are expected to be experimentallyobservable in finite systems, such as an elongated harmonic trap. An excellent analysisof the stability of solitons pinned with impurities is given by Bogdan, et. al. [14]. The
stability of soliton trains is a less certain issue. Such a stability study presents a subject
for future research.
We emphasize that neither the idea of left and right traveling waves nor that ofreflected plus incident waves apply to nonlinear wave equation defined by the nonlinear
Schrödinger equation. This is important since one cannot create wavepackets from linear
combinations of these solutions. Instead, these solutions already contain the wavepacketlike solutions, or solitons, that are necessary to describe the system; moreover, solitons,
unlike wavepackets, are nondispersive. Time dependent nonlinear scattering remains an
open question that can certainly be addressed via numerical studies. In general, thestationary solutions to the NLS give physical insight into its dynamics, without which

49
numerical solutions may be difficult to interpret. Perhaps more importantly, using thegeneral nature of the solutions to the cases of a step function and an impurity, it ispossible to describe all stationary states to piecewise constant potentials.
In conclusion, we have analytically solved for all stationary solutions to the nonlinear Schrödinger equation with a delta function or a step function potential. Thismodels the steady state behavior of the mean field of a Bose-Einstein condensate in
the presence of an impurity, or of a potential step created by, for instance, a laser
passing over the edge of a razor blade. Novel wave function were found, including solitons trapped by the impurity and the nonlinear analog of transmitted and evanescentwaves.

3.6

A Kronig-Penney Potential and Bloch Waves
We now consider the mean-field model of a quasi-one-dimensional BEC in thepresence of a Kronig-Penney potential,
V # (x) = V0#

+∞
$

j=−∞

δ(x − j d) ,(3.62)

where d is the lattice spacing and V0 is the strength of the potential. When the transverse
dimensions of the BEC are on the order of its healing length, and its longitudinal
dimension is much longer than its transverse ones, the 1D NLS [145, 29] which describes
the stationary states of the mean field of a BEC is given by,1− ∂xx Ψ + g|Ψ|2 Ψ + V (x)Ψ = µΨ .
2(3.63)

Here m is the atomic mass, µ is the eigenvalue, g characterizes the short range pairwiseinteraction, and V (x) is an external potential [119].
In Eq. (3.63), the length is scaled according to the lattice spacing d, and theenergy has been rescaled by π 2 /(2E0 ), whereE0 ≡h̄2 π 22md2

(3.64)

50
is the kinetic energy of a particle with a wave vector equal to that at the boundary of
the first Brillouin zone. The variables in Eq. (3.63) are defined by
1 #
x ,
d
π2 #
µ ,
2E0
π2 #
g ,
2E0
π2 # #V (x /d) ,
2E0

x =
µ =
g =
V (x) =

(3.65)(3.66)(3.67)(3.68)

where the primed variables contain the physical units of the system. The renormalized 1D coupling is, g ≡ nas ωmd2 /h̄, where a harmonic oscillator confinement in thetransverse directions has been assumed with frequency ω. Both attractive and repulsive
atomic interactions, i.e., g > 0 and g < 0, shall be considered. The wave function ororder parameter, Ψ(x, t), can be written asΨ(x, t) =&ρ(x, t) exp[−iµt + iφ(x, t)] ,(3.69)

where ρ(x, t) is the line density and v(x, t) = (h̄/m)∂φ(x, t)/∂x is the local superfluid
velocity.
In addition to the NLS, Eq. (3.63), the normalization of the wave function is given
by,
n=

# 1

ρ(x# )dx# ,(3.70)

0where n is the number of atoms per lattice site. The boundary conditions induced by the
Kronig-Penney potential causes a discontinuity in the derivative of the wave function
across each delta function,lim [∂x ρ(j + /) − ∂x ρ(j − /)] = 4V0 ρ(0) ,"→0(3.71)

where j is an integer. By using the complete set of stationary state solutions to theconstant potential case, we can calculate the full set of Bloch solutions for a lattice. The

51boundary conditions are used to determine the appropriate values of the free variables,
B, b, x0 , and k.
In this section, we use the methods from earlier in this chapter and apply them
to a situation in which a lattice of delta functions is present. As is traditionally donein calculating band structures, we seek Bloch wave solutions of the form
Ψ(x) = eiqx fq (x),(3.72)

where q is the wave number, h̄q is the quasi-momentum and fq (x) has the same period as
the lattice, fq (x) = fq (x + 1). By substituting Eq. (3.72) into Eq. (3.69), one finds thatthe density, ρ, must also have the same period as the lattice and that the wavenumber
and energy per particle can be determined from the density profile by,
q =
E[Ψ]
n

=

# 1

α
dx# .
#
0 ρ(x )
(
)
#
1
1 1
g 4
2
2
dx
|∂x Ψ| + |Ψ| + V0 δ(x)|Ψ| .
n 0
2
2(3.73)(3.74)The quasi-momentum is simply the phase jump across each lattice site.
We will examine the quasi-momentum energy bands. The problem reduces to
one in which the density is symmetric about the center of the lattice sites, x = j + 0.5,
where j is an integer. Due to this symmetry, there are only two possible values for thetranslational offset x0 ,bbx0 ∈ − , K(k) −
2
2
,

-

,

(3.75)

where 2K(k) is the period of the density. The offset forces the density in the center of
each site to be either a minimum or a maximum of the site, depending on the sign ofthe interaction.
It is computationally intensive to include the integral in the wavenumber equation,
Eq. (3.73), with a root finding algorithm. Therefore, one of the parameters, b, k, or
B, is varied while the other two are determined from the number equation, Eq. (3.70),
and boundary condition, Eq. (3.71). The offset is then chosen depending on which

52energy band is being examined. The wavenumber and energy are evaluated from theseparameters and can then be plotted parametrically. In the following section, we discuss
the energy bands for repulsive and attractive interatomic interactions.

3.7

Nonlinear Band Structure
The structure of the energy bands is strongly dependent on the strength and signof the atom-atom interactions, g. In this section, the band structure is examined infour parameter regimes: weakly repulsive, strongly repulsive, weakly attractive, andstrongly attractive. Note that in this chapter, a weakly interacting system is defined by
gn/V0 ≤ 1, and a strongly interacting system is defined by gn/V0 ( 1.
In Fig. 3.7 and Fig. 3.8, the energy bands for specific cases of weak and strong
repulsive interactions are presented, gn = E0 and gn = 10E0 , respectively. The condensates are in a repulsive lattice, V0 = E0 . In Fig. 3.7, the interaction strength is small
and deviations from the linear band structure are small as well. The bands are verticallyshifted higher as compared to a linear system due to the repulsive interactions, which
increase the energy of the system. When the interaction strength is further increased theband structure becomes quite different. Swallowtails [175], or loop structures, appear
in the bands, as in Fig. 3.8. The width of these swallowtails grows as the interaction
strength is increased. Swallowtails are a general feature of a nonlinear system in aperiodic potential [174] and appear for both repulsive and attractive interactions.
The presence of swallowtails has been conjectured to be due to the hysteretic
behavior of the superfluid condensate. A thorough discussion of this topic is given
by Mueller [117] (see also references therein). For a completely free, non-interacting
system, the energy has a quadratic dependence on the momentum, shown as the dashedparabolic curves in Fig. 3.9. Since a change in wavenumber of 2π leaves the system
unchanged, the quadratic dependence is repeated centered around integer multiples of
2π. When a periodic potential is added to the system, bands in the energy are formed,

53

E/n (units of E0)

10
8
6
4
2
0
-1

-0.5

0
0.5
q (units of π/d)1Figure 3.7: Energy per particle as a function of the quasi-momentum for the first threebands of a weakly repulsive condensate (gn = E0 ) in a repulsive lattice (V0 = E0 ). The
noninteracting, linear band structure is given by the dashed curves.

shown as the dot-dash sinusoidal curve in Fig. 3.9. An interacting condensate, however,is a superfluid and can therefore screen out the periodic potential [169]. The energy
band will then appear similar to that of a free particle, shown as a solid swallowtail
curve in Fig. 3.9, until a critical point. At the critical momentum, determined by theaverage condensate sound speed, the energy band terminates. If this velocity is greaterthan the edge of the Brillouin zone, there are then two separate energy minima. This
demands that there be a saddle point separating them and hence the three stationary
states. This discussion can be used to consider condensates with repulsive interactions,but not attractive interactions. We find the appearance of swallowtails for both signsof interaction.
For interacting systems in periodic potentials, there is a minimum interaction
strength for which the swallowtails in the energy bands can exist [117]. In general, thiscan be dependent on both the strength of the potential as well as the band that is being
discussed. For a sinusoidal lattice it was shown that the onset of the swallowtail for

54

E/n (units of E0)

20

15

10

5
-1

-0.5

0
0.5
q (units of π/d)1Figure 3.8: Energy per particle as a function of the quasi-momentum for the first threebands of a strongly repulsive condensate (gn = 10E0 ) in a repulsive lattice (V0 = E0 ).
Note the appearance of swallowtails at the edge of the Brillouin zone in the odd bands
and at the center of the even bands.

the lowest band of a repulsive condensate occurs when the interaction strength and thepotential are equal. For higher bands the relationship no longer becomes analytic [107].
For the Kronig-Penney lattice potential, the critical value for the onset of the swallowtails is not dependent on the band under consideration or the sign of the interaction.
Numerically, we are able to determine that the onset occurs when gn = 2V0 . This
holds true in all cases except for the lowest band of an attractive condensate, as will be
explained later in the section.
The energy bands are slightly different when the condensate is in an attractive
potential, i.e. V0 < 0. At the Brillouin zone boundary, the energy gap between bands
is proportional to V0 for a weakly interacting system. The second band then can crossthe first band, since for an attractive potential the gap size is negative. In contrast, thebands are separated by an energy ρ(0)V0 for a repulsive potential. As the interaction
strength increases, the effects of the potential become less noticeable and the bands are

55

E/n (arb. units)

1

0.5

0
-2

-1

0
q (units of π/d)12Figure 3.9: Energy per particle as a function of the quasi-momentum wave number for anoninteracting system with no potential (dashed curve), a periodic potential and small
or no interaction (dot-dashed curve) and a periodic potential with a large interaction(solid curve).

no longer degenerate. Note that repulsive and attractive sinusoidal potentials create
the same band structure. The difference in the current system arises since there are twolength scales associated with a Kronig-Penney potential, the lattice spacing and delta
function width.
The energy bands for an attractive condensate in a repulsive potential, V0 = E0 ,with a small interaction strength, gn = −E0 , have a qualitatively similar form as for aweakly repulsive condensate, see Fig. 3.10. Note that the attractive bands are, however,lower in energy than the repulsive bands due to the attractive interaction strength. In
contrast to the weakly interacting case, a strongly attractive condensate has several
qualitative differences compared to a strongly repulsive condensate.
In Fig. 3.11, the band structure for a large attractive interaction, gn = −10E0 , ina repulsive potential, V0 = E0 , is illustrated. The swallowtails in the bands are now on
the upper band at the band gaps as opposed to the lower bands at the band gaps as they56

10
E/n (units of E0)

8
6
4
2
0
-2
-1

-0.5

0
0.5
q (units of π/d)1Figure 3.10: Energy per particle as a function of the quasi-momentum wave number for
the first three bands of a weakly attractive condensate (gn = −E0 ) in a repulsive lattice(V0 = E0 ). The noninteracting linear band structure is given by the dashed curves.

were for the repulsive case in Fig. 3.8. The first energy band never has a swallowtail
because the swallowtail must be on the lower portion of the band and below the centerof the band there is no quadratic energy dependence for the swallowtail to follow, seeFig. 3.9.
The higher bands of a strongly interacting attractive condensate can look quite
different than those of a strongly interacting repulsive condensate. After the initial
critical value of the interaction strength is reached, a swallowtail in the band starts toform, as in the third band of Fig. 3.11. As the interaction strength increases, anothercritical value is reached where the width of the swallowtail in q reaches π and runs into
the band edge. Eventually, the wave number that should be less than zero becomes
imaginary due to the form of α. This represents a nonphysical solution and is contrary
to the assumption that the phase is real in Eq. (3.69). Therefore, the band appears as
two separate curves, a loop and a separate line. These are both marked as Band 2 inFig. 3.11. It takes on this appearance since the swallowtail cannot extend lower than57

E/n (units of E0)

5

Band 3-5
Band 2-15
Band 1
-25
-1

-0.5

0
0.5
q (units of π/d)1Figure 3.11: Energy per particle as a function of the quasi-momentum wave number for
the first three bands of a strongly attractive condensate (gn = −10E0 ) in a repulsive
lattice (V0 = E0 ). The first band does not have a swallowtail. However, the second andthird bands do have swallowtails. The swallowtail in the second band appears as a loopwith an unattached curve since an attractive condensate has a maximum width for theswallowtails, which in the case of the second band is a width of π.

the minimum in the adjacent quadratic energy dependence. In general, the nth band
will develop this two-part structure when the swallowtail reaches a width of (n − 1)πsince the swallowtail will then have reached the minimum of the corresponding freeparticle energy. Note that this phenomenon does not occur for a repulsive interaction
since there is no extremum to limit the growth of the swallowtail.

3.8

Density Profiles of Bloch Waves
As density is the primary experimental observable for BEC’s, the change in thedensity profile with wave number is important. Weak interactions create similar changes
in the density profile and energy band structure independent of the sign of the interaction. This is because the interaction energies are less than or of the same magnitude as
the potential. For strong interactions, the density is in general more sharply peaked for

58
an attractive condensate and flatter for a repulsive condensate. However, the generalways in which the densities changes are qualitatively similar.
The densities of the weakly attractive condensate in the first band are shownfor three different quasi-momenta in Fig. 3.12. The solid curve, dashed curve, anddotted curve are the density profiles for q = 0, q = π/2, and q = π, respectively. To
understand the energy bands in terms of the density profile, the three terms of theenergy in Eq. (3.74) should be discussed. The kinetic, interaction and potential energy
per particle are given by the first, second, and third terms of Eq. (3.74). Notice thatthe density at the origin decreases as the wave number is increased. Therefore thepotential energy will also decrease due to the delta function at x = j, where j is an
integer. Because the condensate is attractive, the interaction energy decreases as thewave number is increased since the density becomes more peaked. The kinetic energy
increases as the wave number increases since the variations in the density become larger.
The case of the density variations associated with the first band of a stronglyinteracting repulsive condensate is shown in Fig. 3.13 and is qualitatively similar tothat of the weakly attractive case. The solid curve represents the density profile for
q = 0. The dot-dash curve represents the density when q = π at the bottom of theswallowtail. The dashed curve represents the density when q = 0.46π at the end of
the swallowtail. The dotted curve represents the density when q = π at the top of theswallowtail. The density at the origin again decreases as the wave number increases andtherefore the potential energy also increases . The interaction energy increases as theband is traversed since the density becomes more peaked. The kinetic energy follows
the same qualitative path as the total energy. Of the three, it is therefore the kinetic
energy that has the greatest influence on the energy bands.

59

3
Density (units of 1/d)

(c)
(b)

2

(a)
1

E/n (units of E0)

0
0

1

0.5
x (units of d)

0.3
(c)

0.2

(b)

0.1
(a)

0
-0.1
-1

-0.5

0
0.5
q (units of π/d)1Figure 3.12: Changes in the density of a weakly attractive, gn = −E0 , condensate
associated with different positions on the first band. The solid curve represents thedensity when q = 0. The dashed curve represents the density when q = π/2. The
dotted curve represents the density when q = π. The lower plot shows the corresponding
positions on the first energy band.

3.9

Stability of Bloch Waves
We now study the stability of the Bloch states and determine the stable regionsof the bands. In addition to stable solutions, solutions that have instability times much
longer than experimental time scales can be observed in experiments. Recent studiesof the stability of condensates in a periodic potential have focused on linear energetic

60

Density (units of 1/d)

2(c) (d)

(b)
1

(a)

E/n (units of E0)

0
0

1

0.5
x (units of 1/d)

9
(c)

8
7

(a)

(d)
(b)

6
5
-1

-0.5

0
0.5
q (units of π/d)1Figure 3.13: Changes in the density of a strongly repulsive condensate, gn = 10E0 ,associated with different positions on the first band. The solid curve represents thedensity when q = 0. The dot-dashed curve represents the density when q = π at thebottom of the swallowtail. The dashed curve represents the density when q = 0.46π atthe end of the swallowtail. The dotted curve represents the density when q = π at the
top of the swallowtail. The lower plot shows the corresponding positions on the firstenergy band.

and dynamic stability, also called Landau stability [168, 23, 106, 107, 50]. In contrast,
we consider the full response of the condensate to stochastic perturbations. In order tonumerically simulate the NLS with a periodic potential, a ring geometry was used witha quantized phase such that N q = 2πj, where q is the wave number, N is the number

61of sites and j is an integer. To ensure that the phase quantization does not effect thestability properties, enough lattice sites were used to allow for many rotations of thephase. The outcome of the stability analysis is independent of the number of sites for
sufficiently large number of sites. In most cases, j = 4 was found to be adequate to
extract the correct stability properties.
The delta functions were simulated by single point distortions in the potential
grid. They were also implemented by using boxes of varying widths with areas normalized to create the appropriate potential strength. The size of the boxes did not influence
the stability properties until the width became approximately 10% of the healing length,%ξ ≡ h̄/ 2gn .

(3.76)The NLS was evolved using a variable step fourth-order Runge-Kutta algorithm in timeand a filtered pseudo-spectral method in space. The noise introduced into the simulations comes from the round off error associated with numerical simulations. To ensurethat the form of the noise from round off error did not affect the stability propertiesof the system, initial stochastic white noise of various levels was introduced into the
Fourier spectrum. For levels significantly greater than the round off noise, the stability
times approached those from the round off noise. Introduction of white noise at thelevel in the eighth significant digit produced the same instability times as the round off
noise, which effects the sixteenth significant digit. All simulations were performed overtime scales longer than experimental lifetimes of the BEC, which are on the order of
seconds.
The time at which the onset of instability occurs is determined by the effective
variance in the Fourier spectrum,.1// p (f (p, t) − f (p, 0))2
1
σ(t) ≡ 0,2

2

p (f (p, 0))(3.77)

where f (p, t) is the Fourier component of the wave function at momentum, p, and time,t, and the sum is over the momentum. This quantity determines how different the62
Fourier spectrum is compared to the original stationary state. It vanishes when the two
spectra are identical and approaches unity when there are no Fourier components incommon. When σ(t) reaches 0.5, i.e. 50% of the Fourier spectrum is different than theoriginal, the system is considered to have become unstable.
Unless otherwise noted, for the stability analysis the lattice spacing is given by
d = 1 µm, the approximate length scale with which current optical lattices are created.
In addition, all instability time scales will be given for 87 Rb.

3.9.1

Attractive Atomic Interactions

With an attractive interaction of gn = −E0 in a repulsive potential of V0 = E0(see Fig. 3.10), the lowest energy solution, zero quasi-momentum in the lowest band, hasa lifetime greater than experimental time scales. However, when even a slight harmonic
perturbation to the potential is added to the initial time step the condensate becomes
rapidly unstable. For instance, with a harmonic frequency of 120 Hz, which is approximately the experimental trapping frequency [115], simulations show an instability at
1.5 ms. This is short compared to the lifetime of a BEC [115], but still observable.
When the quasi-momentum of the first band is increased, the stability of thesystem becomes dependent on the effective mass. The effective mass, m∗ , is defined as,m∗ ≡1∂ 2 E/∂q 2.

(3.78)The physical meaning of the effective mass is the mass that the particle would appear
to have if the potential was not being considered [4]. The sign of the effective masscan be transferred to the interaction strength, changing an attractive interaction to aneffective repulsive interaction. Therefore, when the quasi-momentum increases and theenergy band becomes concave down, m∗ < 0, the system enters a regime of stability.
The system remains stable even in the presence of a harmonic perturbation.
For zero quasi-momentum in the second band, the system immediately develops

63
temporally periodic variations in the phase and density. There is an additional instability, occurring on the order of 5 ms, that destroys the periodicity of the system. One
might expect that this part of the band be stable since there is a negative effective
mass, but the oscillations due to the two density peaks per lattice site force the system
to be unstable. The oscillations, although periodic in time, create a larger underlying
instability to grow.
The stability properties of a strongly attractive condensate are similar to thoseof a weakly attractive condensate. The stability of the first band is determined by
the effective mass while higher bands always go unstable. Therefore, for an attractive
condensate, the system of Bloch waves is stable only if there is one density peak perlattice site and the effective mass is negative.

3.9.2

Repulsive Atomic Interactions

Like an attractive condensate, a repulsive condensate only has stable regions on
the first band. For a weakly interacting repulsive condensate, gn = E0 and V0 = E0 , theeffective mass in the first band is positive between q = 0 and q = π/2. The effective mass
becomes negative for larger quasi-momentum, since the energy becomes concave down.
Hence the system becomes unstable. For a wave number of q = 9π/16 the instability
time is 10 ms. This decreases to 2 ms for q = π. In this regime, with negative effective
mass, the actual ground state is an envelope soliton that can spread over many latticesites. These types of states are called gap solitons [101, 76, 105, 45] and only occur ininteracting systems. Figure 3.14 presents the unstable evolution of the weakly repulsive
condensate in the first band with a wave number of q = π in Fourier space. Noticethat the instabilities arise from perturbations around the primary Fourier componentsof the wave function. In Fig. 3.15, the effective variance, σ, is plotted as a function ofevolution time. The system becomes unstable around 2 ms. The second band becomes
unstable in 0.5 ms for q = 0 and in 6 ms for q = π. Therefore, the system is stable in the

64

Figure 3.14: Logarithm of the Fourier spectrum during the time evolution. Time is inunits of h̄/E0 and distance in units of the lattice spacing.

first band with positive effective mass and unstable elsewhere. This is consistent withthe effects that the effective mass has on stability in systems described by the lowest
band DNLS. In a work by Fallani et al. [50], the instability time of a condensate in a
lattice was measured by using an RF-shield to remove the hottest atoms produced by
the heating created in the sample by instability. The loss rate, equal to the inverse of the
lifetime, should then be qualitatively similar to the instability time. Our calculations
are consistent with these experimentally observed loss rates of a BEC in an opticallattice [50].
Due to the presence of the swallowtails, the strongly interacting system provides
different stability regimes. For a repulsive condensate with gn = 10E0 , in a repulsive
lattice with V0 = E0 , the main section of the first band, as well as the lower portion of
the swallowtail, have positive effective mass and remain stable. A small portion at theedge of the band still has a negative effective mass and is unstable. The upper portionof the swallowtail, as discussed in Sec. 3.7, is an energy maximum and is not expected

65

Variance σ

1

0.5

0
0

1
2
2
t (units of md /h)3

Figure 3.15: The time evolution of the effective variance of the momentum density, σ.
Time is in units of h̄/E0 .

to remain stable. Our simulations find that the instability time of the upper portion of
the swallowtail is approximately 0.2 ms, independent of the actual quasi-momentum.
The size of the lattice spacing can influence the time scales for which the system
becomes unstable. In Fig. 3.16, the instability time is presented as a function of thelattice spacing for a repulsive condensate, gn = Ea , in a repulsive lattice, V0 = Ea , witha wave number of q = π, where Ea ≡ h̄2 π 2 /2m(1µm)2 .
There is a minimum instability time as the lattice spacing varies that occurs when
the lattice spacing is approximately twice the healing length, ξ. The instability time
is given by half of the interaction strength time, gn/h̄. When the lattice spacing ismuch larger than the healing length, the density becomes extremely flat except at thedelta functions, where the density deformations take the form of pinned dark solitons.
Since the lattice spacing is large, the dark solitons are far apart and are effectively
noninteracting pinned solitons. Dark solitons are known to be robustly stable [95, 24].
Therefore, for a lattice spacing much larger than the healing length the system becomesInstability Time (units of h/gN)

66

1

0
1

2
Lattice Spacing (units of ξ)

3

Figure 3.16: The instability time as the lattice spacing is varied. Time is in unit of√h̄/gn and distance in units of the healing length, ξ = h̄/ 2gn.

stable. For lattice spacing smaller than twice the healing length, the condensate does
not distinguish between the separate delta functions and, therefore, experiences closer
to a constant potential. In this regime, the kinetic energy becomes much greater than
the interaction energy and potential energy since variations in the density occur on thelength scale of d/2, which is less than the healing length. The system becomes effectively
free and noninteracting and, therefore, approaches stability.

3.10

Analytical Methods in Nonlinear Band Theory
The band structure of a BEC in an optical lattice has previously been studied

both analytically and numerically. An array of techniques have been developed. The
sinusoidal potential, which is the experimental case, has been investigated by expanding
the wave function in a Fourier series and minimizing the energy [175, 43, 173, 107, 106].
Two other potentials have been studied in detail. A Jacobi elliptic potential, which
asymptotically approaches the sinusoidal case, has the advantage that it admits a useful

67
class of exact solutions [23, 24, 22]. On the other hand, the Kronig-Penney potential canbe solved completely and exactly, as we have shown previously [153, 152]. A work by
Li and Smerzi followed a similar line of analysis [102]. We compare these three models
with the aid of the Bloch wave representation. For the Jacobi elliptic potential thisrequires a reinterpretation of the original result [23]. In addition, we explain the longstanding open problem of why certain density offsets in this class of solutions lead to
instability. We find that the form of the potential has no effect on the superfluid andother physical properties of the system. Of the closed-form analytical methods available
for the three potentials, only the Kronig-Penney potential admits a full description of
the band structure.
Of the many methods used in nonlinear band theory, analytical methods allow
for greater flexibility in describing solution types, as they result in general expressions
which describe all parameter regimes simultaneously. In this section, we compare threemethods that use different periodic potentials and interpret their solutions in the Blochwave representation. It is found that they all produce the same band structure for thelowest band. Hence, the exact form of the potential is unimportant.

3.10.1

Solution by Cancellation

The first analytical method for obtaining solutions to the NLS with a periodic
potential in the context of a BEC was given by Bronski, et. al. [23] for a Jacobi elliptic
potential of form
V (x) = −V0 (1 − 2 sn2 (2K(k)x, k)) ,(3.79)

where sn is one of the Jacobi elliptic functions [1]. The Jacobi elliptic functions are
generalized periodic functions characterized by an elliptic parameter k ∈ [0, 1]. They
approach circular and hyperbolic trigonometric functions as k → 0 and k → 1, respectively. For k not exponentially close to unity, the potential is similar to the sinusoidal
case. As k approaches unity, the period of the lattice becomes much greater than the

68width of the localized variations of the potential, thereby approaching a lattice of delta
functions, i.e., a Kronig-Penney potential. The period of the lattice is 2K(k) ∈ [π, ∞]where K(k) is a complete elliptic integral of the first kind [1].
Bronski et al. were able to show that it is possible to choose a suitable ansatz
for the density such that the nonlinear term cancels the potential term, thereby leaving
a system that is effectively free. They found a special class of solutions with a period
equal to that of the lattice, where the density was assumed to be of the form
ρ(x) = A sn2 (bx, k) + B .

(3.80)

The following conditions must then be met:A =b2 k2 + 2V0
,
g
2BV0
),
A
1 2
BV0
(b + Ag + 3Bg) −.
2
A(3.81)

α2 = B(B + A)(b2 + Bg −(3.82)µ =(3.83)

Machholm et al. [107] showed that this solution corresponds to the edge of theBrillouin zone for the exactly sinusoidal case, k = 0. If the parameters of the potential,
V0 , k, and b, and the interaction strength gn are set, the wave function is completely
determined without any free parameters. Nevertheless, the solutions are of more generaluse. In fact, although not described in this way by the original authors, this family ofsolutions can be used to map out the whole lowest energy band. It is possible tochange the elliptic parameter k to get a spectrum of solutions. Since the Jacobi elliptic
functions closely approximate the trigonometric functions for all k except exponentiallyclose to unity, there is a wide range of values for k where the potential is approximately
sinusoidal. In addition, the sign of the potential coefficient V0 does not significantly
change the form of the potential and so this is another parameter that can be changed
to determine the energy bands. Therefore, when the elliptic parameter is varied, andthe solutions determined via Eqs. (3.80)-(3.83), the complete lowest energy band can

15

15

E/n (units of E0)

(a)

(c)

10

10

10

5

5

5

0
-1

0
0
q (units of q0)

1

-1

69

15
(b)

0
0
q (units of q0)

1

-1

0
q (units of q0)

1

Figure 3.17: Band structure for a strongly repulsive condensate, gn = 10, in a periodic
potential. Shown are results obtained via the analytical methods of (a) Bronski et al. [23]with a Jacobi elliptic potential (|V0 | = 1), (b) Machholm et al. [107] with a sinusoidal
potential (V0 = 4), and (c) Seaman et al. [153] with a Kronig-Penney potential (V0 = 1).

be analytically determined for arbitrary interaction strength.
An example of the band structure for strong nonlinearity is shown in Fig. 3.17(a).
The interaction strength is chosen to be a factor of ten larger than the strength of thepotential, causing the appearance of swallowtails. Only the lowest energy band canbe extracted with this class of solutions. This shows that these simple exact solutions
are sufficient to describe measurable properties of the condensate such as breakdown ofsuperfluidity for critical values of the nonlinearity [117, 50].
In the work of Bronski et al. [23, 24], the linear stability properties of the solutions
were proved for the case of constant, or trivial phase. No such proof was discovered forthe case of non-trivial phase solutions. The latter are, for example, the only stable onesin the exactly sinusoidal case of k = 0. It was found numerically that non-trivial phasesolutions for repulsive condensates with large offsets B were stable, while for smaller B
they were unstable. No explanation for this stability property was found. After we recast
the solutions in Bloch form and used them to determine the band structure, the stability
properties become immediately apparent. Solutions on the upper edge of the swallowtail
are known to be unstable, as they represent an energy maximum [117, 107, 153]; these70
correspond to small B. Solutions on the main part of the band and the lower edge of theswallowtail represent an energy minimum and are stable; these correspond to large B.
For instance, in Fig. 3.17(a), the main part of the band from q = 0 to q = 1 corresponds
to B ∈ [0.559, 0.509]; the lower edge of the swallowtail near the edge of the Brillouin zone
corresponds to B ∈ [0.509, 0.277]; and the upper edge, which is unstable, corresponds
to B ∈ [0.277, 0]. We note that, in comparing with Bronski et al., they normalized thewave function to gn rather than having gn as the coefficient of the nonlinearity in the
NLS. In their numerical studies of stability, they did not hold their normalization (ournonlinear coefficient) fixed. In Fig. 3.17(a), we fix the nonlinearity to be gn = 10|V0 |.
3.10.2

Solution by Three Mode Approximation

In a subsequent work by Machholm, et. al. [107], a more complete method todetermine the nonlinear band structure was introduced that used a sinusoidal potential
V (x) = 2V0 cos2 (πx) .

(3.84)

The wave function was expanded in a Fourier series for a particular quasi-momentum
and the energy was minimized. This method uses the exact physical form of the typicalexperimental lattice potential that a BEC is held in but requires the root finding of many
free parameters. Machholm et al. showed that one can obtain analytical solutions byusing a three-mode Fourier spectrum:Ψ(x) =√

neiqx (a0 + a1 ei2πx + a−1 e−i2πx ) ,(3.85)

where a0 , a1 and a−1 are real coefficients. Due to the normalization condition on the
wave function, there are two real free parameters of the solution.
Using the three-mode expansion of the wave function, it is possible to extract thefirst two bands and the lower part of the third band, as shown in Fig. 3.17(b). Notice
that, unlike with the previous method, full information from the first two bands can

71
be extracted, instead of only the lowest band. However, the three-mode approximation
overestimates the width of the swallowtails [107]. If more Fourier components are included, for a total between five and ten, the first three bands can be described to within
1% accuracy. For higher bands even more Fourier components are needed. This method
was later extended to include period-doubled states [106].
In addition to determining the band structure, Machholm et al. were able todetermine several analytical expressions associated with the stability properties of thecondensate. In particular, the interaction strength at which swallowtails first appear in
the lowest band, and the width of the swallowtails for a given interaction strength were
determined analytically. Studies of the stability of the condensate to small perturbations
were also performed.

3.10.3

Solution via a Piece-wise Constant Lattice

In the previous section, a lattice of delta functions, a Kronig-Penney potential,
was used
V (x) = V0

+∞
$

j=−∞

δ(x − j) .

(3.86)

We presented the complete set of Bloch waves solutions in analytical form by solving
the piece-wise-constant potential case and using the appropriate boundary conditions
to fix the parameters of the solution. The potential, however, is quite different fromthe experimentally created sinusoidal lattice. The optical lattice potential is composedof a single Fourier component while the Kronig-Penney potential is a comb of equallyweighted Fourier components. We proved that the most general form of the density
over a finite interval of constant potential is given by
ρ(x) = B +k2 b2 2sn (bx + x0 , k) ,g

(3.87)

72
where again sn is one of the Jacobi elliptic functions. The chemical potential µ andphase prefactor α are given by1µ = [b2 (1 + k2 ) + 3Bg] ,2(3.88)

α2 = B(k2 b2 /g + B)(b2 + Bg) .

(3.89)The energy bands can be determined by varying one of the parameters in the wave
function, such as B or b, and determining the quasi-momentum and the energy from
Eqs. (3.73) and (3.74). In particular, the parameter scaling b was varied since it is closely
related to the quasi-momentum. The offset B and elliptic parameter k are determinedby number conservation and the boundary conditions across the delta functions.
The first three bands for the Kronig-Penney lattice are presented in Fig. 3.17(c).
Although only the first three bands are presented, our method can be used to determine
the higher bands with no additional computational intensity, unlike the technique ofSec. 3.10.2. For all energy bands, all that is required computationally is the root findingof two parameters, the offset B and the elliptic parameter k. The stability propertiesof the condensate were then numerically determined by dynamically evolving an initialstate composed of the stationary solution plus a small amount of white noise. It was
found that the stability of the energy bands depended on the interaction strength.
For weak nonlinearity, the bands become unstable for quasi-momentum greater thanapproximately π/2. For strong nonlinearity, the tops of the swallowtails are unstable
but the remainder of the first band is stable. The second and higher bands are always
unstable.
In spite of the differences in the potentials, all three methods find a similar nonlinear band structure. In addition, the stability properties of the three models are also
the same. Therefore, to study the main physical aspects of this system, any method
can be used. All three methods are useful in determining the first band. For the fullband structure, only our method is analytically and numerically tractable.

73

3.11

Periodic Potential Conclusions
The full set of Bloch wave stationary states of a Bose-Einstein condensate ina Kronig-Penney lattice potential, with period commensurate with the lattice, havebeen presented analytically for both repulsive and attractive interactions. The quasimomentum energy bands were found to exhibit a cusp at the critical interaction strength
gn = 2V0 , where g is the interatomic interaction strength, n is the linear density, and
V0 is the lattice potential strength. For larger interaction strengths, swallowtails formin the bands. These swallowtails have the same qualitative form as for a sinusoidal
potental and exhibit the same stability properties.
Both attractive and repulsive condensates were found to be dynamically stable
only in the first band and when the effective interaction, sgn(m∗ )g, was positive. Also,
even in the first band, the upper edges of swallowtails were always unstable. Therefore,
for an attractive condensate, the only stable Bloch states exist in the first band between
the wave number where m∗ becomes negative and q = π. A repulsive condensate is only
stable in the first band from q = 0 to the wave number where m∗ becomes negative.
Higher bands are always unstable, for both attractive and repulsive condensates. Whensolutions became unstable, our numerical studies consistently observed that the instabilities originated around the primary Fourier components of the wave function. This
is in agreement with the formal proof of the instability of constant phase Bloch-wave
solutions by Bronski et al. [24]. The instability time was found to be a function of thelattice constant. If the delta functions are spaced either smaller than or much larger
than the healing length of the condensate, the solutions had instability times longer thanlifetime of the BEC. Thus experiments could access formally unstable sections of theenergy bands, and, by controlling the ratio of the healing length to the lattice constant,directly observe the dynamics of instability. The results of our stability analysis are
consistent with the experimental work performed by Fallani et al. [50], in which the loss

74
rate, the inverse of the lifetime, was determined by removing the hottest atoms with an
RF-shield.
We also found that for a repulsive condensate with a swallowtail, almost the entire
energy band is concave up and, hence, the effective mass is nearly always positive. This
is in contrast to a weakly repulsive condensate, where the concavity of the energy bandchanges, creating a region of negative effective mass. Therefore, there is a maximum
interaction strength for which gap solitons [101, 76, 105, 35, 45] can be formed since
they require a negative effective mass. The maximum interaction energy is given by
gn = 2V0 , the strength at which swallowtails appear.
It should be noted that the solutions given in the chapter map onto both thelinear Schrödinger equation and discrete nonlinear Schrödinger equation (DNLS) limits
if the proper procedure is followed. In the limit that
lim k2 /g → A/b2

g,k→0

(3.90)

the linear dispersion relation for the linear Schrödinger equation is recovered. In thelimit that
V0
(1
gn

(3.91)i∂t ψj = /j ψj + Jj (ψj+1 + ψj−1 ) + Uj |ψj |2 ψj ,(3.92)

one obtains the DNLS,

where /j , Jj , and Uj are obtained via the Wannier formalism [133].
Stationary solutions to the NLS with a Kronig-Penney potential need not take theform of Bloch states. Solutions with a period which is an integer multiple of the lattice
period have been shown to exist for the sinusoidal potential [106] and are expectedalso to be present for the Kronig-Penney potential. Envelope solutions, such as gap
solitons, also play an important role in other systems modelled by the NLS and havebeen observed in BEC’s [45]. The analytic methods which we have described here are
equally applicable to these solution types and form the subject of future study [153].

75

3.12

Two Color Lattice
There has been great interest in the study of Bose-Einstein condensates (BEC’s)from both the experimental and theoretical perspectives since BEC’s were first createdin 1995 [6, 39, 19]. In particular, the examination of BEC’s in periodic potentials inthe superfluid phase has yielded many intriguing phenomena, such as gap solitons [45,7], the appearance of swallowtails or loops in the band structure, pulsed atom lasers
and demonstrations of their phase coherence [5, 70], and matter-wave diffraction [121].
BEC’s in periodic potentials, unlike other solid state systems, have the advantage thatthe lattice geometry and interatomic interactions are highly controllable [136, 81]. We
examine the mean field Bloch states of a BEC in one-color and two-color potentials forarbitrary interaction and potential strengths. The full nonlinear band structure of thesystem is then determined for a two-color potential. We show that period-doubled statesin the usual one-color lattice are directly connected with Bloch waves in the two-color
case. Period-doubled states have a periodicity which is twice that of the underlying
lattice. This allows for a novel interpretation of swallowtails, a key physical property ofthe band structure of a BEC. The study of period-doubled states has been importantin other systems, for instance in optics [144, 158], because it offers experimental accessto the period-doubling route to chaos.
The two-color lattice has already received some attention in the literature. Rothand Burnett showed that when the period of the second lattice is much longer than
the first, so that it forms an envelope, new quantum phases are introduced into thequantum problem [138]. We will consider the superfluid phase only, which is obtained
when the lattice height is on the order of or smaller than the chemical potential. Veryrecently, Louis et al. have studied gap solitons in a two-color lattice for which thesecond lattice has half the period of the first [104]. We study the same potential,
but for Bloch waves, rather than gap solitons. Bloch waves have the same period

76
as the lattice, whereas gap solitons are envelope solutions. Using the Kronig-Penney
potential, we link the period-doubled states of the one-color lattice to Bloch waves inthe two-color case. Adiabatically turning off one of the frequencies in a two-color lattice,
one can clearly observe the origin of swallowtails. Unlike previous explanations of thisintriguing property of nonlinear periodic systems [117], ours encompasses both repulsive
and attractive BEC’s.
Experiments on BEC’s in optical lattice potentials proceed as follows. Alkalimetal bosonic atoms are cooled to the quantum degenerate regime. The interference
pattern of two counter propagating lasers is used to create a sinusoidal potential shift
caused by the ac Stark effect induced by the dipole interaction with the laser fieldon the atoms’ center of mass motion [110, 41]. A small frequency detuning of oneof the lasers allows for examination of the different quasi-momentum states and theirstability properties by creating a traveling wave interference pattern moving at thevelocity v = (λ/2)δν, where λ is the wavelength of the first beam and δν is the detuning [49, 124, 50]. The two-color lattice can then be created by the superpositionof two lasers with frequencies which differ by a factor of two. This can be achieved
using second harmonic generation with nonlinear crystals. Complex lattice configurations with several lasers have already been performed experimentally [68]. We present
a dynamical study of key observables for the two-color lattice, such as the instability
time, that can be experimentally investigated with minor changes to current apparatus.
There have been many theoretical studies of stability properties of BEC’s in one-color
lattices (see [23, 24, 22, 164, 76, 89, 105], to name a few), though not of the two-colorcase.
The band structure is modified from the well-known solutions of the linear Schrödinger
equation with a periodic potential. There are two physical regimes [153]. For |g|n ≤ V0 ,which we term the regime of weak nonlinearity, the linear band structure is simply perturbed up or down, depending on the sign of g. For |g|n ( V0 , which corresponds to

77

(a)

(b)

(c)
Figure 3.18: A sketch of the potential is illustrated for the three cases of Figs. 3.19
and 3.20, with appropriately normalized boxes representing delta functions: (a) a onecolor lattice of period d and ∆ = 0; (b) a two-color lattice with a small difference in thepotential strengths, period 2d and ∆ = 0.2; (c) a one-color lattice with period 2d and∆ = 1. Adiabatically tuning the system from (a) to (c) allows for a novel explanationof the appearance of swallowtails in the band structure.

the regime of strong nonlinearity, the bands wrap back around on themselves to form
loop structures, or swallowtails. These swallowtails have previously been described interms of the superfluid screening properties of the condensate [117]. This explanationis only valid for repulsive condensates, despite the fact that swallowtails also appear in
the attractive case. This motivates the need for another explanation, which we provide
in following section.

3.13

The Two Color Lattice and Formation of Swallowtails
Current experiments with BEC’s utilize only one-color lattices in the weakly non-linear regime. Therefore there has been no experimental observation of the swallowtail
structure, which requires strong nonlinearity. Two-color lattices, on the other hand,allow one to produce swallowtails even for weak nonlinearity, as we shall show in thefollowing. Moreover, adiabatically tuning the two-color lattice allows for an explanationof swallowtails for both repulsive and attractive BEC’s.
A two-color lattice is produced by adding a second frequency component of twice
the fundamental frequency, or half the period. Such a straightforward modificationof existing experiments can be made via second harmonic generation with nonlinear

78
crystals. In Sec. 3.10 we showed that the exact form of the potential is unimportant.
Therefore, to study the two-color lattice, we use a Kronig-Penney-like potential withdelta functions of two different strengths:
V (x) = V0+∞$[(1 − ∆)δ(x − 2j) + (1 + ∆)δ(x − 2j + 1)] ,(3.93)j=−∞where ∆ ∈ [0, 1]. The strengths of the two sublattices are (1 − ∆)V0 and (1 + ∆)V0 ,respectively. The form of Eq. (3.93) allows the potential to be tuned continuously
from the one-color lattice of strength V0 and period d (∆ = 0) to the one-color lattice ofstrength 2V0 and period 2d (∆ = 1). This holds the average potential strength constant.
For 0 < ∆ < 1 the potential is two-color. Figure 3.18 sketches the progression between
the two lattices, where the delta functions are represented by square functions of fixedwidth and varying heights. We will show how period-doubled solutions of the ∆ = 0
lattice map onto Bloch-wave solutions of the ∆ = 1 lattice, including the swallowtails.
The stationary states of the NLS with the potential of Eq. (3.93) are determinedin the same manner as was described in Sec. 3.10.3, with the addition of a second
boundary condition. Period-doubled solutions can be obtained via the transformation
Ψ(x) = fq (x) exp(i2qx), while Bloch solutions were found via Ψ(x) = fq (x) exp(iqx).
In general, one can obtain solutions of the form Ψ(x) = fq (x) exp(imqx), with m an
integer; we focus on the period-doubled case, m = 2, and the Bloch wave case, m = 1.
We first illustrate period-doubled solutions for the lattice of Fig. 3.18(a), ∆ = 0.
There are two types of period-doubled solutions, trivial and non-trivial. The trivial
period-doubled solutions simply reproduce the Bloch waves. That is to say, solutions ofperiod d are also trivially solutions of period 2d. These are illustrated in Fig. 3.19(a)as thin curves for the case of weak nonlinearity. Period-doubled solutions extend acrossa Brillouin zone of half the quasi-momentum-domain as that of Bloch waves. Thus
the domain of Fig. 3.19 (q ∈ [0, 1]) is half that of Fig. 3.17 (q ∈ [−1, 1]). The bands
are required to be symmetric around the center of the Brillouin zone. One reflects

79

2
E/n (units of E0)

a)

1.5

1
0

0.5

1

0.5

1

0.5
q (units of q0)

1

E/n (units of E0)

2
b)

1.5

1
0
2
E/n (units of E0)

c)
1.5

1

0.5
0

Figure 3.19: Band structure for a two-color lattice with weak nonlinearity, gn = V0 .
(a) Shown are period-doubled solutions for a one-color lattice of period d and ∆ = 0(Thin black curve: trivial solutions; thick blue curve: non-trivial); (b),(c) Bloch-wave
solutions with period 2d (Solid curve: lowest band; dashed curve: second band). (b)A two-color lattice with ∆ " 1. Note that the central swallowtail is derived from thenon-trivial period-doubled solution of panel (a). (c) A one-color lattice with ∆ = 1.
The swallowtail has disappeared due to the weak nonlinearity. The lattices associatedwith (a)-(c) are sketched in Fig. 3.18.

80
the band in the right half of the Bloch Brillouin zone around its center at q = 0.5 toobtain the trivial period-doubled solutions. This leads to the two branches which makethe form of an ‘x’ in Fig. 3.19(a) (thin black curves). The non-trivial period-doubled
solutions are shown as a thick blue curve in Fig. 3.19(a). These form a saddle between
the trivial solution branches. The form of the density for the non-trivial period-doubled
solutions is shown as the solid curve in Fig. 3.21. The non-trivial solutions are causedby the nonlinearity in Eq. (3.1); they have no analog in the linear Schrödinger equation.
Comparing Fig. 3.19 to the swallowtails illustrated in Fig. 3.17, we make the following
conjecture: non-trivial period-doubled solutions for a lattice of period d appear as thesaddle of the swallowtail for Bloch waves for a lattice of period 2d.
This conjecture is supported by tuning ∆ from zero to unity. We illustrate thistuning in the subsequent panels of Fig. 3.19. Figure 3.19(b) shows the case ∆ " 1.
Figure 3.19(c) shows the endpoint, ∆ = 1. The shape of the lattice in the panels (a)(c) is sketched in Fig. 3.18. In Fig. 3.19(b), the trivial and non-trivial period-doubled
solutions of Fig. 3.19(a) separate into two bands. The lower band (solid curve) consistsof the lower part of the trivial period-doubled solutions, with the non-trivial perioddoubled solutions forming the upper edge of the swallowtail. The upper band (dashedcurve) consists of the upper part of the trivial period-doubled solutions. One observes
that there is a small gap between the two bands. Moreover, unlike in the case of the
one-color lattice, one obtains a swallowtail even for weak nonlinearity. As ∆ is increased,
this gap increases. Figure 3.19(c) shows the case ∆ = 1. The system is now again in a
weakly nonlinear regime with a one-color lattice. Thus the swallowtail disappears.
One observes that, in general, the first two Bloch-wave bands of a one-color latticeof period 2d (panel (c)) can be derived from the trivial period-doubled solutions of alattice of period d (panel (a)) via a two-color intermediate lattice (panel (b)). This
follows from the fact that the second Bloch-wave band has two density peaks per site.
In the strongly nonlinear regime, the swallowtail forms in a similar manner but

81

3

E/n (units of E0)

(a)

2

1
0

0.5

1

0.5

1

0.5
q (units of q0)

1

3

E/n (units of E0)

(b)

2

1
0
3

E/n (units of E0)

(c)

2

1
0

Figure 3.20: Same as Fig. 3.19 but for strong nonlinearity, gn = 10V0 . In (b), (c) notethat the upper edge of the central swallowtail in the lower band is derived from the nontrivial period-doubled solution of panel (a) (blue thick curve). The lattices associatedwith (a)-(c) are sketched in Fig. 3.18.

82

Density (units of 1/d)

2

1

0
0

1
q (units of q0)

2

Figure 3.21: The condensate density profiles are illustrated for three cases: ∆ = 0, nontrivial period-doubled solutions (solid curve); ∆ = 1/2, upper edge of the swallowtail(dashed curve); and ∆ = 1, the lowest band (dotted curve). The parameters are q = 0.5
and gn = V0 , which is the regime of weak nonlinearity, as in Fig. 3.19.

does not disappear for ∆ = 1. We illustrate this sequence in Fig. 3.20, again for ∆ = 0,
∆ " 1, and ∆ = 1. One clearly sees that the upper edge of the swallowtail in the
one-color lattice of period 2d is derived from the non-trivial period-doubled solutions of
the one-color lattice of period d. The lower edge of the swallowtail is derived from thetrivial period-doubled solutions.
This introduces a novel perspective on the origin of the swallowtails in the nonlinear band structure. In the figures we illustrated the case of a repulsive condensate,as that is the most common experimentally. The same argument holds for attractive
condensates, where the swallowtails form below, rather than above the bands. In particular, in our study of the one-color Kronig-Penney potential [153], we found it intriguingthat no swallowtail formed on the lowest band for the attractive case.
Previous interpretations of swallowtails were based on the superfluid screening
properties of the condensate [117]. In this argument, the condensate sees the quadratic
free particle dispersion up to the sound speed. When written in the form of the Blochansatz, these quadratic curves repeat in each Brillouin zone. When the sound speed issuch that the maximum quasi-momentum of each curve overlaps with the curve from83
the adjoining Brillouin zone, one obtains swallowtails. This does not apply to attractive
condensates, since the swallowtails form on the lower edge of the bands and not at allon the lowest band.
In contrast, our argument based on period-doubling applies to both repulsive
and attractive condensates. For attractive nonlinearity, the non-trivial period-doubled
states form a convex saddle on the lower edge of the ‘x’ of Fig. 3.19(a), rather than theupper edge. As ∆ is tuned from zero to unity, the upper part of the ‘x’ separates fromthe lower part, carrying the swallowtail with it. Thus no swallowtail can form on thelowest band.
Thus, we understand swallowtails to be adiabatically connected to period-doubled
solutions of a lattice of half the period. They originate in non-trivial period-doublingbrought about by the nonlinear term in the NLS.
Finally, since the density is a key observable in experiments on BEC’s, in Fig. 3.21
density profiles are illustrated for weak nonlinearity and ∆ = 0, 1/2, 1. The plot is madefor the quasi-momentum q = 0.5 solution that lies at the center of the non-trivial perioddoubled solution (thick blue curve) in Fig. 3.19(a), at the top of the central swallowtails
on the lowest band in Fig. 3.19(b), and on the lower band in Fig. 3.19(c).

3.14

Stability Properties and Comparison with Experiments
Several important experiments have been performed with superfluid BEC’s in

optical lattices. In particular, an experiment by Fallani et al. [49, 42] studied the loss
rate of the BEC as the quasi-momentum was varied. This loss rate is expected to be
monotonically related to the instability time of the condensate. This observable canbe modelled theoretically by adding white noise to an initial stationary state on thelattice in simulations. In this way, we investigate the nonlinear stability properties ofBloch-wave solutions to the two-color lattice. Note that unstable solutions which have
lifetimes longer than experimental timescales will appear experimentally stable.

84

-1

γ (units of ms )

0.75

0.5

0.25

0
0

0.2

0.4
0.6
q (units of q0)

0.8

1

Figure 3.22: Shown is the inverse of the instability time γ = Ti−1 of stationary states in
the lowest band, as a function of the quasi-momentum. White noise was added to theinitial state and the wave function was evolved numerically for the NLS with a one-color
lattice and weak nonlinearity, gn = V0 .

In our simulations, Eq. (3.1) was dynamically evolved using a variable step fourthorder Runge-Kutta algorithm in time and a pseudospectral method in space. We chose
periodic boundary conditions in one dimension, with a sufficient number of sites so thatthe results were independent of the ring circumference. Two schemes for delta functions
on a grid were considered: point defects, and narrow square potentials covering several
grid points. For square potentials less than one tenth the lattice period d, differences
between the two schemes were negligible. Simulations were performed over time scales
on the order of those relevant to the experiments, namely, hundreds of milliseconds.
The lattice spacing d is taken to be 1 µm and the atomic mass to be that of 87 Rb.
The stability of the system is quantified by the variance in the Fourier spectrum,
σ(t) ≡1(f (p, t) − f (p, 0))2
1,
2 (f (p, 0))2

(3.94)

where f (p, t) is the Fourier component of the wave function at momentum p and time t.
The variance approaches unity when the wave function has no component in common
with the initial state and zero when then the wave function it is identical to the initialstate. We define the instability time Ti to occur when the variance has increased toσ = 0.5.

85
A clear experimental observable is the rate of loss of atoms from the system, since
the density of the sample can be imaged at successive times. This loss rate is expected
to be monotonically related to the inverse time in which instabilities arise. Therefore,
we present two relevant studies of the instability times of the condensate. First, thegrowth rate of the instability,γ ≡ 1/Ti ,(3.95)

is determined for a one-color lattice and weak nonlinearity. This relates directly to theexperiment of Fallani et al.. Second, the lowest positive quasi-momentum for which thesystem becomes unstable is evaluated as a function of ∆ for the two-color lattice.
Figure 3.22 presents the results of the first study. The solutions are stable until
the quasi-momentum where the nontrivial and trivial period-doubled solutions connect
in the band structure (see Fig. 3.19(a), as well as Ref. [106]). Stability can also be
determined from the effective mass of the system,m∗ ≡1∂ 2 E/∂ 2 q.

(3.96)

If the effective mass becomes negative, the NLS becomes effectively attractive. Since theground state of an attractive condensate is localized, solutions with negative effective
mass m∗ < 0 are unstable. The point where the nontrivial and trivial period-doubled
solutions connect corresponds to a change in the sign of the effective mass from positive
to negative. Our figure qualitatively replicates the results found by Fallani et al. [50].
Figure 3.23 presents the second study. In particular, the smallest quasi-momentum
for which the periodic system first becomes unstable in a time less than 100 ms is plotted. The two end points at ∆ = 0 and ∆ = 1 are easily determined, since those systems
are one-color lattices. As per our discussion in the previous paragraph, the Bloch wave
solutions become unstable at the point where the non-trivial period-doubled solutions
connect to the Bloch-wave band. The two one-color lattice systems have the same ratio
of nonlinearity to potential strength, since for ∆ = 0 the potential strength is V0 and86

q (units of q0)

0.6

0.5

0.4

0.3
0

0.5

Δ (units of E0)

1

Figure 3.23: The lowest positive quasi-momentum for which the system becomes unstable. The parameter ∆ is varied, tuning from a lattice with period d to one of period 2d.
The results of simulations of the NLS with weak nonlinearity, gn = V0 , are shown.

the number of atoms is n, while for ∆ = 1 they are 2V0 and 2n. Thus the existence ofthe non-trivial period-doubled states scales with the length of the Brillouin zone. The
absolute quasi-momentum where instability occurs then differs by a factor of two, sincefor a lattice of period 2d the Brillouin zone is half the length as for d. In Fig. 3.23,
this is q ∼ 0.6 for ∆ = 0 and q ∼ 0.3 for ∆ = 1. In the case of the linear Schrödinger
equation, the endpoints would be 0.5 and 0.25, respectively. Note that our simulations
are accurate to within a few percent.
The intermediate points in Fig. 3.23, where 0 < ∆ < 1 and one obtains a twocolor lattice, can be understood by again considering where the effective mass becomes
negative. As the energy band from the ∆ = 0 lattice separates into the first two bandsof the ∆ = 1 lattice, the upper edge of the swallowtail shrinks. Recall that stateson the lower edge of swallowtails are stable [153]. As ∆ is increased, the swallowtail
disappears, since the nonlinearity is weak. The quasi-momentum at which the bandchanges from concave up to concave down, i.e. where the effective mass changes sign,
moves from right to left, as can be seen in the descending data points of Fig. 3.23. A
similar argument holds for the case of strong nonlinearity.

87

3.15

Two-Color Lattice Conclusion
We have presented a novel way of understanding the appearance of swallowtails

in nonlinear band structure. Swallowtails in a lattice of period 2d are adiabaticallyconnected to the non-trivial period-doubled solutions of a lattice of period d. Nontrivial period-doubled solutions are caused by the nonlinearity in the mean field of thecondensate [106]. They do not appear as solutions to the linear Schrödinger equation andtherefore have no analog in linear band theory. Our way of understanding swallowtails is
valid for both repulsive and attractive condensates, unlike previous explanations [117].
We used the two-color lattice to adiabatically connect the one-color lattice of
period d and 2d. We showed that swallowtails appear in the band structure of thetwo-color lattice even for weak nonlinearity, unlike in the one-color case. We then madeexplicit predictions for the onset of instability in the band structure for the one- and twocolor lattices based on numerical simulations. Instability is experimentally observableas an increased loss rate due to heating of the condensate [50], and weak nonlinearity
is the present regime of experimental investigation.
Our studies utilized a two-color Kronig-Penney potential. In order to justify theuse of this model we compared analytic solutions to the one-color Kronig-Penney [153],
Jacobi elliptic [23], and sinusoidal potentials [107]. We showed that each of these models
resulted in a similar band structure. Moreover, by putting the exact solutions to theJacobi elliptic potential in Bloch wave form, we were able to resolve a longstanding
question concerning density offsets and stability in the work of Bronski et al [23, 24].
We note that there are other solution types than the ones we have considered. For
instance, there are gap solitons [105, 104, 109] and time-periodic solutions [129, 128].
These can have novel features in the two-color lattice [104]. We emphasize that we havetreated only the superfluid phase of cold bosons on a lattice [52, 64].
Nonlinearity often introduces solution types which do not appear in linear equa-

88tions [84]. For instance, the plane wave solutions to the particle-on-a-ring problem,
which are well known from linear quantum mechanics, have 2j nodes, where j is a positive definite integer. When attractive nonlinearity is added to this problem, so that one
has the attractive NLS, the constant density, or j = 0 solution undergoes a bifurcation
for a critical value of the nonlinearity. The ground state spontaneously develops a peaksomewhere on the ring, while the constant density solution remains a stationary but
excited state [27]. The NLS with a periodic potential has a similar feature. For a criticalvalue of the nonlinearity, the bands form a cusp which then bifurcates into an upper andlower portion of the band, i.e., the swallowtail. The non-trivial period-doubled solutions
are caused by this property.
In our previous work [153], we gave an instance of the same effect for the attractive
NLS with a single repulsive impurity, i.e., a delta function with positive coefficient.
There, for zero or weak nonlinearity there is no bound state and the solution is a pairof solitons propagating away from the impurity. At a critical value of the nonlinearity,
this bifurcates into a travelling wave solution and a bound state. The latter occurs
due to the dominance of the attractive nonlinearity over both the kinetic energy term
and the repulsive effect of the impurity. This is the idea underlying the non-trivial
period-doubled solutions on the lattice, which, in our Kronig-Penney formulation, havean “impurity” located in the center of each density peak.
Interestingly, the symmetry breaking in the above mentioned case of the nonlinear
analog of the particle-on-a-ring problem occurs at exactly the same point even in the
many body quantum field theory [88]. In this sense, it would be intriguing to determine
in the many body quantum theory of weakly interacting bosons with an external periodic
potential (a) whether swallowtails occur, and (b) if they do so at the same critical valueof the nonlinearity as predicted by the mean field approximation.

Chapter 4

Atomtronics

A collection of ultracold atoms subject to a spatially periodic potential can exhibit behavior analogous to electrons in a lattice. This fact has been established in
an impressive series of experiments with Bose-Einstein condensates and Fermi gases in
optical lattices [5, 31, 62, 116, 41, 114, 135, 126, 54, 97]. The analogy between ultracold
atoms in lattice potentials and electrons in crystals is manifestly a rich one. It extends
to strongly interacting ultracold Bose gases which exhibit both superfluid and insulating
behavior, and feature a gapped many-body energy spectrum. The tunability of interactions in optical lattices has led to the spectacular demonstration of these properties[64, 163].
Other interesting condensed matter systems that can be modeled in atomic gas
experiments are electronic devices. Ruschhaupt and Muga [140] have described an
atom device with diode-like behavior, and Micheli et. al. [111] have proposed a singleatom transistor that serves as a switch (see also [139, 141, 142, 37, 112]). Both of thesedevices depend on control and coherence at the single atom level. Moreover, Stickney
et. al. [162] have recently demonstrated that a Bose-Einstein condensate in a triple wellpotential can exhibit behavior similar to that of a field effect transistor.
Our intent is to establish ultracold atom analogs of electronic materials and semiconductor devices that can be used to leverage the vast body of electronic knowledge
and heuristic methods. From semiconductor materials and fundamental devices, the90analogy expands into what can be referred to as atomtronics. With diodes and transistors in hand, it is straightforward to conceive of atom amplifiers, oscillators, flip-flops,
logic gates, and a host of other atomtronic circuit analogs to electronic circuits. Such a
set of devices can serve as a toolbox for implementing and managing integrated circuits
containing atom optical elements [137, 15, 157] or quantum computation components[20, 21, 161] and might be of particular interest in the context of rapidly advancing
atom chip technologies [53, 132].
In this chapter, we introduce analogs of electronic materials, including metals, insulators, and semiconductors, in the context of ultracold strongly interacting
bosons [156]. We use lattice defects to achieve behavior similar to doped P-type andN-type semiconductors. The interest is to adjoin P-type and N-type lattices to create
diodes and then NPN or PNP structures to achieve behavior similar to that of bipolarjunction transistors. We show that such heterogeneous structures can indeed be madeto mimic their electronic counterparts.
Atoms in periodic structures and electrons in solid state crystals have much incommon. In both systems, particle motion occurs by tunneling through the potential
barrier separating two lattice sites. A particle, an electron or an atom, can delocalize
over the entire lattice and sustain currents. In both systems, currents are created when
there is a potential gradient which causes the particles to move from a region of higher
potential to a region of lower potential. In electronics, potentials arise from electricfields. In atomtronics, potential gradients can be understood in terms of chemical
potential gradients. The characteristics of both electronic and atomtronic devices are
examined by using a battery to apply a potential difference across the system andobserving the response in the current.
Different types of electronic conductors exist because electrons in a crystal structure occupy states of an energy spectrum that feature a band structure. The materials
can carry a current, i.e. are conducting, only if the highest occupied energy band is91
only partially filled with electrons. On the other hand, the system is an insulator if
all occupied bands are full. These properties of electronic materials can be directly
reproduced with weakly interacting fermionic atoms in periodic potentials.
Using fermionic atoms is not the only way in which analogs of electronic materials
can be created. In this chapter, we focus on ultracold strongly interacting bosonic atoms
in periodic structures. These systems can also be made to behave similarly to theirelectronic counterparts. Strong repulsive interactions prevent atoms from occupying
the same lattice site, mimicking the fermionic behavior of electrons. Hence, a current
can flow easily as long as there are empty sites available. However, once the filling
reaches one atom per site, the system becomes an insulator. A large energy gap givenby the repulsive onsite interaction must be overcome in order to add another particle
to this configuration. Particles added in excess of a filling of one atom per site canagain carry a current since they can move around freely above the filled layer of one
atom per site. The system remains a conductor until it arrives at a filling of two atoms
per site and becomes an insulator again. This is the same mechanism that drives the
Mott-insulator superfluid transition.
These properties show that a strongly interacting Bose gas features a close analog
of the band structure of electronic materials. The filling of the bands determines whether
a material is a conductor or an insulator. There is, however, an important difference
between atomtronic and electronic systems with respect to the band structure. The band
structure of strongly interacting bosons does not describe states that can be occupied
by a single particle independently of the configuration of others, but represents theenergies of many-body states. While the energy gap is due to the Pauli exclusion
principle combined with the single-particle band structure in the case of electrons, thegap arises from the onsite repulsive interaction between atoms in the bosonic case. In
this chapter, we will speak of many-body energy bands in the case of bosonic systems totake account of this important difference. The many-body band structure of strongly

92

Figure 4.1: Schematic of the atomtronic many-body energy band structure of stronglyinteracting bosons in a lattice. The first band consists of all states with a filling of lessthan one atom per site while the second band contains all states with filling between one
and two atoms per site. The two bands are separated by the onsite repulsive interaction
U.

interacting bosons in a lattice is depicted schematically in Fig. 4.1. The lowest band ismade up of states with between zero and one atom per site. The next band contains allstates with one to two atoms per site. Higher bands are formed analogously. The highest
occupied band of a conductor is only partially filled while insulators are characterized
by full bands.
Apart from the many-body character of the band structure, a second importantfeature of the atomtronic system is that the atomtronic conductor exhibits superfluid
rather than normal flow and in that aspect resembles an electronic superconductor.
In this chapter, the behavior of atomtronic materials in simple circuits is presented and we show how to use these materials to build more complex circuit devices,
specifically diodes and bipolar junction transistors. The chapter is structured as follows. Section 4.1 introduces the Bose-Hubbard formalism that will be used to describe
the atomtronic systems. In particular, it discusses the zero temperature phase diagram and the properties of an atomtronic battery. Doped atomtronic materials andthe current-voltage behavior of atomtronic wires will be discussed in Sec. 4.3. A diodeobtained by combining a P-type and an N-type atomtronic material with a voltage bias
applied by a battery is examined in Sec. 4.4. Section 4.5 presents the atomic analog

93
of a semiconductor bipolar junction transistor. Finally, Sec. 4.6 contains remarks on
possible applications, on the differences between atomtronic devices and their electronic
counterparts and on future perspectives. Appendix C contains details of the calculation
methods that were used throughout the chapter.

4.1

Quantum Phase Transition in the Bose-Hubbard Formalism
The zero temperature quantum phases of bosons in a lattice are key to under-

standing the different kinds of atomtronic materials. A system of repulsive bosons in a
one-dimensional chain of lattice sites can be modeled with the Bose-Hubbard Hamiltonian
Ĥ =

$ †
$
U$
n̂i (n̂i − 1) − Jâi âj +(/i − µ)n̂i ,
2 i
i
$ij%

(4.1)

where âi is the annihilation operator for a particle at site i, n̂i ≡ â†i âi is the number
operator at site i, U is the onsite repulsive interaction strength, J is the hopping matrix
element between nearest neighbors, $ij# labels nearest neighbors, /i is the externalpotential at site i and µ is the chemical potential of the system. The Bose-Hubbard
Hamiltonian is obtained by retaining only the contributions of the lowest single particle
Bloch band to the Hilbert space and by making a tight binding approximation (for areview see Ref. [179]). This yields an accurate description of an ultracold dilute Bosegas in a periodic potential at low energies. The zero temperature phase diagram of thisHamiltonian was first studied by Fisher, et. al., [52].
For very large onsite repulsion, U ( J, the system enters the regime of fermionization where bosons are impenetrable and only two Fock states, |ni # and |ni + 1#,
are needed at each site to accurately describe the system (two-state approximation).
Note that this is equivalent to mapping bosonic operators onto fermionic ones via the
Jordan-Wigner transformation [86, 143]. The data presented in this chapter is obtainedby considering the system described by the Hamiltonian, Eq. (4.1), in the fermionized

94

2
SF

µ/U

MI: 2

Band Width

1
MI: 1
0
0

Band Gap

0.1

J/U

0.2

0.3

Figure 4.2: Zero-temperature phase diagram of a Bose gas obtained from Eq. (4.2) for an
infinite one-dimensional lattice using the two-state approximation (fermionized regime,
U ( J). At large values of J/U the gas is superfluid (SF). Below a critical value of
J/U the system enters a Mott-insulator phase (MI) for integer filling, while it remains
superfluid for non-integer filling. The M I : 1 Mott-insulator region has one atom persite and the M I : 2 region has two atoms per site. The width of a Mott lobe at a given
J/U gives the size of the band gap in the many-body band structure while the widthof a band is given by the width of the superfluid region. The Mott lobe boundaries arelinear due to the use of the two-state approximation.

regime. We have verified that the error resulting from the two-state approximation
becomes negligible for U ≥ 100J (see also [12]).
The phase diagram of the Bose-Hubbard Hamiltonian contains information about
the many-body band structure of the system. At zero temperature, the Bose-Hubbard
model has two distinct phases, a Mott-insulating phase and a superfluid phase. Figure
4.2 presents the boundary between the conducting and insulating phases as a functionof the hopping parameter J/U and the chemical potential µ/U . The Mott-insulatingphase is entered below a critical value of J/U for an integer number of particles per site.
In this phase, strong interactions completely block particle motion rendering the gas
incompressible, that is ∂n/∂µ = 0, where n is the average filling of a site. The two lobes
shown limit the Mott-insulator zones with one atom per site (lower lobe, MI:1) and two

95
atoms per site (upper lobe, MI:2). The remainder of the depicted part of the phase
diagram is in the conducting superfluid phase, labeled SF . The superfluid phase is
obtained for non-integer filling. No insulating phase exists for values of J/U larger than∼ 1. Note that the triangular, non-rounded, shape of the Mott lobes in Fig. 4.2 is dueto the two-state approximation becoming increasingly inaccurate as J/U is increased.
Each value of µ/U and J/U maps onto a particular lattice filling, n. Figure 4.2
is obtained from the relation between these parameters in the superfluid phase for an
infinite number of lattices sites and strong interactions (J " U ) [123],µ = U (m − 1) + (−1)m 2mJ cos(πn) ,(4.2)

where m = 1, 2, . . . is the band index and the filling n satisfies (m − 1) ≤ n ≤ m. The
superfluid-insulator phase boundaries presented in Fig. 4.2 are obtained by setting n to
an integer and m = n for the lower boundary and m = n + 1 for the upper boundary.
The expression Eq. (4.2) is derived using a Jordan-Wigner transformation approachwithin the two-state approximation. A more accurate model of the phase boundary canbe calculated using a mean field model [98] to find that
µ
J
=
U
U

(

(2m − 1/2) − J/U ±&(J/U )2 − J/U (2m + 1) + 1/4

)

.

(4.3)

This is of a parabolic form and rounds the tips of the lobes.
Plotting the phase diagram as a function of µ/U rather than the number ofparticles is useful because in this way the gaps between many-body bands become
visible. The size of the gaps are given by the widths of the insulating zones, whereas
the sizes of the bands are given by the widths of the superfluid zones. Note that thesize of the gaps depend on the value of J/U . In order to have access to both insulating
and conducting phases, the ratio J/U must be small. If this ratio is too large there
is no well defined gap and therefore no transition to an insulating phase for integer
filling. The basic ideas presented in this chapter rely on this condition being satisfied.

96
They do not require the stronger condition U ( J for fermionization. The condition offermionization is assumed merely to facilitate calculations.
Most atomic systems do not have a constant density throughout the system due tothe ways in which the atoms are trapped. For instance, a harmonic trapping potential
will cause more atoms to be in the center then at the edges. In this case, a local
chemical potential, µi = µ − Vi , becomes important when determining the effects of
the lattice. This can lead to insulating regions and superfluid regions all in the same
system, traditionally in shell structures [64].

4.2

Atomtronic Battery
Just as for electronic circuits, the atomtronic version of a battery is crucial but it

is also subtle. In this section, we identify the basic properties and actions of the atom
analog of a battery.
Energy for electronic circuits is supplied by sources of electric potential. Furthermore, electric potentials are used to set the bias points of circuit elements to achieve
their desired behavior. For simplicity, we will use the term battery to refer to a device
that provides a fixed potential difference and can supply an electric current or atomic
flux.
In the electronic case, one is interested in the electric potential difference, or thevoltage, between two points. A specific potential difference between two points withina device or circuit is achieved by connecting those two points to the terminals of abattery. In an atomtronic circuit, the function of the battery is to hold the two contacts
at different values of the chemical potential, say µL on the left and µR on the right.
The applied voltage is then defined byV ≡ µL − µR .

(4.4)

The current flows from higher to lower chemical potential. Chemical potential difference

97
in atomtronic system is analogous to electric potential difference in electronic systems.
Note that the equilibrium value of the chemical potential is important since it sets theaverage filling of the lattice.
To understand the physics underlying this concept, note that bringing the systemwith chemical potential µ in contact with a battery pole of chemical potential µL > µleads to the injection of ∆n particles. The magnitude of ∆n is given by the differencein filling of states with chemical potentials µ and µL and can be determined from thephase diagram. The particle transfer increases with increasing µL within a superfluid
region and becomes constant as µL is moved into a Mott insulating zone where thesystem is incompressible. In the fermionization regime, the magnitude of ∆n is fixed by
Eq. (4.2). Analogous reasoning applies to the removal of particles at the battery polewith µR < µ. While µR remains in superfluid region of the phase diagram, changes
in the chemical potential will lead to changes in the atom transfer numbers. However,once µR enters the Mott insulator region of the phase diagram, small changes do notcreate changes in particle number transfer.
Feeding atoms into a circuit element through a contact at one end and removing
them through a contact at the other end generates a current. After time-averaging, thiscurrent reaches a constant value if the carrier excess at one end is replenished throughone contact with the battery at the same rate at which the deficit is maintained at theother end through the contact with the other pole of the battery. In addition to thistime averaged current, there is also the presence of the shot noise associated with thetime distributions of individual atoms entering or leaving the battery terminals.
Experimentally, a battery can be created by establishing two separate large systems which act as reservoirs, each with its own constant chemical potential. These may
also be lattices or other experimentally plausible systems, such as large harmonic traps,
containing a large number of atoms. Changing the frequency of the harmonic traps, oradjusting the lattice height, can be used to tune the chemical potentials of these reser-

98

Figure 4.3: Schematic of atoms in a lattice connected to an atomtronic battery. A
voltage is applied by connecting the system to two reservoirs, one of higher (left) and
one of lower (right) chemical potential. An excess (deficit) of atoms is generated in theleft (right) part of the system, giving rise to a current from left to right.

voirs. Each of these reservoirs can be connected to one end of the atomtronic systemand current is then possible from the higher chemical potential system to the lower one.
This configuration is displayed in Fig. 4.3. From a practical point of view, the chemical
potentials of the battery poles can be maintained by transferring atoms, possibly classically, between the two poles. In addition to the structure and chemical potential of
the reservoirs, the temperature of the reservoirs can play a role in the transfer of atoms
between the wire and the battery terminals. Since a high temperature reservoir canwash out the effects of the quantum phase transition in the lattice, a low temperature
reservoir is required. The temperature of the reservoir should be small enough so that
excitations on the order of U are highly unlikely.

4.3

Atomtronic Conductors
An attractive feature of atomtronic materials is that their conductivity properties

can be easily modified. The primary conductor for an atomtronic system is an opticallattice with no other external potentials. This corresponds to a wire in an electronicsystem. The properties of the wire are derived from the atomic filling and associatedchemical potential. For chemical potentials in the superfluid region of the phase diagram,the wire conducts well. However, when the chemical potential is in the Mott zones, the99

Figure 4.4: (a) Schematics of an N-doped lattice. The donor sites feature a level right
below the first empty many-body band. An atom which occupies this level can easily be
excited and move throughout the lattice. (b) Schematics of a P-doped lattice. Acceptor
sites have a level right above the highest full band. Atoms can easily be excited into
this level and allow for a hole to move throughout the lattice.

wire becomes an insulator. In this section, we first discuss the possibility of varying the
conductivity of a wire by lattice doping as is done in semiconductor systems. The idea
is that what is originally an insulator can be coaxed into behaving like a conductor justas in the electronic case. Secondly, we focus on the current-voltage characteristics ofthe different types of atomtronic materials.

4.3.1

Doped Materials

New materials with interesting properties can be designed by modifying the latticein which the atoms are confined, and hence the many-body band structure. This canbe done in various ways. In the case of an optical lattice, the periodicity of the lattice

100
can be modified by superimposing a periodic potential of different wavelength [124].
Disorder can be introduced by randomly modifying individual sites [78, 56, 149]. A
further handle on the properties of the material is the symmetry of a two-dimensionalor three-dimensional lattice [67, 147, 146]. Finally, it is conceivable to introduce a second
atomic species, either bosonic or fermionic, and to modify the properties of the material
via inter-species interactions.
The atom analog of the doping of a semiconductor is particularly interesting.
The aim of doping is to create energy levels in the energy gap between two bands.
N-type doping is associated with energy levels located close to the lower edge of thefirst empty band while P-type doping gives rise to levels close to the upper edge ofthe highest full band. Both are accomplished by modifying the potential at individual
lattice sites. N-type doping is achieved by replacing some lattice sites with donor sites.
These correspond to potential wells which are slightly deeper than those of the unmodified lattice. Analogously, P-type doping requires introducing acceptor sites of slightly
shallower potential. The two potential configurations are shown in Fig. 4.4.
The advantage of doping is that it can turn an insulator into a conductor without
having to excite particles across the band gap into the empty conduction band. In
atomtronics, this means that doping shifts the insulator zone boundary in the phase
diagram. Figure 4.5 compares the phase diagram of the undoped lattice with that ofa N-doped and a P-doped lattice. N-type doping shifts the insulating zone downwards
such that states that were previously in the insulating zone come to lie right above
the insulating zone where the lattice has a full valence band and a few free carriers in
the conduction band. Similarly, P-type doping shifts the insulating zone upwards suchthat states previously in the insulating region come to lie right below the insulatingzone where the lattice has an almost full valence band with a few free hole carriers. In
Fig. 4.5 the new insulating zones were derived for an infinite lattice where every thirdsite was doped by ∆/ = ±5J using a modified Jordan-Wigner transformation [123].

101

1
SF
P-type dopedµ/U

undopedMI: 100

N-type doped
0.1
J/U

0.2

Figure 4.5: Zero-temperature phase diagram of a Bose gas in an undoped (blue), an Ndoped (red) and a P-doped (green) lattice. The boundary of the lowest Mott insulator
lobe with a filling of one atom per site is displayed as calculated using the two-state
approximation. The results for the doped lattices were obtained by adding an energy of±5J to every third site in an infinite lattice using a modified Jordan-Wigner transformation. N-type doping turns an insulating state into a superfluid state with more than
one atom per site. Similarly, P-type doping turns an insulating state into a superfluid
state with less than one atom per site.

4.3.2

Material Current Properties

This section discusses the current response of atomtronic wires to an applied
voltage. The magnitude of the current depends on the properties of the material and on
the nature of the contact with the battery. The material can be in the Mott insulatingphase or in the superfluid phase. In the Mott phase, a small voltage does not yielda current. Only at a large voltage, of the order of the gap U , is the battery able togenerate a current by feeding particles into the next unoccupied band. If the material
is in the superfluid phase, a small voltage is enough to generate a current. Due to thesuperfluid nature of the atomic carriers, this current is not slowed by friction. Hence, theratio of voltage to current does not have the physical meaning usually associated withresistance, but reflects the limits on the current at a given voltage due to factors other102
than dissipation. A primary limit is set by the hopping parameter J which quantifies
how fast atoms can move from one site to the next. The value of the hopping parameter
depends on the shape and the depth of the lattice. A further limit on the current is dueto the interaction between atoms. Repulsion forces them to move in a highly correlated
fashion. For this reason, the current does not grow linearly with the number of carriers
at fixed voltage. Instead, the current per particle drops as more carriers are added andeventually goes to zero when the filling is one atom per site and the system becomes aninsulator.
Apart from material dependent factors, the nature of the contact with the battery
can set the maximum achievable current. If it is more difficult for atoms to pass throughthe contact than to hop from one lattice site to the next then the current is not limited
by J but by the rate at which the battery can feed in and remove particles. Thissituation is encountered when operating a battery in a regime of very weak coupling.
The same properties are present in electronic systems where there can be different typesof contacts, such as rectifying and ohmic contacts [8]. In this chapter we focus on theregime of maximum currents where currents are limited by the hopping parameter J
and are not influenced by the properties of the battery contact.
Our calculation of the current response is performed in several steps. First, thechemical potentials of the left and right batteries are chosen. These chemical potentials
correspond to an average site filling. For a wire or a diode, a two well system is analyzed
corresponding to a junction in the system. An initial state of NL atoms on the left and
NR atoms on the right, as determined by the chemical potentials and Eq. (4.2), is set
with a phase corresponding to the lowest energy state with this density configuration.
The state is propagated in time. The calculated current is determined by the time average current during half of an oscillation. This calculation is motivated by properties oftraveling waves in which many of these oscillations occur in sequence. See the Appendix
for a further discussion of the calculation of currents.

103
Figure 4.6 presents the current as a function of voltage for lattices, or wires,of average filling n = 1.1, 1.3, 1.5, 1.7, 1.9. The calculation performed to obtain thesecurrent-voltage characteristic curves are described in the Appendix. At fixed filling,
the current increases monotonically with increasing voltage. The maximum current
attainable for a half-filled second band is ∼ 1.4 J/h̄. The voltage in Fig. 4.6 is given
in units of ∆µmax . This quantity denotes the chemical potential differences that yield
the maximum currents. It corresponds to the chemical potential difference at which thesystem, for a given filling, enters an insulating regime at one of the two battery contacts.
The curves for fillings with an equal number of free particles and holes coincide,demonstrating the equivalence of hole and particle motion for U ( J. The inset in
Fig. 4.6 presents the currents at different average fillings, corresponding to different
materials, at ∆µ = ∆µmax . The plot is symmetric around its maximum at half fillingreflecting particle-hole symmetry. Note that the average current per free particle, i.e.
holes at n > 1/2, is constant.

4.4

Atomtronic Diode
A diode is a circuit element that features a highly asymmetric current-voltage

curve. It allows a large current to pass in one direction, but not in the other. The analog
of an electronic diode is an atomtronic circuit element that lets an atomic current passthrough when applying a voltage, V = µL − µR , while allowing no current or only asmall saturation current for a voltage, V = −(µL − µR ).
In solid state electronics, diodes are built by setting up a PN-junction in which a
P-type semiconductor is brought into contact with an N-type semiconductor. Electrons
move through the junction until an equilibrium is reached. This process depletes thejunction region of free charges, leaving behind the static charges of the donor andacceptor impurities. As illustrated in Fig. 4.7(a), this creates an effective potential
step across the junction. When a reverse bias voltage is applied, the energy barrier is

104

2

Max Current

Current (Units of J/h̄)

2

1

1

0
0

0.5
1
Average Filling

n=1.5
n=1.3, 1.700

n=1.1, 1.9
0.5
1Voltage (Units of ∆µmax)

Figure 4.6: Current as a function of chemical potential difference (voltage) for different
materials (different fillings n) in the fermionization regime where bosons are impenetrable. A wire with a certain number of atoms carries the same current as a wire with that
number of holes (particle-hole symmetry). The chemical potential difference is given
in units of ∆µmax . This quantity denotes the chemical potential differences that yield
the maximum currents and corresponds to the chemical potential difference at which
the system enters an insulating regime at one of the two battery contacts for a given
n. Inset: Current as a function of filling at ∆µmax . Maximum currents are attained at
half filling. Note the symmetry in the current due to the particle-hole symmetry.

increased, reducing the flux of electrons from N-type to P-type. At the same time thenumber of electrons that can fall down the step remains constant giving rise to a reverse
bias saturation current that is independent of voltage (see Fig. 4.7 (b)). However, if thediode is forward biased, more electrons are able to move from the N-type to the P-typematerial than at equilibrium since the potential barrier is decreased by the forward bias
voltage (see Fig. 4.7(c)). A detailed discussion of the diode behavior of a semiconductor
PN-junction can for example be found in [8].

4.4.1

Diode PN-junction configuration

As discussed previously, we can design atomtronic wires whose conductivity properties can be described by locating the material’s chemical potential in its phase diagram.

105

Figure 4.7: Schematic of the conduction band of an electronic (left panels) and atomtronic (right panels) PN-junction diode in (a) equilibrium, (b) reverse bias and (c)forward bias with N-type materials on the left and P-type on the right of each configuration. Left panels: The electronic system features a voltage dependent energy barrierat the junction leading to an increasing (decreasing) flux from the N-type to the P-typematerial as the junction is forward (reverse) biased while the current from P to N is
independent of voltage. Right panels: The operation of an atomtronic diode is based on
the existence of insulating phases where ∂n/∂µ = 0, i.e. a change in voltage does notlead to a change in particle transfer between battery and system. As a consequence,
only a small current can flow from N-component to P-component through an atomtronic diode in reverse bias while in forward bias the particle transfer between battery
and system can be varied over a large range.

Composite materials, such as PN-junctions, can be produced by connecting lattices ofdifferent doping. Another possibility of building a junction is to superimpose additional
external potentials, for example a simple potential step. This has the effect of shifting
the phase diagram of a part of the lattice upwards or downwards with respect to that of
the rest. Since phase boundary effects are small, the state of the different components
can be accurately described by the phase diagrams of the individual materials (see Ap-

106pendix for details). This means that a local chemical potential can be associated witheach component of the conductor that can be located in each material’s phase diagram.
Thereby, the local conductivity properties can be identified. Of course, at equilibrium,
i.e. at zero voltage, the composite material is actually described by a single chemical potential µ but it is µ relative to the zero point energy of each lattice site that determines
the filling of each site.
In the following, we focus on a junction created by an external potential step since
this may be the easiest experimental realization of a diode. The potential step could be
generated experimentally by exposing one part of the system to off-resonant laser light.
The main characteristics of the diode behavior are not affected by this choice.
The conduction band of a semiconductor PN-junction features a small thermal
electron population on the P-side and a considerably larger electron filling on the Nside. An example of an atomtronic equilibrium configuration with analogous features isrepresented by the squares in Fig. 4.8. The small (large) population of the second bandin the P-type (N-type) material yields the analog of the small (large) thermal electron
population of the conduction band in a semiconductor.

4.4.2

Diode current-voltage characteristics

To achieve diode characteristics, we exploit the possibility of undergoing a quantum phase transition between insulating and superfluid phases. The materials are configured such that the chemical potentials of the battery poles remain in the superfluid
regime when hooking up the battery in one direction, but they easily enter insulating
regimes when the voltage is applied in the opposite direction.
The effect of applying a voltage is illustrated in Fig. 4.8. Forward bias is achievedby connecting the P-side to the low voltage pole and the N-side to the high voltage poleof the battery. In this situation, the chemical potentials of the battery poles are locatedin the superfluid region of the phase diagram and atoms can flow from the P-component

Reverse Bias

Forward Bias

(a)

MI: 2

(b)

MI: 2

N-typeSF

1

N-typeSF

1P-typeµ/U

P-type

MI: 1

0
0

0.1
J/U

107

MI: 1

0.2

0
0

0.1
J/U

0.2

Figure 4.8: Phase diagram of the PN-junction configuration of an atomtronic diode. The
small (large) population of the second many-body band in the P-type (N-type) material
yields the analog of the small (large) thermal electron population of the conduction bandin a semiconductor. The states of the P-type and N-type materials at zero voltage arerepresented by the squares. Arrows indicate the chemical potential difference (voltage)
imposed to obtain a reverse (left) and a forward biased junction (right).

to the N-component. The larger the applied voltage, the larger the generated current.
When the battery contacts are switched, the diode is reverse biased. As the voltage isincreased a small current starts flowing. However, as soon as the voltage is large enoughto make the battery chemical potentials enter the insulating zones, the current can notincrease any further. As a consequence, the current-voltage curve is asymmetric. The
remnant current obtained in reverse bias, the saturation current, becomes smaller as
the components’ initial states are moved closer to the insulating phase.
Figure 4.7 presents a schematic comparison of the conduction band of an electronic
and an atomtronic diode obtained using a step potential. For the atomtronic diode, thefact that only a small current can flow in reverse bias is not due to the presence of a
voltage-dependent energy barrier at the junction as in the electronic case. Instead, it
arises from the battery chemical potentials moving into insulating zones corresponding
to a full conduction band on the N-side and an empty conduction band on the P-side.

108
An important difference between electronic and atomtronic case is the opposite directionof current flow. In forward bias, atoms flow from the P-type to the N-type material as
opposed to the other way around for electrons in a semiconductor.
Figure 4.9 displays the highly asymmetric current-voltage curve obtained from ourcalculation. The potential step is chosen such as to yield an equilibrium configuration
with a filling of 1.99 and 1.01 atoms on the N-side and P-side respectively. In thisconfiguration, we obtain a reverse saturation current of 0.14 J/h̄ while in forward bias,
currents can exceed 1.4 J/h̄. Note that the current changes strongly in the vicinity of
V = 0. This is reminiscent of the behavior of an electronic diode as the temperatureapproaches zero. Reducing the potential step lowers the population difference between
P-type and N-type materials at equilibrium and leads to an increase in saturation current
and to a decrease in the slope of the current-voltage curve around V = 0.
The diode currents have been calculated in the same manner as the currents
carried by atomtronic wires except for the addition of the potential energy step between
the two sites (see Appendix for details of calculations).

4.5

Atomtronic Transistor
As in electronics, the highly asymmetric current-voltage curve of atomtronic

diodes can be exploited to build a transistor. Bipolar junction transistors (BJT) arecircuit elements that can serve as amplifiers and switches. In electronics, they consist
either of a thin P-type layer sandwiched between two N-type components (NPN) or athin N-type layer between two P-type components (PNP). For our discussion, we consider an NPN configuration. A detailed discussion of semiconductor bipolar junction
transistors can, for example, be found in [8]. The basic circuit schematic is displayedin Fig. 4.10. The voltage VBC that is applied to the PN-junction formed by the middle
component, the base, and one of the outer components, the collector, puts this junction
into reverse bias. At the same time, the other junction formed by the base and the109

Current (Units of J/h̄)

1.5
1
0.5
0

−0.5
−4

0

Voltage (Units of J)

4

Figure 4.9: The characteristic current-voltage curve for an atomtronic diode. The larger
the forward bias (voltage > 0), the higher the current of atoms flowing from the Pcomponent and to the N-component. In reverse bias (voltage < 0) the current saturates
since the particle transfer between battery and system can not be increased beyond acertain small value when the battery pole chemical potentials enter insulating zones
where ∂n/∂µ = 0.

other outer component, the emitter, is put into forward bias by applying a voltage VEB .
The key idea is to use the voltage VEB to control the current IC leaving the collector
element and thereby achieve differential gain in IC relative to the base current IB . At
VEB = 0 one is simply dealing with the reverse biased base-collector junction. In thiscase, the currents IC and IB both equal the small reverse bias saturation current of
the base-collector junction. The collector current IC grows drastically when VEB is
increased such that the emitter-base junction is forward biased. This effect relies on thebase region being very thin. The forward bias gives rise to a flow of electrons from theemitter into the base region, thereby significantly increasing the number of electrons at
the base-collector junction. Recall that at VEB = 0, the base is depleted of electronsby the reverse bias VBC . The emitter-base junction thus serves to greatly modify the
number of carriers in the base that are subjected to the base-collector reverse bias. This
leads to an increase of IC beyond the saturation current. Since the base is extremely

110

Figure 4.10: Circuit schematic of an electronic or atomtronic bipolar junction transistorof the NPN-type. A thin P-type component (base) is sandwiched between two N-type
components (emitter and collector). The key feature is that the voltage VEB can beused to obtain gain in the collector current IC relative to the base current IB .

thin there is less opportunity for the electrons to leave the base compared to entering
the collector. Because of this, most of the current that enters the base from the emitter
moves on to the collector instead of it leaving out the base terminal. Therefore therelative changes in the current from the base IB and from the collector IC yields a large
differential gain dIC /dIB .
These key features of an electronic transistor can be translated into atomtronics
with atoms taking the place of electrons. The important point is to mimic the carrier
densities and thus the fillings of the three components. This can be achieved by dopingor, equivalently, by using potential steps. As with the diode, we focus on the latter
implementation. The atomtronic transistor can be created by setting up a configurationsuch that for zero voltage the conduction bands of the left and right regions have a large
filling compared to the sandwiched thin base region. The phase diagram for such anarrangement is depicted in Fig. 4.11(a), where the squares represent the equilibrium
configuration, while the arrows indicate the way the battery chemical potentials are

111
MI: 2

(a)

Collector

1

Emitter

µ/U

Base

SF

MI: 1

0
0

0.1

0.2

0.3

J/U
200

Gain |dIC /dIB |

100

50

0
−4

0.2

0
−4

0.2

0
−2
0
VEB (Units of J)−2
VEB (Units of J)

IB (Units of 10−2 J/h̄)150

IC (Units of J/h̄)

( b)

0

Figure 4.11: (a) The squares indicate the equilibrium configuration of an atomtronic
bipolar junction transistor in the phase diagram while the arrows show the direction inwhich the chemical potentials of the battery contacts are varied. (b) Differential gaindIC /dIB as a function of emitter-base voltage VEB . Inset: Base current and collector
current as a function of emitter base voltage. The voltage is changed by varying theemitter contact potential from µ(nE = 1) to µ(nE = 1.5). The base contact is kept atµ(nB = 1) right below its equilibrium value µ(nB = 1.005) while the collector contact
is kept at the equilibrium value µ(nC = 1.5).

tuned. A small voltage is applied to the collector-base junction such that a small current
flows from the collector into the base. When the emitter battery chemical potential µE
is lowered, atoms move from base to emitter and leave through the emitter, giving rise
to a non-zero emitter current IE . As a back-effect, this leads to an increase in IC since
the fast removal of atoms from the base through the emitter allows more atoms to moveinto the base from the collector.

112
Meanwhile, the base current IB becomes smaller the further the emitter chemical
potential µE is lowered. This is due to the base being very thin relative to the emitterand the collector, so atoms traverse preferentially from collector to emitter rather thanleaving out the base contributing to IB . The effect of the forward bias VEB on the basecurrent IB is thus opposite to its effect on IC . Therefore, our atomtronic transistorfeatures an inverted amplification in which a small decrease in the base current goesalong with a large increase in the collector current (negative gain).
An example for the behavior of the currents in an atomtronic transistor is given
in Fig. 4.11(b). The equilibrium configuration has 1.5 atoms per site in emitter andcollector and 1.005 atoms per site in the base. In Fig. 4.11(b) we plot data obtained
for the individual currents IC and IB upon variation of the emitter chemical potentialµE while keeping the base battery contact at a chemical potential µ(nB = 1) in the
n = 1 Mott-insulating zone and the collector contact at µ(nC = 1.5). To demonstratedifferential gain we display the quantity |dIC /dIB | as a function of VEB . The details of
the calculation are presented in the Appendix.

4.6

Atomtronics Conclusions
We have shown how strongly interacting ultracold bosonic gases in periodic po-tentials can be used as conductors and insulators in a circuit and how they can be
employed to build atomtronic analogs of diodes and bipolar junction transistors. From
here, the implementation of an atom amplifier is immediate. An atom amplifier is a
device that allows control of a large atomic current with a small one. The transistor
presented above directly serves this purpose since small changes in the base current
bring along large changes in the collector current. It is straight-forward to conceiveof more complex devices such as a flip flop, a bistable device that uses cross negative
feedback between two transistors.
The similarity in qualitative behavior goes along with a number of significant

113
differences in the underlying physics. Firstly, in the atomic case the energy gap resultsfrom interactions rather than from statistics as in electronics. Secondly, the atomic
currents are superfluid. As a consequence, the ratio between voltage and current has
the meaning of a dissipationless resistance. Further differences arise in both diodes
and transistors. Our atomtronic diode does not feature a depletion layer, i.e. it does
not exhibit a voltage-dependent energy barrier at the junction. The asymmetry in thecurrent-voltage curve results from voltage-sign dependent quantum phase transitions toan insulating phase. As a consequence, atoms flow from P-type to N-type in forward bias
rather than flowing, as in electronics, from N-type to P-type. Note that this difference
can not be resolved by drawing the analogy between atom holes and electrons rather
than atoms and electrons since this would also require relabeling the N-type material
as P-type and vice versa and hence the current direction would again be reversed incomparison to the electronic case. Due to the difference in diode behavior the atomic
collector current in a transistor flows from collector into base and the emitter current
flows from base into emitter, i.e. opposite to electronic flow in a NPN transistor. A
significant difference in the qualitative behavior of electronic and atomtronic transistor isgiven by the gain being negative in the atomtronic case. The collector current increases
as base current decreases. We expect that this does not affect the functionality ofdevices based on the operation of bipolar junction transistors. Yet, an adaptation of
their design will be necessary.
The data presented in this chapter is obtained from calculations for a one-dimensionallattice. This choice is of an entirely practical nature. The basic ideas also hold for twoand three-dimensional cubic lattices and extend to other lattice geometries that maketransitions between superfluid and insulating phases upon changes of the chemical potential.
Working with higher-dimensional lattices increases the magnitude of the currents
that can be achieved. In an one-dimensional optical lattice with an experimentally

114
realizable depth of 10 ER the hopping parameter J ∼ 0.02ER , where ER = h̄2 π 2 /2md2
is the recoil energy which is fixed by the lattice period d. For 87Rb and a lattice periodof d = 400nm, the recoil energy is ER = (2πh̄)3.55kHz yielding currents on the scale
J/h̄ ∼ 2π71Hz. This implies, that the currents through a three-dimensional latticeof the same period and depth with a cross-section of 10 × 10 sites are of the order
100 × 2π71Hz.
Experiments aimed at atomtronic devices will be challenging given the current status of ultracold atom technology. Among other challenges is the need to combine atomic
potentials that have typically been demonstrated in isolated systems, and, moreover,is the need to spatially cascade regions in which atoms undergo nonlinear interactions.
An optical lattice, for example, might well serve as the basis for fabricating atomtronic
P-type and N-type materials. Alternatively, lattice potentials could be created by passing current through micro-fabricated wires on a chip surface. And perhaps coupling an
atom magnetic waveguide to an optical or non-optical lattice as a means of transporting
atoms to and from the lattice can set the stage for atom diode experiments. In principle,
atom chips provide the means for supporting different types of atomic potentials on asingle substrate. The lattice systems used in atomtronic circuits do not necessarily haveto support coherent currents. The important feature is that the system’s tendency toaccept or not accept particle transfer from the battery exhibits the incompressibility
associated with the Mott insulating phase. We have provided a conceptual frameworkof semiconductor material and device analogs that can serve as building blocks to moresophisticated atomtronic devices and circuits. A discussion of actual atomtronic device
implementation necessarily goes beyond the scope of this work. Pepino, et. al., havefurther developed methods for constructing devices from multiple transistors, such as
quantum logic gates [125].
Finally, it is important to keep in mind that this chapter develops atomtronics
within the Bose-Hubbard model. This model provides an excellent description of ongo-

115
ing experiments with ultracold bosonic atoms in optical lattices. Other Hamiltonians
might offer alternative ways of drawing the analogy with electronics. A natural choicefor further study are Hamiltonians describing bosons with beyond onsite interactionsand Hamiltonians for fermionic gases.

Chapter 5

Photons Interacting in a Cavity

The interactions between light and matter can be examined from several perspectives. The previous chapters showed how the properties of atoms can be altered by thepresence of an electromagnetic field. The intensity of the electromagnetic field, in theform of an optical lattice, is able to enhance the effects of atomic two-body interactions
relative to the kinetic energy of the atoms. In this chapter, the interactions between
light and atoms are examined from the perspective of the properties of the photons.
In the presence of an atomic field, photons can exhibit many interesting properties,
the most famous of which is lasing. Effective photon-photon interactions mediated byatoms can also arise. This chapter presents the Hamiltonian that is used to describe theproperties of a collection of photons confined to a cavity filled with a nonlinear medium.
The nonlinear medium could be created in several ways. For instance, a Rydberg gas of
atoms or the deposition of a nonlinear crystal can be used. For the description of thissystem, it is useful to consider the medium as a collection of two-level atoms.
In this chapter, a model of a photon fluid confined to a nonlinear Fabry-Perotcavity will be presented. First, the bare modes of the cavity are briefly discussed. Thenthe Hamiltonian of the photon fluid is derived by examining the nonlinear polarizabilityof the nonlinear atomic medium. Several properties of a photon fluid are then examined.

117

5.1

Photon Modes of the Cavity
In deriving the Hamiltonian that describes the evolution of a collection of pho-tons confined to a cavity interacting via atomic coupling, it is important to accurately
describe the bare modes of the cavity. These photon modes which are associated witha vacuum cavity with perfectly conducting walls will be used as the basis states of thesystem. The system is then an electromagnetic field in a cavity, which we may take to bea box of lengths lx , ly , and lz . The walls of the cavity will be considered perfect conductors and the interior of the cavity will have no electric polarizability or magnetization.
Maxwell’s equations and the boundary conditions completely define the system.
Each mode of the cavity is identified by three quantum numbers, ni , one for eachdirection, and a polarization direction. As long as none of the n’s are zero, there are
two independent polarizations for each set of quantum numbers. If an n is zero , thenthere is only one possible polarization. Note that if two or more of the n’s are zero,then the field vanishes. If ni is zero, then the polarization is in the i direction.

5.2

Quantum Hamiltonian
Here we construct the quantum Hamiltonian for an electromagnetic field inter-acting with a nonlinear medium. The energy density of an electromagnetic field in a
medium is given byU=+
1 *# #
# ·H
# ,
E·D+B
2

(5.1)

# ≡ /0 E# + P# and H# ≡ B/µ# 0−M# , P# is the polarization, and M# is the magnetiwhere D
zation. If the the medium has zero magnetization, the Hamiltonian is given by
H=

1
2

# *
V

+

# ·E
# +B
# · B/µ# 0+E# · P#nl ,
/0 E

(5.2)

# is the nonlinear part of the polarization. It will be this polarization
where P#nl ≡ P# −/0 E
that must be calculated by examining the interactions of the photons and the medium in

118
the cavity. It is then straight forward to use traditional second quantization techniqueson the photon modes to obtain the Hamiltonian
Ĥ =

$(

(

h̄ω â†i âi +

i

1
2

)

1
4

+

# *
V

Êi · P̂nl + P̂nl · Êi

+)

,

(5.3)

where
# i = âi E# − e−iωi t + â† E# + iωi t ,
E
i
i i e

(5.4)

is the electric field operator of mode i with energy h̄ωi , âi is the annihilation operator#− = E# +∗ contains the spatial dependence of the mode.
for a photon in mode i, and E
i
i
5.2.1

Polarization of Atoms

The dipole moment characterizes the polarization of the medium. We will therefore calculate the dipole moment of an atoms that creates the nonlinearity polarization.
It is useful to return momentarily to a classical description of the multi-mode light field,
# =
E

$

# i cos(ωi t) .
E

(5.5)

i

where mode i has energy h̄ωi , and later we will quantize the photons,# i → âi E# i + â† E
#∗
E
i i ,

(5.6)

The polarization of a two level atom is given byP#≡ T r(er̂ρ̂)

(5.7)

$$

(5.8)

= e

ij

C

$j|C|i#ρij Ĉ ,where C ∈ {x, y, z} represents each direction and ρ̂ is the density matrix decribing theatomic degrees of freedom. This formula simplifies since only states of differing parity
are connected. We are therefore left with,
P#

= e

$$

$g|C|e#ρeg Ĉ + c.c.

(5.9)

ρeg

$

e$g|C|e#Ĉ + c.c.

(5.10)+ c.c. ,(5.11)

C

=

$
e

=

$
e

e

C
#e∗

ρeg d

119
where d#e ≡ $e|e#r |g# is the dipole moment between ground state g and excited state e.
It is now necessary to calculate the components of the density matrix. The
Hamiltonian describing the interaction between the atom and the photons, without
spontaneous emission, is given byĤ = h̄ωe$e|e#$e| + h̄ωg |g#$g| +$e

Veg |e#$g| +$e

Vge |g#$e|(5.12)

where ωe is the energy of the excited state, ωg is the energy of the ground state, and#.
Veg = −d#e · E(5.13)The dipole approximation can be used since the wavelength of the electromagnetic field
is much larger than the atomic size and so the electric field can be assumed to havethe value at the center of the atom. Note that this is the same Hamiltonian that wasused to describe the creation of an optical lattice earlier in this thesis and the same
approximations used still hold such as the use of the dipole approximation and therotating-wave approximation.
The density matrix equation of motion in which spontaneous emission from theexcited states to the ground state is allowed is given by
3 Γ $2
3
i2
i i
i i
i
i
ρ̂˙ = − Ĥ, ρ̂ −σ̂+σ̂− ρ̂ + ρ̂σ̂+σ̂− − 2σ̂−
ρ̂σ̂+
,
h̄
2 i

(5.14)

i is the raising operator from the groundwhere Γ is the spontaneous emission rate, σ̂+state g to the excited state i and σ̂− is the lowering operator. The effects of spontaneousemission is found by following a traditional master equation formalism and ending withterms of the Lindbladt form as seen in Eq. (5.14).
The density matrix equation of motion, Eq. (5.14), can be rewritten in termsof the individual components. The matrix elements can be broken up into three sets,
whose equations are given by
ρ̇eg





$
i
Γ
+ iΩ ρeg − Veg ρgg −Vf g ρef  ,
= −
2
h̄
f
(

)

(5.15)

120
ρ̇ef
ρ̇gg

i
(Veg ρgf − Vgf ρeg ) ,
h̄ 8
(
9)
$
i
=
Γρee + Veg ρge + c.c.
,h̄e= −Γρef −(5.16)(5.17)

where Ω ≡ ωe − ωg is the atomic transition frequency and both e and f refer to excitedstates.
Just as the electric field is composed of a set of frequency dependent terms, theelements of the density matrix can be expanded in a Fourier series,
ρeg =

$

ρneg e−iωn t ,(5.18)nfor example. Since the interaction terms Veg are linear in the electric field, expanding
in a Fourier series removes the temporal components of the density matrix equations of
motion. It is then possible to write the density matrix elements in a series expansion of
the electric field modes, e.g.,
ρneg = ρn(0)eg +$

ijf1 f2*# i · d#f1
ρn(ij)
Eeg+*+# j · d#f2 + · · · .
E(5.19)

As long as the population of the atomic excited state is small, terms beyond second
order can be neglected. This corresponds to an electric field that is far-detuned fromthe atomic transition. The elements of the density matrix can be solved for and replacedback into the equation for the polarization, Eq. (5.11). Under these assumptions, thepolarization of the atoms up to the first two nonvanishing orders of electric field is given
by
P#

= in

$

D(ωc )

i

#2
|d|
#i
e−iωi t E
2h̄

(5.20)−in$D(ωk + ωj − ωi )# 4 D ∗ (ωi ) + D(ωj )|d|# i∗ · E# j )E#ke−i(ωk +ωj −ωi )t (E8h̄3 Γ − i(ωj − ωi )−in$

D(ωk + ωj − ωi )# 4 D ∗ (ωi ) + D(ωk )|d|#∗ · E# j )E#ke−i(ωk +ωj −ωi )t (Ei8h̄3 Γ − i(ωk − ωi )

ijk

ijk

+c.c. ,

121
where n is the atomic density,
1,
2 + i(Ω − ω)D(ω) ≡ Γ(5.21)

is a detuning factor, and Ω is the atomic transition frequency. Note that this is still aspatially dependent polarization.

5.2.2

Fourth Order Hamiltonian

It is now possible to create the Hamiltonian describing the atom-mediated photonphoton interactions. In the rotating-wave approximation, only energy conserving terms
remain. By replacing the polarization from Eq. (5.21) into the original Hamiltonian ofEq. (5.3) and using the quantized form of the light from Eq. (5.6), a new Hamiltonian
up to the two lowest orders of the electric field is found. The Hamiltonian is then given
by
H =

$

h̄ωi (1 + ∆i ) â†i âi +i$

Γijkl â†i â†j âk âl(5.22)

ijkl

where the Hamiltonian has been written in proper normal ordering. The factor ∆i is anadjustment of the mode energies due to the linear polarizability of the medium and canbe approximated by ∆i ≈ ∆0 and incorporated into the bare mode frequencies since(ωi − ω0 )/Γ " 1. The process of a single photon temporarily interacting with an atomand quickly leaving creates this polarization and is dependent on the parameters of themedium. For the fourth order term, we need only consider the zeroth order contribution
to Γijkl , since again (ωi − ω0 )/Γ " 1. This term becomesΓijkl ≈ −n#4
L20 Ω0 ω02 |d|
16/2 h̄#V#i · E# k∗ )(E# j∗ · E# l) ,(E(5.23)

where
L0 ≡ Γ2
4

1+ (Ω − ω0 )2,(5.24)

122and depends on the overlap integral of photons wave functions. Note that the approximations to the coupling terms do not have to be explicitly made and can for a generalsystem be frequency dependent.
The presence of the overlap integrals in Eq. (5.23) leads to several properties ofthe system since they describe the rates at which photons can interact, if at all. Only
if the addition or subtraction of the quantum numbers for each mode can be set equal
to zero in each direction is an interaction allowed. This can be understood from theperspective that a standing wave can be thought of as two counter-propagating waves
with opposite momentum. Therefore, even though the full mode has zero momentum
since it is trapped in a cavity, one of the waves can still interact with a wave fromanother mode. Only if this collision satisfies momentum conservation is the interaction
allowed. This is how momentum conservation plays a role in the photonic collisions.
The quantum number can also be thought of being related to the standard deviation of
the momentum of the mode.
There are three types of interactions that occur. The first kind are nonlinear
self-interaction terms of the form,V ∝ â†i â†i âi âi .

(5.25)

These interactions create a mean-field energy shift in the photon modes. Interactionsof the form,V ∝ â†i â†j âj âi ,(5.26)

conserve the number of photons in each mode and again lead to a mean field energy shift
but are due to the interactions of two different modes. Self interactions are a subset ofthese interactions. Interactions of the form,V ∝ â†i â†j âk âl ,(5.27)

actually change the photon mode populations, and are the only ones to do so. Thesetypes of interactions are those associated with effects such as four wave mixing.

123

5.3

Photon Fluid
A weakly interacting photon gas has become known as a photon fluid [33]. This

fluid possesses the same Bogoliubov dispersion relation for elementary excitations as thatof a weakly interacting Bose gas of atoms. In addition, the photonic system possesses
a well defined chemical potential, which noninteracting photons lack. These properties
imply that the photon fluid system admits a superfluid state of light.
The Bogoliubov excitations can easily be seen by examining the Hamiltonian of
the photon fluid, which can be expressed as,Ĥ =$h̄ω(p)â†p âp +p$

V (k)â†p+k â†q−k âp âq ,(5.28)

kpq

where
&

ω(k) = c kx2 + ky2 + kz2 .

(5.29)

This describes the scattering of photons with momentum transfer k. The free energy ofa photon is then given by&

E = c p2⊥ + m2 c2(5.30)&

where p2⊥ ≡ h̄ kx2 + ky2 and mc ≡ h̄kz = h̄nπ/L. If p2⊥ " p2z , then the energy can beapproximated by a Taylor expansion to yield,
E ≈ mc2 +

p2⊥
.
2m

(5.31)The zero-point energy can be offset to remove the constant mc2 factor. For a large
number of photons in the ground state mode, N0 ( Ni , a mean-field approximation√
can be made in which â0 = N0 . In this case, the Hamiltonian becomes,Ĥ =$/# (p)â†p âp + N0

$

p>0

p>0

*

+

V (p) â†+p â†−p + â+p â−p + V (0)N02(5.32)

where#/ (p) ≡!

"

p2⊥
+ N0 V (p) + N0 V (0)2m,(5.33)

124
and collisions of two high energy photons has been neglected since their populations arevery small, N02 ( N0 Ni .
It is possible to diagonalize this quadratic form of the Hamiltonian using traditional Bogoliubov methods. With the canonical transformation,α̂p = up âp + vp â†−p(5.34)α̂†p = up â†p + vp â−p(5.35)

in which u2p − vp2 = 1, the Hamiltonian becomesĤ = constant +$(ω # (p) α̂†p α̂p +

p

1
2

)

(5.36)withh̄ω # (p) ='p4
p2 N0 V (p)
+
.
m
4m2

(5.37)

This is the same form as the traditional Bogoliubov dispersion relation. When the momentum is small, the excitation spectrum is linear and, therefore, excitations behaveas phonons. For larger momenta, the quasi-particles have a quadratic dispersion relation and behave as massive, nonrelativistic free particles. This transition occurs at a
momentum of approximately&

pc = 2 mN0 V (pc ) .

(5.38)

Figure 5.1 presents the excitation spectrum with the assumption of a constant interaction strength, V (p) = V (0). For small momenta, the spectrum is given byh̄ω = p'NV,m(5.39)

and for large momenta the spectrum becomes quadratic
h̄ω =

p2
+ NV .
2m(5.40)

Large momentum excitations are essentially massive free particles with a mean-fieldinteraction energy offset.

1255h̄ω (p2c /2m)

Bogoliubov
Quadratic
Linear001p(pc )

2

Figure 5.1: The solid line is the dispersion relation for excitations of the ground stateof a photon fluid. The dashed curve is the quadratic free-particle spectrum associatedwith high energy excitations. The dash-dot curve is the characteristic linear componentof the spectrum which describes phonon-like excitations.

The Bogoliubov dispersion relation implies that the photon fluid should exhibit
the property of superfluidity in which dissipationless flow occurs. The transition between
dissipation less flow and turbulent flow should occur at a well defined critical velocity.
In addition, the occurrence of quantized vortices would be a sign of superfluidity.
The photon fluid has another property that distinguishes it from traditional statesof light. The chemical potential of a photon fluid is nonvanishing, since,µ=∂$Ψ|H|Ψ#)= 0 .
∂N(5.41)In the limit of large photon number and low temperature, the energy of the condensate
is dominated by the photon-photon interactions in the ground state mode,E ≈ N 2 V (0) ,(5.42)

which leads to a chemical potential of,µ ≈ N V (0) .

(5.43)

126
The chemical potential is due to the atom-mediated photon-photon interactions which
are not present in traditional systems of light.

Chapter 6

Evaporative Cooling of a Photon Fluid to Quantum Degeneracy

The phenomenon of coherence has played a crucial role in many areas of physics,
in particular with the various states of light [59, 58, 159]. The extraordinarily longcoherence time is the fundamental property that distinguishes laser light from ordinary
light. The extensive coherence properties of matter have recently been investigated intrapped ultracold atoms by studying de Broglie matter waves in Bose-Einstein condensates [6, 39, 19]. Coherence is a ubiquitous phenomenon: many other systems, such asa collection of polaritons [11], also show similar coherent attributes.
In the case of atoms, coherence presents itself when the atomic gas has been condensed to quantum degeneracy. Condensation is often achieved using the technique of
evaporative cooling in which the highest energy atoms are forced to escape the trap,
and those that remain rethermalize through binary collisions at lower and lower temperature. In a nonequilibrium system, continuous evaporative cooling may generatea coherent atom laser [77]. In this chapter, we demonstrate how evaporative coolingand Bose stimulated emission can be used to condense a photon fluid into a quantum
degenerate superfluid [155]. In a general sense, a photon fluid is a collection of spatiallylocalized photons which, through their interactions via a nonlinear medium, exhibithydrodynamic or similar fluidic behavior.
A condensed photon fluid should possess superfluid properties such as coherence,phase rigidity, quantized vortices, a critical Landau velocity, and the usual Bogoliubov

128
dispersion [33]. This is distinct from a normal laser, where the effective photon interactions are so weak that they are unimportant. Another significant difference with usual
lasers is that it is not necessary to have a population inversion of an internal atomic
state, and instead the phase-space compression to produce a macroscopically occupied
mode relies on an inversion of photon cavity mode populations with respect to the traditional Planck distribution. Since the relevant modes of the system are determined by
the structural properties of the cavity and not atomic quantities, the frequency of thecoherent field could be highly tunable.
Two high-reflectivity mirrors placed close together form a Fabry-Perot cavity
and support many photon modes. The closeness of the mirrors creates a cavity mode
volume which has a high aspect ratio, with large energy gaps between modes which haveadjacent longitudinal quantum numbers and small energy gaps between modes whichhave adjacent transverse quantum numbers. We utilize this to pump the cavity suchthat only one longitudinal quantum number is relevant, and the modes considered are
distinguished only by their transverse degrees of freedom. Photons in vacuum interact
weakly, therefore, it is necessary to incorporate a nonlinear medium, such as an atomic
Rydberg gas or nonlinear crystal, into the cavity to allow for atom-mediated photonphoton interactions. Since the longitudinal degree of freedom is frozen out by having
a small cavity length, an intriguing link exists between this systems and condensationin two-dimensions with the well-known physics of the Berezinskii-Kosterlitz-Thoulesstransition, which has recently been investigated in atomic gases [69, 150].

6.1

System Dynamics
The photon fluid in a Fabry-Perot cavity is governed by the Hamiltonian
Ĥ =

$
i

h̄ωi â†i âi +

$
ijkl

† †
h̄ΓC
ijkl âi âj âk âl ,

(6.1)

129

Figure 6.1: The five primary processes that occur in the Fabry-Perot cavity are presented. Photons enter the system from a pumping source. While in the cavity, thephotons recombine and self-interact. Finally, the photons can leave the cavity throughevaporation or output coupling. The combination of these processes allows for a conditions suitable for the evaporative cooling of a photon fluid.

where âi is the annihilation operator for a photon in mode i with energy h̄ωi , and ΓC
ijkl
is the scattering or collision (C) rate for photons from modes k and l into modes i andj. Note that the mode energy h̄ωi contains the bare mode energy of the photons as wellas the energy due to the linear polarization of the medium in the cavity. The collision
rates are due to the nonlinear polarizability of the medium in the cavity. This is the
Hamiltonian that was derived and examined in the previous chapter. There are twokinds of interactions, as illustrated in Fig. 6.1. Terms which do not change populations
give rise to a self-interaction energy and leave the photons in the same mode they enteredin, e.g. â†0 â†1 â1 â0 . Population transfer is due to recombination terms, such as â†0 â†2 â1 â1 .
The self-interaction terms will later be shown to have a dispersive effect on the coherenceof the photons. The recombination terms will have a crucial role in developing a systemcapable of evaporative cooling methods.
Photons can also enter and leave the cavity through irreversible pumping anddecay, which is not described by the Hamiltonian in Eq. (6.1). Pumping and decay
can be incorporated into the system dynamics through the quantum master equation

130
formalism [110]
$
$
∂ ρ̂
i
ΓD
= − [Ĥ, ρ̂] +ΓPi L̂[â†i ] +i L̂[âi ] ,
∂t
h̄
i
i

(6.2)

where ρ̂ is the density matrix operator of the system, ΓPi is the pumping rate (P ) andΓD
i is the decay rate (D) for mode i, and the Lindbladt superoperator isL̂[Ô] ≡ 2Ô ρ̂Ô† − Ô† Ôρ̂ − ρ̂Ô† Ô .

(6.3)

The pumping rates across the spectrum of system modes are dependent on the pumping
source, which for simplicity we take to be constant across all relevant modes, ΓPi → ΓP .
It is possible to consider a pumping source that has a particular width or broadness
in the modes being pumped. This can correspond to a narrow or broad bandwidthlaser that pumps the system. The pumping width corresponds to a particular energy
spectrum of modes being pumped. This can be implemented by considering all modes
below a given energy pumped with the same pumping strength, ΓP . The remaining
modes can then be considered to remain unpumped.
Evaporative cooling requires a decay rate from the cavity which is strongly dependent on the energy. A mirror with a frequency dependent reflectivity, which allows
modes to decay at different rates, can be utilized. A reflectivity which drops for large
frequencies can create conditions where high energy photons leave the cavity faster thanlow energy photons. Since a high energy photon carries away more than the average
photon energy, subsequent rethermalization through photon interactions should reduce
the system temperature, in much the same way as in the evaporative cooling of atoms.
We implement energy dependent reflectivity by dividing the decay rates into twoclasses. High energy photons (H) are described by a decay rate of ΓH and low energy
photons (L) by a decay rate of ΓL with ΓH ( ΓL . The large decay rate of the high energy
modes means they can be eliminated in the following way. Consider an interaction inwhich photons from modes i and j create a photon of low energy in mode k and a photonof high energy in mode l. The rapid decay of mode l causes this process to be effectively

131
irreversible. This is described most simply by extracting such processes which involve ahigh energy mode from the reversible Ĥ term in Eq. 6.2, and adding them back to thedensity matrix evolution as an irreversible term of superoperator form
ΓE L̂[â†k âi âj ] ,(6.4)

where ΓE ≡ (ΓC )2 /ΓH is the evaporation rate (E), taken to be constant. This evaporation rate comes from the master equation formalism in which the correlationR$âl (t)â†l (0) = e−Γ t ,(6.5)

expresses the fast recombination rate involved. The system Hilbert space is now reduced
and the high energy modes no longer appear in the theory explicitly.
Even with these simplifications, the problem is still intractable in general. The
dimensionality of the Hilbert space grows exponentially with the number of photons in
the lower energy modes of the system; the number of states of a system with N photons
in M modes scales as(N + M )!
.
N !M !

(6.6)

Most of this complexity is uninteresting from the point of view of evaporative cooling ofphotons since it arises from all possible pathways to redistribute and entangle photons
amongst the lowest energy modes. Consequently, we divide the system into what wewill refer to as plaquettes, each plaquette consisting of a pair of modes i and j, which
may interact to produce a photon in the lowest energy mode k = 0, and a photon in ahigh energy mode l which is lost through the mirrors rapidly. Adding the contributions
from D distinct plaquettes as an incoherent sum of pumping rates into the ground state,
the density matrix equation of motion in the interaction picture is
∂ ρ̂
∂t

= +ΓP

$

L̂[â†i ] + ΓL

i

−iΓC

$
i

$

L̂[âi ](6.7)

i[â†0 â†i âi â0 , ρ̂] + ΓE$$ij%L̂[â†0 âi âj ] ,132
where the sum over $i, j# implies modes i and j are from the same plaquette.
Amongst the populations, the evolution time scales are not all equivalent. If only
one or no atoms are in a plaquette, a slow evolution takes place on time scales governedby ΓP and ΓL . If two atoms are in a plaquette, however, and providing ΓE ( ΓP , ΓL ,they will rapidly collide and form a ground state photon and a high energy photon
which will evaporate away. Consistent with this is that we may neglect the possibilityof three or more photons in a plaquette.
The basis states needed then describe the number of photons in each plaquetteand the number of photons in the lowest energy mode. The precise occupation of themodes is not important, only how many plaquettes have a given number of photons.
For instance, an arbitrary basis state |Ψ# can be uniquely identified with a reduced
computational state described by|Ψ# → |n, abcd# ≡ |n, (00)a (01)b (11)c (02)d # ,(6.8)

where n is the number of ground state photons, a is the number of plaquettes withno photons, b is the number of plaquettes with one photon, and c and d are the number of plaquettes with two photons in different modes or the same mode, respectively.
When considering the normalization of the wave function it is necessary to consider themultiplicity of basis states that all correspond to the same computational basis state.

6.2

Photon Number Distribution
The photon number distribution is an important quantity to characterize thephoton state since a coherent state has a Poissonian distribution while a thermal state
has an exponentially decaying distribution. The populations of all states |Ψ# that canbe expressed as |n, abcd# are given by
Pn,abcd ≡ $Ψ|ρ|Ψ# .

(6.9)

133
We now take advantage of the natural separation of time scales by adiabaticallyeliminating the fast processes. This is implemented by solving for the fast variables
in steady state and substituting back into the equations for the slow variables. The
populations of states with a plaquette with two photons are then approximately givenby their steady state values
Pn,ab10 =Pn,ab01 =2 ΓPP,n + 1 ΓE n,a(b+1)00
1 ΓPP,n + 1 ΓE n,a(b+1)00(6.10)(6.11)

as determined by Eq. 6.7 and Eq. 6.9. The evolution of the remaining populations isthen determined by
∂Pn,a
∂t

= −2(n + D − a)ΓL Pn,a(6.12)−2(n + 1 + 3D − a))ΓP Pn,a+2(n + 1)ΓL Pn+1,a + 2nΓP Pn−1,a+2(D − a)ΓP Pn,a+1 + 4aΓL Pn,a−1+12aΓP Pn−1,a−1 ,where Pn,a ≡ Pn,a(D−a)00 . When normalizing the populations it is necessary to properly
weight the individual computational basis states,
D−1$

2D−a D!
1 = Pn,D +a!(D − a − 1)!
a=0

!

1
2 ΓP
+
D − a n + 1 ΓE"

Pn,a ,(6.13)

due to the many basis states that correspond to the same computational basis state.
The number distribution of photons in the low energy mode is highly dependent
on the incoherent pumping rate and the number of plaquettes pumped. The pumping
rate can be designated by an effective temperature, T ,e−h̄ω/kB T ≡

n̄
ΓP
=
,
L
Γ
n̄ + 1(6.14)

where n̄ is the number of photons that would be in each mode if the modes were isolated
and non-interacting and ω is the bare mode energy. This temperature corresponds to

134
0.1

kB T = 0h̄ω(Poissonian)kB T = 0.8h̄ωProbability

kB T = 1.5h̄ωkB T = 5.6h̄ωkB T = 11h̄ωkB T = 21h̄ω

0
0

10

20
30
Photon Number

40

Figure 6.2: Photon number probability distribution of the lowest energy photon mode.
As the pumping power per mode is decreased and more modes are pumped, denotedby decreasing T , the distributions change from an exponential decay to Poissonian,indicating a transition from a thermal to a superfluid state. The systems are describedby D = {0, 1, 3, 20, 60} pumped plaquettes and a pumping rate that sets 20 photons in
the lowest energy mode.

the occupancy of the external pumping modes. In order to make effective comparison
between different temperatures, we hold constant the average number of photons in thelowest energy mode, and vary both the temperature and correspondingly the energy
width of the pump in order to maintain this condition. It would be expected that for
an intense, narrow bandwidth pump connected only with the lowest energy mode of
the cavity, the high temperature and incoherent characteristics of the pump would be
imprinted on the photons in that mode. However, for a weaker but broader pump, in
which a wider spectrum of modes are pumped but at a slower rate so that the average
photon number remains the same, the lowest energy mode is not as influenced by thetemperature of the pumping modes. The occupation of the lowest energy mode is due
mainly to the recombination and stimulated emission which should create a coherent
system.
Figure 6.2 presents the number distributions of photons in the lowest energy
mode during steady state operation for systems with 20 photons in the low energy

135
mode. With a large effective temperature, as denoted in the figure by kB T = 21h̄ω, thenumber distribution is essentially what would be expected for a thermal source. As theeffective temperature decreases, by decreasing the pumping strength and increasing thenumber of pumped plaquettes, the effects of bosonic amplification dramatically increase.
This leads to a nearly Poissonian number distribution associated with zero temperaturecondensates. The key to converting the incoherent thermal pump into a coherent systemis that the occupation of the lowest energy mode comes from stimulated emission fromhigher energy modes and not the thermal pump. With over half of the photons in
the lowest energy mode, the photons have entered quantum degeneracy and formed acondensate.

6.3

Photon Spectrum
The spectrum of the photons is also an important quantity to examine as itprovides the coherence time and spectral width. Since the mode separation is greaterthan the collisional interaction energy and the system is in steady state, the fluctuation
spectrum is determined byS(ν) =

# ∞

−∞

:

;

â†0 (τ )â0 (0);,dτ e−iντ : †â0 (0)â0 (0)(6.15)

which is the Fourier transform of a two-time correlation function. The correlation
function can be expressed as$â†0 (t)â0 (0)# =$√nCn,{i} (t) ,(6.16)n,{i}with
Cn,{i} (t) ≡ T r[ρ̂|n, {i}#$n − 1, {i}|(t)â0 (0)] ,(6.17)

where n is the number of ground state photons and {i} represents the configuration of
the remaining modes. According to the quantum regression theorem [110], the equation

136
of motion for the quantity Cn,{i} is given by∂Cn,{i}∂t=$

Mn{i}m{j} Cm,{j} ,(6.18)

m{j}where the matrix elements are determined from the off-diagonal density matrix equationof motion∂ρn{i},(n−1){i}
∂t

=

$

Mn{i}m{j} ρ(m−1){j},m{j} ,(6.19)

m{j}as given by Eq. 6.7. Since this is a reduced set of the full density matrix equationsof motion it is merely necessary to replace these matrix elements into Eq. (6.18) with√
the initial condition Cn,{i} (0) = nPn,abcd (0) where |n, {i}# corresponds to the reduced
basis state |n, abcd#. As with the populations, it is convenient to work in the reduced
computational basis states. Therefore, we consider quantities of the form
Cn,abcd (t) ≡ T r[ρ̂|n, abcd#$n − 1, abcd|(t)â0 (0)] .

(6.20)

States with two excited photons in one plaquette are given by their steady state values
Cn,ab10 =Cn,ab01 =4ΓP Cn,ab+100,(2n + 1)ΓE + i(2n + b)ΓC
4ΓP Cn,ab+100.
2(2n + 1)ΓE + i(2n + b)ΓC

(6.21)(6.22)The remaining terms are found by the equation of motion
∂Cn,a
∂t

= −ΓP (2n − 1 + 6D − 2a)Cn,a−ΓL (2D − 2a + 2n − 1)Cn,a−i(2n − 2 + D − a)ΓC Cn,a&+2 n(n − 1)ΓP Cn−1,a&+2 n(n + 1)ΓL Cn+1,a+2ΓP (D − a)Cn,a+1 + 4ΓL aCn,a−1

8a n(n − 1)ΓE ΓP Cn−1,a−1+(2n − 1)ΓE + i(2n + D − a − 2)ΓC%
32a n(n − 1)ΓE ΓP Cn−1,a−1+,2(2n − 1)ΓE + i(2n + D − a − 2)ΓC%

(6.23)

137

Figure 6.3: The spectrum is presented with vanishing collision rates. As the photon
number is increased, the spectral width narrows. In the inset, the inverse of the spectral
width is proportional to average photon number as expected from the Schawlow-Townes
linewidth.

where Cn,a ≡ Cn,a(D−a)00 and initial condition Cn,a (0) =√nPn,a (0). With this termi-

nology,$â†0 (t)â0 (0) =D$
2D−a−1 (D − 1)!

a=0

a!(D − a − 1)!

Cn,a (t) ,(6.24)

which deals with the multiple number of basis states that correspond to the same computational basis state. Note that collisional terms, with prefactor ΓC , and evaporation
terms, with prefactor ΓE , are now present, contrary to the equations of motion for the
populations.
There are two driving forces of decoherence, incoherent transfer through the mirrors and phase dispersion due to collisions, which both limit coherence time. The
decoherence due to transfer through the mirrors can be examined by setting a vanishingcollision rate. Due to bosonic amplification, the photons are able to support a statewith a spectral width much narrower than that given by the linewidth of the photon
decay.
Figure 6.3 presents the spectrum for vanishing collision rates for systems with 20
pumped plaquettes and different average photon numbers. As the photon number is138
increased due to a higher pumping rate, the spectral width decreases due to Bosonic
amplification. The spectral width can be much narrower than the width associatedwith the decay rate through the mirrors. As shown in the inset of Fig. 6.3, the spectral
width, defined as half width at half maximum, is inversely proportional to the average
photon number as expected by the Schawlow-Townes linewidth [148].
The effects of non-vanishing collision rates are of course an important aspect of
the system. Figure 6.4 presents the fluctuation spectrum for a system with an averageof 100 photons as the collision rate is increased. There is a mean-field interaction energy
shift of the frequency since the energy of the photons is no longer just the bare modeenergy. The spectral width rapidly grows as the collision rate becomes larger. One
possible way to remedy this problem is to create a system that has a small collision rate
for the low energy photons but still possesses a large recombination rate for the highenergy photons. Such a system can be manufactured by creating a spatially dependent
nonlinear medium. The low energy photons are more localized to the center of the cavity
than the high energy photons. Therefore, if the medium is less dense in the center of
the trap than toward the edges, the cavity will allow for fast recombinations of highenergy photons to promote rapid evaporative cooling but small collision rates for low
energy photons to promote coherent evolution.

6.4

Summary and Outlook
In this chapter, we examined the application of evaporative cooling to a photon

fluid trapped in a Fabry-Perot cavity. The evaporation occurs due to an energy dependent reflectivity of the cavity mirrors. The presence of atom-mediated photon-photon
interactions creates a situation in which rethermalization can occur. The amplificationof the lowest energy mode due to stimulated emission can create a coherent system,
even though the cavity is pumped with an incoherent source.
The condensation of a photon fluid using evaporative cooling provides a new

139

Figure 6.4: The spectrum is presented for systems with different collision rates. Larger
the collision rates lead to wider spectra. The spectrum is also increased from the baremode energy due to mean-field interaction shifts.

mechanism for creating a high-intensity narrow-spectrum coherent beam of light. Unlikea laser, the photon fluid, due to the presence of interactions, is expected to be superfluid
and should possess for example classic properties such as phase rigidity and the ability
to support quantized vortices. The key to the evaporative cooling mechanism in aphoton fluid is that the population of the lowest energy mode be due to the stimulated
emission from higher modes. The photon number distribution of the lowest energy mode
then becomes Poissonian instead of exponentially decaying, a signature of a transition
from a thermal to a coherent field. The spectrum of the photons also points towards along coherence time. The width of the spectrum narrows as the inverse of the photon
population as expected from the Schawlow-Townes linewidth.
Future work includes examining larger systems and models that more completely
describe experimental systems. The specific energy structure of the cavity modes canalter the precise coherence times. After creating a condensed photon fluid, it would be
instructive to examine the many superfluid properties that it should possess.

Appendix A

Jacobian Elliptic Functions

A brief review of the Jacobian elliptic functions [1, 113] is given. Of the 12 elliptic
functions, there are only six that are normalizable. Of these six, only three represent a
different physical form, sn, cn and dn. However, they are still related by,cn2 (x, k) = 1 − sn2 (x, k) ,(A.1)

dn2 (x, k) = 1 − k2 sn2 (x, k) ,(A.2)

where k is the elliptic modulus. The six non-normalizable elliptic functions can also be
reduced through a phase shift to three with different forms, ns, ds, and cs, which arealso related by,cs2 (x, k) = ns2 (x, k) − 1 ,(A.3)

ds2 (x, k) = ns2 (x, k) − k2 .

(A.4)

The normalizable and nonnormalizable functions can be related throughsn2 (iK(1 − k2 ) + x, k) = k2 ns2 (x, k) ,(A.5)

where K(x) is the complete elliptic integral of the first kind. Therefore, the square ofany elliptic function can be related linearly to sn2 (x, k).
In the limit that k → 0, the Jacobi elliptic functions approach the traditional
trigonometric functions. The sn function becomes the sine function. The cn function

141
approaches the cosine function. The dn function approaches unity. As k → 1, the Jacobi
elliptic functions become the hyperbolic functions. sn becomes the hyperbolic tangent
function. Both the cn and dn functions approach the hyperbolic secant function. The
period of these functions, as given by 4K(k), goes from 2π to infinity as k goes fromzero to one.

Appendix B

Completeness of Constant Potential Solution Set

It is possible to prove that Eq. (3.9) is a Jacobian elliptic integral of the first kindand can therefore be inverted to produce the Jacobian elliptic functions. If arbitrary
parameters are used, Eq. (3.9) becomes,#

1%

A3 ρ3 + A2 ρ2 + A1 ρ + A0dρ = x + x0 ,(B.1)

where the Ai ’s are real constants. The cubic polynomial can be factored to give,#1dρ = x + x0 ,(ρ + B1 )(ρ + B2 )(ρ + B3 )

%

(B.2)

where at least one of the constants Bi must be real. Without loss of generality we may
take the real constant as B1 . The substitution ρ = y 2 − B1 is then made in Eq. (B.2)
to yield#2%(y 2 + (B2 − B1 ))(y 2 + (B3 − B1 ))dy = x + x0 .

(B.3)

This is the general form of the elliptic integral of the first kind [1] and therefore gives
C1 el−1 (C2 y, C3 ) = x + x0 ,(B.4)

where the Ci ’s are constants and el is one of the twelve elliptic functions. This can beinverted and ρ replaced to produceρ + B1 = C2−1 el(C1−1 (x + x0 ), C3 ) ,

%

(B.5)

143
or finally,ρ = C2−2 el2 (C1−1 (x + x0 ), C3 ) − B1 .

(B.6)

Since the square of any elliptic function can be related linearly to the square of sn, only
one independent solution of the form
ρ = A sn2 (bx + x0 , k) + B ,need be considered.

(B.7)

Appendix C

Atomtronics Calculational Details

In the following, we discuss the calculation of the current-voltage characteristics
presented in Figs. 4.6, 4.9 and 4.11. The role of the system size is considered in section
C.1. The results of this discussion justify the use of small systems in the calculationof the current response to a given voltage. Our approach is developed in sections C.2,
C.3 and C.4. These three sections discuss the effect of the battery on the system,
the character of the resulting current carrying state and the way these two aspects
are used as a starting-point in formulating our method to calculate the current-voltage
characteristics. The details of the algorithm used for wires and diodes are laid out insection C.5. Finally, the calculation for the transistor is discussed in C.6.

C.1

Small System Considerations
There are three primary concerns about the system size when performing cal-culations with a small lattices: effects on the phase diagram, effects at the boundary
between two phases and effects on the current. All of these effects we find to be small.
Although phase diagrams of ultracold bosons in lattices, strictly speaking, referto infinite systems, it is possible to describe phase diagrams for much smaller systems
in a similar manner. For example, phase diagrams can be created for systems thatare only four or six sites large. For a small system, the Mott insulator phase can be
defined as a state with an integer filling of the lattice. The superfluid phase is defined

145
by a state with noninteger filling of the lattice. Note that adding an extra atom to asmall system does require slightly more energy than adding an extra atom to an infinite
system because in a large system the atom can distribute itself among more sites thanfor a small system. This, however, does not appreciably affect the characteristics of thedifferent phases. The phase diagrams for small systems are similar to that presented
for an infinite system in Fig. 4.2. The tip of the lobe is moved only slightly and thewidth of the superfluid region, where the lattice has non-integer filling and number
fluctuations are large, only changes from 2J for a system of only two sites to 4J for
an infinite lattice. The width of the superfluid region is already 3.2J for six sites. In
addition, the number fluctuations characteristic of superfluid and insulating phase do
not depend on the system size and only vary on the order of a few percent between small
and large systems of the same filling. Similarly, the doped phase diagrams presented
in Fig. 4.5 do not change considerably when calculated for small systems. The phase
diagrams of these systems all possess the same qualitative feature that the insulatingzone of the N-type doped lattice is pushed downwards while the insulating zone of the
P-type doped lattice is pushed upwards. In conclusion, we find that small systems canbe used to replicate the infinite lattice relations µ(n) plotted in Figs. 4.2 and 4.5.
Adjoining P-type and N-type materials creates a junction at which boundary
effects could alter each component. However, we have found that the boundary only
minimally affects one to two sites of each component in the limit U ( J. We haveexamined the ground state of a PN configuration that extends over 20 lattice sites andcompared each side with P-type and N-type materials in the absence of a junction. The
difference in the occupation numbers and their fluctuations between the two systems ison the order of a few percent for the first site at the junction and is down to a tenthof a percent by the third site. In conclusion, using the expression µ(n) of P-type andN-type materials to describe each side of a PN-junction individually in a small system
is justified because boundary effects due to the interface between the two materials are

146
small.
Finally, the third reason why small systems can be used to model larger systems isthat the net current is not significantly affected by the size of the system. We have verified this numerically by comparing the exact dynamical evolution in one-dimensionallattices of between two and ten sites at constant filling with periodic boundary conditions. In the initial state, the density is modulated with period 2d and there is aconstant phase difference between neighboring sites. This initial state is an exampleof the kind of running wave that would be created by a battery (as will be discussed
below). The evolution of this state with the Hamiltonian Eq. (4.1) within the two-state
approximation yields a net current that does not significantly depend on the number oflattice sites. This property is reminiscent of the current carried by a Bloch state whichis independent of the system size at fixed quasi-momentum. A good estimate of thecurrent response can therefore be obtained by considering a very small number of links.
In this way, it is possible to circumvent calculating the dynamics of a large lattice. Notethat even within the two-state approximation such a dynamical calculation is feasible
only for systems with a small number of sites.

C.2

Battery Contacts
To power a circuit, the system is connected to an atomtronic battery by bringing

two reservoirs of different chemical potential µL and µR into contact with the two endsof the lattice. The values of µL and µR determine the particle transfer ∆n which is
injected through one contact and removed through the other. The transfer ∆n follows
from the relations µL (n + ∆n) and µR (n − ∆n) with n the filling of the system at zerovoltage. Note, the chemical potentials µL and µR should be chosen appropriately toensure that the average filling is kept at n.
The net current obtained by applying a voltage is given by the average numberof particles passing through the system per unit time after a transient time has passed.

147
The net current not only depends on the particle transfer ∆n as determined by thevalues of µ, µL and µR , but also on the strength of the coupling between battery andsystem. The latter is set by the properties of the battery contact and determines therate at which transfers of magnitude ∆n take place.

C.3

Maximum current solution
The maximum current is obtained when subsequent transfers of the battery areseparated by the time it takes for one transfer of ∆n to free up the site to which thebattery is connected. Such a situation is achieved for transfer rates that are close to
J/h̄. This is the regime we focus on. It yields the largest currents and reveals the limits
on the currents that are due to the properties of the material rather than due to theproperties of the battery contact.
The state reached after a time of transient behavior takes the form of a runningdensity wave. The maxima (minima) of the wave take the value n + ∆n (n − ∆n) and,hence, the wave amplitude depends on the voltage. The wavelength is given by twice thelattice period at maximum current. Note that in the limit of small voltages, the wave
created by the contact with the battery can be described by the elementary excitationsof the gas.

C.4

Current Response Calculation
The computation of the full dynamics of a realistic system would require a simu-

lation of a large lattice interacting with several reservoirs, one at each battery contact.
Instead, we use a simplified small-system approach. This approach does not allow us todescribe the transient behavior. We expect the duration of the transient to be set by
the time it takes for the particle transfer ∆n from the battery to traverse the (M − 1)links of a lattice with M sites. Hence, the transient should occur on the time scale of(M − 1)(h̄/J).

148
Our calculation of the current response is performed in several steps. First, thechemical potentials of the left and right batteries are chosen. These chemical potentials
correspond to a average site filling. For a wire or a diode, a two well system is analyzed
corresponding a junction in the system. An initial state of NL atoms on the left and NR
atoms on the right, as determined by the chemical potentials and Eq. (4.2), is set witha phase corresponding to the lowest energy state with this density configuration. The
state is propagated in time. The calculated current is determined by the time average
current during half of an oscillation. This calculation corresponds to the current of a
traveling wave with min and max densities of NL and NR .

C.5

Atomtronic Wires and Diodes
As discussed above in C.1, currents can be calculated by considering small sys-tems. As a consequence, for a simple atomtronic wire the second step of the calculation,i.e. the calculation of the net current carried by a wave of amplitude ∆n and wavelength
2d, can be carried out by considering a single link.
In an atomtronic diode, it is the interface region rather than the bulk P-type or
N-type material that poses the biggest obstacle for atom flow and thus determines thenet current. The current calculation can be reduced down to just this single interface
link since the evolution is little affected by adding further surrounding sites.
For a given transfer from the battery ∆n, the current response of a single link iscalculated using an initial state with nR = n+∆n particles on the right and nL = n−∆nparticles on the left. The initial phase relation between the two sites is chosen so as toproduce the lowest energy configuration that has the given filling on each site. We thencompute the dynamical evolution yielding the time-dependent current i(t). The net
current I is obtained by calculating the time average of i(t) over the time interval ∆tand maximizing it with respect to ∆t. With ν the frequency of the density oscillation,
the time it takes for 2∆n particles to move from the initially higher populated site to149
its neighbor is given by 1/2ν. Hence, the result for the net current is approximately
given by I = 4∆nν.
To calculate the currents through an atomtronic wire (see Sec. 4.3) we consider
a single link in the absence of additional external potentials with an equilibrium occupation 1 < n < 2. The evolution of the two site system yields the time-dependent
currenti(t) =

2J(nR − nL ) sin(4Jt/h̄) .
h̄

(C.1)

To obtain the net current we maximize the time average of Eq. (C.1) over a time interval∆t. This yields
J
I = 1.45 (nR − nL )
h̄

(C.2)

with ∆t = 0.58h̄/J. This result is close to the value I = 4∆nν with the Rabi frequency
ν = 2J/h̄π. From the result for the relation I(∆n) we obtain the current-voltage
dependence I(∆µ) using Eq. (4.2) for an infinite lattice.
The current characteristics of a diode (see section 4.4) is calculated in the same
way with a potential step included between the two lattice sites. This potential step
ensures that the two sides have the desired difference in occupation at equilibrium. In
this configuration, the calculation yields a slightly more complicated current relation.

C.6

Transistor
The calculation for the transistor are performed in a similar way as that for a wireor a diode. Similar to the case of the diode, the current response to particle transferfrom the battery is determined by the atom flux through the two interfaces between
emitter/base and base/collector. Moreover, as discussed earlier, the current response
has a small dependence on the number of links. Therefore, the current can be calculated
by considering a three-site system, the two links representing the two links at the twointerfaces of the NPN system. An external potential is added to raise the middle site

150
playing the role of the base while the outer sites represent the emitter and collector.
The initial state is prepared as the lowest energy state with nE atoms in theemitter and nC = 1.5 atoms in the collector. The filling on the base is determinedfrom the height of the potential on the base and the condition of being in the lowest
energy state with fillings nE and nC on the emitter and collector. At equilibrium, thebase is occupied by nB = 1.005 atoms. The current of atoms passing through thebase from collector to emitter is given by iC = (iCB + iBE )/2 with iCB the current
from collector to base and iBE the current from base to emitter. As in the previous
sections, we calculate the net collector current IC by maximizing the time average of iC .
Since we are interested in the µE -dependence of IC at fixed VBC , the initial collector
occupation is kept fixed at nC = 1.5 while the emitter occupation nE is varied in therange 1 ≤ nE < 1.5.
The base current is set to be IB = Γ(nB − 1), where Γ must satisfy the conditionΓ " J/h̄. The quantity Γ describes the weak coupling with the battery at the basecontact which is due to the thinness of the base. The data for IB displayed in Fig. 4.11(b)is obtained using Γ = 0.01J/h̄. The number of current carriers involved in IB is givenby (nB − 1) because the chemical potential of the base contact is set to be in the n = 1
insulator region. Note that our calculation neglects the effect of IB on IC . This isjustified because the weak coupling at the base contact implies that the reduction ofnB by IB is small. Overall, the transistor current gain relies on two factors. The baseregion must be thin and the emitter must have an equilibrium filling significantly larger
than the base.

Bibliography[1] M. Abramowitz and I. A. Stegun, editors. Handbook of Mathematical Functions.
National Bureau of Standards, Washington, D.C., 1964.
[2] S. A. Alexander and R. L. Coldwell. Solution of the kronig-penney model using
variational monte carlo. Int. J. Quantum Chem., 86(4):325, Feb. 2002.
[3] J. F. Allen and A. D. Misner. Flow of liquid helium ii. Nature, 141:75, 1938.
[4] S. L. Altmann. Band Theory of Metals: The Elements. Pergamon PRess, Oxford,
1970.
[5] B. P. Anderson and M. A. Kasevich. Macroscopic quantum interference fromatomic tunnel arrays. Science, 282(5394):1686–1689, 1998.
[6] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell. Observation of bose-einstein condensation in a dilute atomic vapor. Science,
269:198, July 1995.
[7] T. Anker, M. Albiez, B. Eiermann, M. Taglieber, and M. K. Oberthaler. Linear and nonlinear dynamics of matter wave packets in periodic potentials. Opt.
Express, 12(1):11, Jan. 2004.
[8] P. D. Ankrum, editor. Semiconductor Electronics. Prentice-Hall, EnglewoodCliffs, N.J., 1971.
[9] N. W. Ashcroft and N. D. Mermin. Solid State Physics. W. B. Saunders Company,
1976.
[10] G. E. Astrakharchik and L. P. Pitaevskii. Motion of a heavy impurity through abose-einstein condensate. Phys. Rev. A, 70:013608, July 2004.
[11] R. Balili, V. Hartwell, D. Snoke, L. Pfeiffer, and K. West. Bose-einstein condensation of microcavity polaritons in a trap. Science, 326:1007, 2007.
[12] R. Bhat, M. Holland, and L. Carr. Bose-Einstein Condensates in Rotating Lattices. Physical Review Letters, 96(6):60405, 2006.
[13] R. Bhat, B. Peden, B. Seaman, M. Krämer, L. Carr, and M. Holland. Quantizedvortex states of strongly interacting bosons in a rotating optical lattice. PhysicalReview A, 74(6):63606, 2006.

152[14] M. M. Bogdan, A. S. Kovalev, and I. V. Gerasimchuk. Dynamics and stabilityof localized modes in nonlinear media with point defects. Low Temp. Phys.,23(2):145–152, Feb. 1997.
[15] K. Bongs and K. Sengstock. Physics with coherent matter waves. Reports on
Progress in Physics, 67(6):907–963, 2004.
[16] S. N. Bose. Plancks gesetz und lichtquantenhypothese. Z. Phys., 26:178, 1924.
[17] F. Bowman. Introduction to Elliptic Functions, with Applications. Dover, New
York, 1961.
[18] C. C. Bradley, C. A. Sackett, and R. G. Hulet. Analysis of in situ images of
bose-einstein condensates of lithium. Phys. Rev. A, 55:3951, 1997.
[19] C. C. Bradley, C. A. Scakett, J. J. Tollett, and R. G. Hulet. Evidence of boseeinstein condensation in an atomic gs with attractive interactions. Phys. Rev.
Lett., 75(9):1687–1690, Aug 1995.
[20] G. Brennen, C. Caves, P. Jessen, and I. Deutsch. Quantum Logic Gates in Optical
Lattices. Physical Review Letters, 82(5):1060–1063, 1999.
[21] H. Briegel, T. Calarco, D. Jaksch, J. Cirac, and P. Zoller. Quantum computing
with neutral atoms. Journal of Modern Optics, 47(2-3):415–451, 2000.
[22] J. C. Bronski, L. D. Carr, R. Carretero-González, B. Deconinck, J. N. Kutz,
and K. Promislow. Stability of attractive bose-einstein condensates in a periodic
potential. Phys. Rev. E, 64:056615–1–9, 2001.
[23] J. C. Bronski, L. D. Carr, B. Deconinck, and J. N. Kutz. Bose-einstein condensatesin standing waves: The cubic nonlinear schrodinger equation with a periodic
potential. Phys. Rev. Lett., 86(8):1402, Feb. 2001.
[24] J. C. Bronski, L. D. Carr, B. Deconinck, J. N. Kutz, and K. Promislow. Stability of repulsive bose-einstein condensates in a periodic potential. Phys. Rev. E,
63:036612, Feb. 2001.
[25] P. Capuzzi and E. S. Hernandez. 3 he impurity in a bose-einstein condensate. Phys.
Rev. S, 62:023603, 2000.
[26] L. D. Carr, C. W. Clark, and W. P. Reinhardt. Stationary solutions of the onedimensional nonlinear schrodinger equation. i. case of repulsive nonlinearity. Phys.
Rev. A, 62:063610, December 2000.
[27] L. D. Carr, C. W. Clark, and W. P. Reinhardt. Stationary solutions of the
one-dimensional nonlinear schrodinger equation. ii. case of attractive nonlinearity.
Phys. Rev. A, 62:063611, December 2000.
[28] L. D. Carr, J. N. Kutz, and W. P. Reinhardt. Stability of stationary states in the
cubic nonlinear schrdinger equation: Applications to the bose-einstein condensate.
Phys. Rev. E, 63:066604, May 2001.

153[29] L. D. Carr, M. A. Leung, and W. P. Reinhardt. Dynamics of the bose-einsteincondensate: quasi-one-dimension and beyond. J. Phys. B, 33:3983, 2000.
[30] L. D. Carr, K. W. Mahmud, and W. P. Reinhardt. Tunable tunneling: An application of stationary states of bose-einstein condensates in traps of finit depth.
Phys. Rev. A, 64:033603–1–9, 2001.
[31] F. Cataliotti, S. Burger, C. Fort, P. Maddaloni, F. Miardi, A. Trombettoni,
A. Smerzi, and M. Inguscio. Josephson junction arrays with bose-einstein condensates. Science, 293(5531):842, 2001.
[32] R. M. Cavalcanti, P. Giacconi, G. Pupillo, and R. Soldati. Bose-einstein condensation in the presence of a uniform field and a pointlike impurity. Phys. Rev. A,
65:053606, 2002.
[33] R. Y. Chiao and J. Boyce. Bogoliubov dispersion relation and the possibility ofsuperfluidity for weakly interacting photons in a two-dimensional photon fluid.
Phys. Rev. A, 60(5):4114–4121, November 1999.
[34] A. P. Chikkatur, A. Gorlitz, D. M. Stamper-Kurn, S. Inouye, S. Gupta, and
W. Ketterle. Suppression and enhancement of impurity scattering in a boseeinstein condensate. Phys. Rev. Let., 85(3):483, 2000.
[35] D. N. Christodoulides, F. Lederer, and Y. Silberberg. Discretizing light behavior in
linear and nonlinear waveguide lattices. Nature (London), 424(6950):817, August
2003.
[36] C. Conti, M. Peccianti, and G. Assanto. Route to nonlocality and observation ofaccessible solitons. Phys. Rev. Let., 91(7):073901, 2003.
[37] A. Daley, S. Clark, D. Jaksch, and P. Zoller. Numerical analysis of coherent
many-body currents in a single atom transistor. Physical Review A, 72(4):43618,
2005.
[38] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari. Theory of bose-einsteincondensation in trapped gases. Rev. Mod. Phys., 71:463, April 1999.
[39] K. B. Davis, M. O. Mewes, M. R. Andrews, N. J. Vandruten, D. S. Durfee, D. M.
Kurn, and W. Ketterle. Bose-einstein condensation in a gas of sodium atoms.
Phys. Rev. Lett., 75(22):3969–3973, Nov. 1995.
[40] B. Deconinck, B. A. Frigyik, and J. N. Kutz. Dynamics and stability of boseeinstein condensates: The nonlinear schrodinger equation with periodic potential.
J. Nonlinear Sci., 12:169–205, 2002.
[41] J. Denschlag, J. Simsarian, H. Haeffner, C. McKenzie, A. Browaeys, D. Cho,
K. Helmerson, S. Rolston, and W. Phillips. A bose-einstein condensate in an
optical lattice. Journal of Physics B Atomic Molecular and Optical Physics,
35(14):3095–3110, 2002.

154[42] L. DeSarlo, L. Fallani, J. E. Lye, M. Modugno, R. Saers, C. Fort, and M. Inguscio.
Unstable regimes for a bose-einstein condensate in an optical lattice. e-printcond-mat/0412279v1, 2004.
[43] D. Diakonov, L. M. Jensen, C. J. Pethick, and H. Smith. Loop structure of thelowest bloch band for a bose-einstein condensate. Phys. Rev. A, 66:013604, 2002.
[44] E. A. Donley, N. R. Claussen, S. T. Thompson, and C. E. Wieman. Atom-molecule
coherence in a bose-einstein condensate. Nature, 417:529, 2002.
[45] B. Eiermann, T. Anker, M. Albiez, M. Taglieber, P. Treutlein, K. P. Marzlin,
and M. K. Oberthaler. Bright bose-einstein gap solitons of atoms with repulsive
interaction. Phys. Rev. A, 92(23):230401, 2004.
[46] B. Eiermann, P. Treutlein, T. Anker, M. Albiez, M. Taglieber, K.-P. Marzlin, and
M. K. Oberthaler. Dispersion management for atomic matter waves. Phys. Rev.
Lett., 91(6):060402, 2003.
[47] A. Einstein. Quantentheorie des einatomigen idealen gases. Sitzungsber. Kgl.
Preuss. Akad. Wiss., page 261, 1924.
[48] T. Esslinger and K. Molmer. Atoms and molecules in lattices: Bose-einsten condensates built on a shared vacuum. Phys. Rev. Lett., 90(16):160406, Apr. 2003.
[49] L. Fallani, F. S. Cataliotti, J. Catani, C. Fort, M. Modugno, M. Zawada, and
M. Inguscio. Optically induced lensing effect on a bose-einstein condensate expanding in a moving lattice. Phys. Rev. Lett., 91(24):240205, Dec. 2003.
[50] L. Fallani, L. De Sarlo, J. E. Lye, M. Modugno, R. Saers, S. Fort, and M. Inguscio.
Observation of dynamical instability for a bose-einstein condensate in a moving
1d optical lattice. Phys. Rev. Lett., 93(14):140406, Oct. 2004.
[51] J. Felber, R. Gahler, and C. Rausch. Matter waves at a vibrating surface: Transition from quantum-mechanical to classical behavior. Phys. Rev. A, 53(1):319,
1996.
[52] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher. Boson localization and the superfluid-insulator transition. Phys. Rev. B, 40:546–570, 1989.
[53] R. Folman, P. Kruger, J. Schmiedmayer, J. Denschlag, and C. Henkel. Microscopic
atom optics: From wires to an atom chip. Adv. At. Mol. Opt. Phys, 48:263, 2002.
[54] C. Fort, F. Cataliotti, J. Catani, L. De Sarlo, L. Fallani, J. Lye, M. Modugno,
R. Saers, and M. Inguscio. Dynamics of bloch matter waves and its breakdown: abose–einstein condensate in a moving 1d optical lattice. Laser Physics, 15(3):447–
453, 2005.
[55] Y. B. Gaididei, P. L. Christiansen, K. O. Rasmussen, and M. Johansson. Twodimensional effects in nonlinear kronig-penney models. Phys. Rev. B, 55(20):55,
May 1997.

155[56] H. Gimperlein, S. Wessel, J. Schmiedmayer, and L. Santos. Ultracold atoms in
optical lattices with random on-site interactions. Phys Rev Lett, 95(17):170401,
2005.
[57] M. Girardeau. Relationship between systems of impenetrable bosons and fermions
in one dimension. J. Math. Phys., 1:516, 1960.
[58] R. J. Glauber. Coherent and incoherent states of the radiation field. Phys. Rev.,
131(6):2766–2788, September 1963.
[59] R. J. Glauber. The quantum theory of optical coherence. Phys. Rev., 130(6):2529–
2539, June 1963.
[60] J. P. Gordon. Interaction forces among solitons in optical fibers. Opt. Lett.,
8:596–598, 1983.
[61] J. P. Gordon and H. A. Haus. Random walk of coherently amplified solitons in
optical fiber transmission. Opt. Lett., 11:665–667, 1986.
[62] M. Greiner, I. Bloch, O. Mandel, T. Hänsch, and T. Esslinger. Exploring
phase coherence in a 2d lattice of bose-einstein condensates. Phys. Rev. Lett.,
87(16):160405, 2001.
[63] M. Greiner, I. Bloch, O. Mandel, T. W. Hansch, and T. Esslinger. Bose-einsteincondensates in 1d and 2d optical lattices. Appl. Phys. B, 73, 2001.
[64] M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and I. Bloch. Quantum
phase transition from a superfluid to a mott insulator in an ultracold gas of atoms.
Nature, 415:39, Jan. 2002.
[65] E. P. Gross. Structure of a quantized vortex in boson systems. Nuovo Cimento,
20:454, 1961.
[66] G. Grynberg, B. Lounis, P. Vererk, J. Y. Courtois, and C. Salomon. Quantizedmotion of cold cesium atoms in two- and three-dimensional optical potentials.
Phys. Rev. Lett., 70:2249, 1993.
[67] G. Grynberg and C. Robilliard. Cold atoms in dissipative optical lattices. Physics
Reports, 355(5):335–451, 2001.
[68] L. Guidoni, C. Triche, P. Verkerk, and G. Grynberg. Quasiperiodic optical lattices.
Phys. Rev. Lett., 79(18):3363–3366, Nov 3 1997.
[69] Z. Hadzibabic, P. Kruger, M. Cheneau, B. Battelier, and J. Dalibard. Berezinskiikosterlitz-thouless crossover in a trapped atomic gas. Nature, 441:1118–1121, June
2006.
[70] E. W. Hagley, L. Deng, M. Kozuma, J. Wen, K. Helmerson, S. L. Rolston,and W. D. Phillips. A well-collimated quasi-continuous atom laser. Science,
238(5408):1706–1709, Mar 1999.
[71] V. Hakim. Nonlinear schrodinger flow past an obstacle in one dimenstion. Phys.
Rev. E, 55(3):2835, March 1997.

156[72] J. Harms and J. P. Toennies. Experimental evidence for the transmission of 3heatoms through superfluid 4he droplets. Phys. Rev. Let., 83(2):344, 1998.
[73] H. Hasimoto and H. Ono. Nonlinear modulation of gravity waves. J. Phys. Sco.
Jpn, 33:805, 1972.
[74] H. Haus and W. S. Wong. Solitons in optical communications. Rev. Mod. Phys.,
68:423–444, 1996.
[75] A. Hemmerich and T. W. Hänsch. Two-dimesional atomic crystal bound by light.
Phys. Rev. Lett., 70:410, 1993.
[76] K. M. Hilligsoe, M. K. Oberthaer, and K.-P. Marzlin. Stability of gap solitons ina bose-einstein condensate. Phys. Rev. A, 66:063605, Dec. 2002.
[77] M. Holland, K. Burnett, C. Gardiner, J. I. Cirac, and P. Koller. Theory of an
atom laser. Phys. Rev. A, 54(3):R1757, September 1996.
[78] P. Horak, J. Courtois, and G. Grynberg. Atom cooling and trapping by disorder.
Physical Review A, 58(5):3953–3962, 1998.
[79] E. Infeld and G. Rowlands. On the stability of electron plasma waves. J. Phys.
A, 12(11):2255–2262, 1979.
[80] E. Infeld and G. Rowlands. Nonlinear Waves, Solitons and Chaos. Cambridge
University Press, United Kingdom, 2nd edition, 2000.
[81] A. Inouye, M. R. Andrews, J. Stenger, H.-J. Miesner, D. M. Stamper-Kurn, and
W. Ketterle. Observation of feshbach resonances in a bose-einstein condensate.
Nature, 392:151, March 1998.
[82] D. Jacksch, A. Brude, J. I. Cirac, C. W. Gardiner, and P. Zoller. Cold bosonicatoms in optical lattices. Phys. Rev. Lett., 81(15):3108, Oct. 1998.
[83] J. D. Jackson. Classical Electrodynamics. Wiley, New York, 1999.
[84] R. Jackson and M. Weinstein. Geometric Analysis of Bifurcation and Symmetry Breaking in a GrossPitaevskii Equation. Journal of Statistical Physics,
116(1):881–905, 2004.
[85] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller. Cold bosonicatoms in optical lattices. Phys. Rev. Lett., 81(15):3108, Oct. 1998.
[86] P. Jordan and E. Wigner. Uber das PaulischeAquivalenzverbot. Z. Physik, 47:631–
651, 1928.
[87] P. Kapitza. Viscosity of liquid helium below the λ-point. Nature, 141:74, 1938.
[88] G. Kavoulakis. Bose-Einstein condensates with attractive interactions on a ring.
Physical Review A, 67(1):11601, 2003.
[89] P. G. Kevrekidis, R. Carretero-Gonzalez, G. Theocharis, D. J. Frantzeskakis, and
B. A. Malomed. Stability of dark solitons in a bose-einstein condensate trappedin an optical lattice. Phys. Rev. A, 68:035602, 2003.

157[90] H. Kimble, M. Dagenais, and L. Mandel. Photon Antibunching in ResonanceFluorescence. Physical Review Letters, 39(11):691–695, 1977.
[91] H. Kimble and L. Mandel. Theory of resonance fluorescence. Physical Review A,
13(6):2123–2144, 1976.
[92] Y. S. Kivshar. Soliton scattering by impurities in hydrogen-bonded chains. Phys.
Rev. A, 43(6):3117, March 1991.
[93] Y. S. Kivshar, T. J. Alexander, and S. K. Turitsyn. Nonlinear modes of a macroscopic quantum oscillator. Phys. Lett. A, 278(4):225–230, 2001.
[94] Y. S. Kivshar, Z. Fei, and L. Luis Vzquez. Resonant soliton-impurity interactions.
Phys. Rev. Lett., 67(10):1177, September 1991.
[95] Y. S. Kivshar and B. Luther-Davies. Dark optical solitons: Physics and applications. Physics Reports, 298:81–197, 1998.
[96] A. G. Klein, P. Prager, P. J. Wagenfeld, H.and Ellis, and T. M. Sabine. Diffractionof neutrons and x-rays by a vibrating quartz crystal. Appl. Phys. Lett., 10(10):294,
1967.
[97] 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, andinteractions. Physical Review Letters, 94(8):80403, 2005.
[98] W. Krauth, M. Caffarel, and J. Bouchaud. Gutzwiller wave function for a modelof strongly interacting bosons. Physical Review B, 45(6):3137–3140, 1992.
[99] P. Leboeuf, N. Pavloff, and S. Sinha. Solitonic transmission of bose-einstein matterwaves. Phys. Rev. A, 68:063608, Dec 2003.
[100] A. J. Leggett. Bose-einstein condensation in the alkali gases: Some fundamental
concepts. Rev. Mod. Phys., 73(2):307, April 2001.
[101] G. Lenz, P. Meystre, and E. M. Wright. Nonlinear atom optics: General formalism
and atomic solitons. Phys. Rev. A, 50(2):1681, 1994.
[102] W. Li and A. Smerzi. Nonlinear kronig-penney model. Phys. Rev. E, 70(1):016605,
July 2004.
[103] F. London. The λ-phenomenon of liquid helium and the bose-einstein degeneracy.
Nature, 141:643, 1938.
[104] P. J. Y. Louis, E. A. Ostrovskaya, and Y. S. Kivshar. Dispersion control for matterwaves and gap solitons in optical superlattices. Phys. Rev. A, 71(2):023612, Feb
2005.
[105] P. J. Y. Louis, E. A. Ostrovskaya, M. Savage, Craig, and Y. S. Kivshar. Boseeinstein condensates in optical lattices: Band-gap structure and solitons. Phys.
Rev. A, 67:013602, Jan. 2003.

158[106] M. Machholm, A. Nicolin, C. J. PEthick, and H. Smith. Spatial period doublingin bose-einstein condensates in an optical lattice. Phys. Rev. A, 69:043604, 2004.
[107] M. Machholm, C. J. Pethick, and H. Smith. Band structure, elementary excitations, and stability of a bose-einstein condensate in a periodic potential. Phys.
Rev. A, 67(5):053613, May 2003.
[108] P. Massignan and M. Modugno. One-dimensional model for the dynamics andexpansion of elongated bose-einstein condensates. Phys. Rev. A, 67:023614, 2003.
[109] I. M. Merhasin, B. V. Gisin, R. Driben, and B. A. Malomed. Finite-band solitonsin the kronig-penney model with the cubic-quintic nonlinearity. Phys. Rev. E,
71(1):016613, Jan 2005.
[110] P. Meystre and M. Sargent III. Elements of Quantum Optics. Springer, Berlin,
3rd edition, 1998.
[111] A. Micheli, A. Daley, D. Jaksch, and P. Zoller. Single Atom Transistor in a 1D
Optical Lattice. Physical Review Letters, 93(14):140408, 2004.
[112] A. Micheli and P. Zoller. Single-atom mirror for one-dimensional atomic latticegases. Physical Review A, 73(4):43613, 2006.
[113] L. M. Milne-Thomson. Jacobian Elliptic Function Tables. Dover Publications
Inc., New York, 1950.
[114] G. Modugno, F. Ferlaino, R. Heidemann, G. Roati, and M. Inguscio. Productionof a fermi gas of atoms in an optical lattice. Physical Review A, 68(1):11601, 2003.
[115] H. Moritz, T. Stöferle, M. Köhl, and T. Esslinger. Exciting collective oscillations
in a trapped 1d gas. Phys. Rev. Lett., 91:250402, 2003.
[116] O. Morsch, J. Müller, M. Cristiani, D. Ciampini, and E. Arimondo. Bloch oscillations and mean-field effects of bose-einstein condensates in 1d optical lattices.
Physical Review Letters, 87(14):140402, 2001.
[117] E. J. Mueller. Superfluidity and mean-field energy loops: Hysteric behavior inbose condensates. Phys. Rev. A, 66(6):063603, Dec 2002.
[118] J. G. Muga and M. Buttiker. Time dependence of evanescent quantum waves.
Phys. Rev. A, 62:023808, 2000.
[119] M. Olshanii. Atomic scattering in the presence of an external confinement and agas of impenetrable bosons. Phys. Rev. Lett., 81:938–941, 1998.
[120] E. A. Ostrovskaya and Y. S. Kivshar. Localization of two-component bose-einsteincondensates in optical lattices. Phys. Rev. Lett., 92(18):180405, May 2004.
[121] Y. B. Ovchinnikov, J. H. Müller, M. R. Doery, E. J. D. Vredenbregt, K. Helmerson, S. L. Rolston, and W. D. Phillips. Diffraction of a released bose-einsteincondensate by a pulsed standing light wave. Phys. Rev. Lett., 83(2):284, July
1999.

159
[122] N. Pavloff. Breakdown of superfluidity of an atom laser past an obstacle. Phys.
Rev. A, 66:013610, July 2002.
[123] B. Peden, R. Bhat, M. Krämer, and M. Holland. Quasi-angular momentum ofBose and Fermi gases in rotating optical lattices. Journal of Physics B: Atomic,
Molecular and Optical Physics, 40(18):3725–3744, 2007.
[124] S. Peil, J. V. Porto, B. L. Tolra, J. M. Obrecht, B. E. King, M. Subbotin, S. L.
Rolston, and W. D. Phillips. Patterned loading of a bose-einstein condensate into
an optical lattice. Phys. Rev. A, 67(5):051603, May 2003.
[125] R. Pepino, J. Cooper, D. Anderson, and M. Holland. Atomtronic circuits of diodes
and transistors. cond-mat/0705.3268, 2007.
[126] L. Pezzè, L. Pitaevskii, A. Smerzi, S. Stringari, G. Modugno, E. de Mirandes,
F. Ferlaino, H. Ott, G. Roati, and M. Inguscio. Insulating behavior of a trapped
ideal fermi gas. Physical Review Letters, 93(12):120401, 2004.
[127] L. P. Pitaevskii. Vortex lines in an imperfect bose gas. Zh. Eksp. Teor. Fiz.,
40:646, 1961.
[128] M. A. Porter and P. Cvitanovic. Modulated amplitude waves in bose-einsteincondensates. Phys. Rev. E, 69(4):047201, Apr 2004.
[129] M. A. Porter, P. G. Kevrekidis, and B. A. Malomed. Resonant and non-resonant
modulated amplitude waves for binary bose-einstein condensates in optical lattices. Physica D, 196:106–123, Sep. 2004.
[130] A. Radouani. Soliton and phonon production by an oscillating obstacle in aquasi-one-dimensional trapped repulsive bose-einstein condensate. Phys. Rev. A,
70:013602, July 2004.
[131] C. A. Regal, M. Greiner, and D. S. Jin. Observation of resonance condensationof fermionic atom pairs. Phys. Rev. Lett, 92:040403, 2004.
[132] J. Reichel. Microchip traps and Bose–Einstein condensation. Applied Physics B:
Lasers and Optics, 74(6):469–487, 2002.
[133] A. M. Rey. Ultracold bosonic atoms in optical lattices. PhD thesis, University of
Maryland at College Park, 2004.
[134] A. M. Rey, K. Burnett, R. Roth, M. Edwards, C. J. Williams, and C. W. Clark.
Bogoliubov approach to superfluidity of atoms in an optical lattice. J. Phys. B,
36(5):825–841, Mar. 2003.
[135] G. Roati, E. de Mirandes, F. Ferlaino, H. Ott, G. Modugno, and M. Inguscio. Atom interferometry with trapped fermi gases. Physical Review Letters,
92(23):230402, 2004.
[136] J. L. Roberts, N. R. Claussen, J. P. Burke, Jr., C. H. Greene, E. A. Cornell, andC. E. Wieman. Resonant magnetic field control of elastic scattering in cold 85rb.
Phys. Rev. Let., 81:5109, 1998.

160[137] S. Rolston and W. Phillips.
416(6877):219–224, 2002.

Nonlinear and quantum atom optics.

Nature,[138] R. Roth and K. Burnett. Phase diagram of bosonic atoms in two-color superlattices. Phys. Rev. A, 68(2):023604, Aug 2003.
[139] A. Ruschhaupt, F. Delgado, and J. Muga. Velocity selection of ultra-cold atoms
with Fabry–Perot laser devices: improvements and limits. Journal of Physics B:
Atomic, Molecular, and Optical Physics, 38(15):2665–2674, 2005.
[140] A. Ruschhaupt and J. Muga. Atom diode: A laser device for a unidirectionaltransmission of ground-state atoms. Physical Review A, 70(6):61604, 2004.
[141] A. Ruschhaupt and J. Muga. Adiabatic interpretation of a two-level atom diode,a laser device for unidirectional transmission of ground-state atoms. PhysicalReview A, 73(1):13608, 2006.
[142] A. Ruschhaupt, J. Muga, and M. Raizen. One-photon atomic cooling with an
optical Maxwell demon valve. Journal of Physics B: Atomic, Molecular and Optical
Physics, 39(18):3833–3838, 2006.
[143] S. Sachdev, editor. Quantum Phase Transitions. CAmbridge Univ. Press, 1999,
Cambridge.
[144] J. Sacher, D. Baums, P. Panknin, W. Elsasser, and E. O. Gobel. Intensity instabilities of semiconductor-lasers under current modulation, external light injection,and delayed feedback. Phys. Rev. A, 45(3):1893–1905, Feb 1 1992.
[145] L. Salasnich, A. Parola, and L. Reatto. Effective wave equations for the dynamics
of cigar-shaped and disk-shaped bose condensates. Phys. Rev. A, 65(4):043614,
Apr 2002.
[146] L. Sanchez-Palencia and L. Santos. Bose-Einstein condensates in optical quasicrystal lattices. Physical Review A, 72(5):53607, 2005.
[147] L. Santos, M. Baranov, J. Cirac, H. Everts, H. Fehrmann, and M. Lewenstein. Atomic Quantum Gases in Kagomé Lattices. Physical Review Letters,
93(3):30601, 2004.
[148] A. L. Schawlow and C. H. Townes. Infared and opitcal masers. Phys. Rev.,
112(6):1940–1949, Dec. 1958.
[149] T. Schulte, S. Drenkelforth, J. Kruse, W. Ertmer, J. Arlt, K. Sacha, J. Zakrzewski,
and M. Lewenstein. Routes towards Anderson localization of Bose-Einstein condensates in disordered optical lattices. Phys. Rev. Lett., 95:170411, 2005.
[150] V. Schweikhard, S. Tung, and E. A. Cornell. Vortex proliferation in the berezinskiikosterlitz-thouless regime on a two-dimensional lattice of bose-einstein condensates. Phys. Rev. Lett., 99:030401, 2007.
[151] M. M. Scott, Y. K. Fetisov, V. T. Synogach, and C. E. Patton. Suppressionof microwave magnetic envelope solitons by continuous wave magnetostatic wave
signals. J. Appl. Phys., 88:4232, 2000.

161[152] B. T. Seaman, L. D. Carr, and M. J. Holland. Effect of a potential step or impurity
on the bose-einstein condensate mean field. Phys. Rev. A, 71(3):033609, March
2005. cond-mat/0410345.
[153] B. T. Seaman, L. D. Carr, and M. J. Holland. Nonlinear band structure in bose
einstein condensates: The nonlinear schrödinger equation with a kronig-penney
potential. Phys. Rev. A, 71(3):033622, March 2005.
[154] B. T. Seaman, L. D. Carr, and M. J. Holland. Period doubling, two-color lattices, and the growth of swallowtails in bose-einstein condensates. Phys. Rev. A,
72(3):033602, March 2005.
[155] B. T. Seaman and M. J. Holland. Evaporative cooling of a photon fluid to quantum
degeneracy. submitted to Phys. Rev. Lett., 2008. condmat.
[156] B. T. Seaman, M. Krämer, D. Z. Anderson, and M. J. Holland. Atomtronics:
Ultracold-atom analogs of electronic devices. Phys. Rev. A, 75(2):023615, 2007.
[157] C. P. Search and P. Meystre. Nonlinear and quantum optics of atomic and molecular fields. PROGRESS IN OPTICS, 47:139, 2005.
[158] T. B. Simpson, J. M. Liu, A. Gavrielides, V. Kovanis, and P. M. Alsing. Perioddoubling route to chaos in a semiconductor-laser subject to optical injection. Appl.
Phys. Lett., 64(26):3539–3541, Jun 27 1994.
[159] R. E. Slusher, L. W. Hollberg, B. Yurke, J. C. Mertz, and J. F. Valley. Observationof squeezed states generated by four-wave mixing in an optical cavity. Phys. Rev.
Lett, 55(22):2409–2412, November 1985.
[160] K. Stasiewicz, P. K. Shukla, G. Gustafsson, S. Buchert, B. Lavraud, B. Thide,and Z. Klos. Slow magnetosonic solitons detected by the cluster spacecraft. Phys.
Rev. Let., 90(8):085002, 2003.
[161] A. Steane. Quantum computing. Reports on Progress in Physics, 61(2):117–173,
1998.
[162] J. Stickney, D. Anderson, and A. Zozulya. Transistorlike behavior of a BoseEinstein condensate in a triple-well potential. Physical Review A, 75(1 B), 2007.
[163] T. Stöferle, H. Moritz, C. Schori, M. Köhl, and T. Esslinger. Transition from a
Strongly Interacting 1D Superfluid to a Mott Insulator. Physical Review Letters,
92(13):130403, 2004.
[164] A. A. Sukhorukov and Y. S. Kivshar. Nonlinear localized waves in a periodic
medium. Phys. Rev. E, 87(8):083901, Aug. 2001.
[165] A. A. Sukhorukov and Y. S. Kivshar. Spatial optical solitons in nonlinear photonic
crystals. Phys. Rev. E, 65(3):036609, Mar. 2001.
[166] C. Sulem and P. L. Sulem. Nonlinear Schrödinger Equations: Self-focusing
Instability and Wave Collapse. Springer-Verlag, New York, 1999.

162[167] V. I. Talanov. Self-focusing of wave beams in nonlinear media. Zh. Eksp. Teor.
Fiz., 2:31, 1965.
[168] D. Taras-Semchuk and J. M. F. Gunn. Superfluid flow past an array of scatterers.
Phys. Rev. B, 60(18):13139, November 1999.
[169] E. Taylor and E. Zaremba. Bogoliubov sound speed in periodically modulated
bose-einstein condensates. Phys. Rev. A, 68(5):053611, Nov. 2003.
[170] S. Theodorakis and E. Leontidis. Bound states in a nonlinear kronig-penney
model. J. Phys. A, 30(13):4835–4849, July 1997.
[171] R. A. Usmanov and L. B. Ioffe. Theoretical investigation of a protected quantum
bit in a small josephson junction array with tetrahedral symmetry. Phys. Rev. B,
69(21):214513, Jun 2004.
[172] J. Villavicencio, R. Romo, and S. S. y. Silva. Quantum-wave evolution in a step
potential barrier. Phys. Rev. A, 66:042110, 2002.
[173] B. Wu, R. B. Diener, and Q. Niu. Bloch waves and bloch bands of bose-einsteincondensates in optical lattices. Phys. Rev. A, 65:025601, 2002.
[174] B. Wu and Q. Niu. Nonlinear landua-zener tunneling. Phys. Rev. A, 61:023402,
2000.
[175] B. Wu and Q. Niu. Landau and dynamical instabilities of the superflow of boseeinstein condensates in optical lattices. Phys. Rev. A, 64:061603, 2001.
[176] H. C. Yuen and B. M. Lake. Nonlinear deep water waves: Theory and experiment.
Phys. Fluids, 18(8):956–960, Aug. 1975.
[177] V. E. Zakharov and A. B. Shabat. Exact theory of two-dimensional self-focusing
and one-dimensional self-modulation of waves in nonlinear media. Zh. Eksp. Teor.
Fiz., 61:118, 1971.
[178] V. E. Zakharov and A. B. Shabat. Interaction between solitons in a stable medium.
Zh. Eksp. Teor. Fiz., 64:1627, 1973.
[179] W. Zwerger. Mott-hubbard transition of cold atoms in optical lattices. Journal of
Optics B: Quantum and Semiclassical Optics, 5(2):S9–S16, 2003.