The Road Less Traveled: Resonant Bose-Einstein
Condensates via a Hyperspherical Lowest-Order
Constrained Variational ApproachbyMichelle Wynne Ching Sze
B.S., Ateneo de Manila University, 2010
M.S., University of the Philippines, Diliman, 2012

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
2019This thesis entitled:
The Road Less Traveled: Resonant Bose-Einstein Condensates via a Hyperspherical Lowest-Order
Constrained Variational Approach
written by Michelle Wynne Ching Sze
has been approved for the Department of Physics

Prof. John L. Bohn

Prof. Jose P. D’ Incao

Date

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

iii
Sze, Michelle Wynne Ching (Ph.D., Physics)
The Road Less Traveled: Resonant Bose-Einstein Condensates via a Hyperspherical Lowest-Order
Constrained Variational Approach
Thesis directed by Prof. John L. Bohn

In this work, we study the ground state properties of a system of N harmonically trappedbosons of mass m interacting with two-body contact interactions, from small to large scattering
lengths. This is accomplished in a hyperspherical coordinate system that is flexible enough todescribe both the overall scale of the gas and two-body correlations. By adapting the lowestorder constrained variational (LOCV) method, we are able to semi-quantitatively attain BoseEinstein condensate ground state energies even for gases with infinite scattering length. In thelarge particle number limit, our method provides analytical estimates for the energy per particlep
E0 /N ≈ 2.5N 1/3 ~ω and two-body contact C2 /N ≈ 16N 1/6 mω/~ for a Bose gas on resonance,
where ω is the trap frequency. Further, by considering only two-body correlations, we note thata sudden quench from small to large scattering lengths leads to out-of-equilibrium resonant BEC.
As an alternative, we propose a two-step scheme that involves an intermediate scattering length,
between 0 and ∞, which serves to maximize the transfer probability of N bosons in a harmonic
trap with frequency ω to the resonant state. We find that the intermediate scattering length shouldpbe a ≈ 3.16N −2/3 ~/(mω), and that it produces an optimum transition probability of 1.03N −1/6 .

Dedication

To my invisible partner.

v

Acknowledgements

I am grateful to my thesis adviser, John Bohn. Dr. Bohn has trained me to think about
physical problems critically and outside the box, and to tell a story behind every set of data ormodel equation that I gather or derive. I have aspired, with little success, to emulate his exceptional
ability to communicate science in simple pictures. Most importantly, I have learned from him thatphysics is even more fun with humor!
I wish to extend my gratitude to Dr. Doerte Blume. She has been very generous in sharing her
expertise in stochastic methods. “Monte-Carlo method” now sounds less intimidating to me after
she gave a crash course on it. She also gave a handful of invaluable advice and recommendations
to improve my work (part of this thesis).
I thank Dr. Jose D’ Incao for reading this thesis in detail and giving helpful suggestions tomake it more readable. If there is any flaw in this work, I take full responsibility for it. Jose has
also been generous in sharing his knowledge in few body physics. While I was working on the firstpart of this project, he would patiently listen to my updates during group meetings and give hishonest thoughts on it.
I would also like to thank Dr. Victor Gurarie. Through his lectures from some of the graduate
school courses that I took, I learned a few tricks that I was able to apply in this work.
I wish to thank Dr. Eric Cornell and Dr. Deborah Jin (and their group) for working on theexperimental unitary Bose gas. Without their amazing work, my work would seem less relevant.
Special thanks go to John Corson and Andrew Sykes for their direct and indirect input on
and contributions to this work. Their mentorship has made a positive impact on my research. And,vi
thanks to the rest of the Bohn group whom I have had the pleasure to interact with – Eli Halperin,
Lucie Augustovicova, Brandon Ruzic, Goulven Quemener and James Croft.
Thanks to the JILA Computing team – Jim McKown and Corey Keasling – for accommodating all my questions about and problems with linux and the JILA cluster.
To my officemates – Peiru He, Felix Leditzky, and Bjorn Sumner, thanks for making S370 a
conducive work environment.
To all the friends who helped me get through graduate school – Pengxiao Hao, Grace
Grigereit, Winnie Zhuang, Xin Xie, Kia Ng, Hao Wu, Athreya Shankar, Elaine Lim, and Kerim
Kaylan to name a few – by playing badminton, hiking, and fun conversations over good food, orjust by checking up on me occasionally, thank you! I would like to especially acknowledge Javier
Orjuela, Dmitriy Zusin, and Kayla Kohake who were there for me during the most trying times. I
owe you big time.
And lastly, I am forever indebted to my Mommy, Dad, Achi, Dichi and Sachi for their constantlove and support despite my failures.

vii

Contents

Chapter
1 Introduction

1

2 Background8

2.1

Low energy scattering . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.2

Contact interactions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

2.3

Weakly-interacting BEC in a trap

2.4

Review of previous studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16

. . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

3 Toolbox
3.1

3.2

18

Hyperspherical formalism . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
3.1.1

Hyperspherical coordinates . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3.1.2

Hamiltonian and wave function in hyperspherical coordinates . . . . . . . . . 22

Basis functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2.1

Jastrow form and pair wave function . . . . . . . . . . . . . . . . . . . . . . . 25

3.2.2

Boundary conditions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28

3.2.3

Renormalized scattering length . . . . . . . . . . . . . . . . . . . . . . . . . . 31

3.2.4

Hyperangular kinetic energy integral hν|Λ2N −1 |νi . . . . . . . . . . . . . . . . 31

4 Ground-State Properties of the BEC in the Hyperspherical LOCV Approximation
4.1

8

38

General features of the ground state . . . . . . . . . . . . . . . . . . . . . . . . . . . 38

viii
4.2

Positive scattering length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40

4.3

Limiting cases . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.3.1

Small scattering length. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.3.2

Infinite scattering length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5 Two-Step Production of Resonant Bose-Einstein Condensates

50

5.1

Motivation

5.2

The two-step scheme . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.3

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

5.2.1

Franck-Condon factors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

5.2.2

The optimum intermediate state . . . . . . . . . . . . . . . . . . . . . . . . . 60

5.2.3

Wave packet dynamics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

Large N limit . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.3.1

Overlap between hyperradial wave functions using the reflection formula . . . 64

5.3.2

Overlap between LOCV hyperangular wave functions . . . . . . . . . . . . . 66

5.3.3

Transition amplitude and probability . . . . . . . . . . . . . . . . . . . . . . . 69

6 Conclusions and Prospects

73

Bibliography

76

Appendix
A Jastrow-Jackson-Feenberg-LOCV approach to homogeneous Bose gas

85

B Mean Density Calculations at a → ∞88

C Weak Interaction in the Large-N Limit89

D Strong Interaction in the Large-N Limit91

ix
E Saddle Point Approximation

94

x

TablesTable
4.1

Ground state energy per particle, computed by various methods. The uniform gas
energies are given in units of ~2 n2/3 /2m, while the result of Ref. [1] for a trappedsystem is given in units of ~2 hn2/3 i/2m with the mean density hn2/3 i determined byaveraging over a Thomas-Fermi profile in a trap. This same averaging allows thisenergy to be written in terms of the trap frequency ω in the final column.

4.2. . . . . 48

Contact densities, computed by various methods. For the uniform gas, the intensive
contact densities C2 are given in units of n4/3 , while the extensive contact C2 /N , inunits of hn1/3 i, is given in Ref. [1] for a trapped gas. Averaging over a Thomas-Fermi
profile in a trap allows this C2 /N to be written in terms of the characteristic trap
length aho in the final column. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49

xi

FiguresFigure
3.1

The angular wave function φν as a function of r12 for N = 10 at (a) small scatteringplength and on (b) resonance a = ∞. All length scales are in units of aho = ~/mω.
The insets show the zero-crossings of the curves. The crossing occurs at ≈ a forsmall scattering length and at some finite value in the unitary limit. Note that φν
√is actually a function of the hyperangle α. The relation r12 = 2ρ sin α is used toconvert the angle α to the pair distance r12 . Here ρ has its value at the minimum of
the corresponding hyperradial potential.

3.2. . . . . . . . . . . . . . . . . . . . . . . . 32

The zero-crossing ac of the angular wave function φν for N = 10 as a function of
the scattering length a. The inset shows the value ρmin that is used to compute ac ;ρmin is taken to be the value of ρ where the effective potential V eff , discussed in thenext section, is at its minimum.

3.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

The effective potential V eff as a function of the hyperradius ρ for a = 0 (noninteracting case), a > 0 (repulsive interaction), and a < 0 (attractive interaction).
For any a > 0 and N , a local minimum always exists. However, for a given a < 0,there exists a maximum N when the local minimum of V eff starts to disappear. The
highest V eff curve corresponds to the a → +∞ case. . . . . . . . . . . . . . . . . . . 37

4.1

Ground state energy E0 as a function of the scattering length a for N = 10. The
left and right panels consider negative and positive scattering lengths. . . . . . . . . 39

xii
4.2

Ground state energy per particle E0 /N as a function of the scattering length a for
N = 4.

4.3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

Ground state energy per particle E0 /N as a function of the scattering length a for
N = 10. The inset shows a blow-up of the weakly interacting regime. . . . . . . . . . 42

4.4(a) Ground state energy per particle, in units of ~ω. The inset shows the various
perturbative versions (small a region). (b) Two-body contact per particle, in unitsof a−1ho , as a function of the scattering length a for N = 1000. . . . . . . . . . . . . . 45

5.1

The scale of the problem. Each curve represents an effective potential energy surface
for a BEC with a = 0 (bottom) and a = ∞ (top), in our hyperspherical representation. A BEC having a = 0 (Gaussian centered at ρ = 12.2 aho ) has essentially nooverlap with a resonant BEC having a = ∞ (Gaussian centered at ρ = 33.8 aho ). . . 53

5.2

The two-step scheme from non-interaction to small a then to resonance. . . . . . . . 56

5.3

Franck-Condon factors from the (a) non-interacting to intermediate states |h0, 0|a, ni|2 ,(b) intermediate to resonant states |ha, n|∞, 0i|2 , and (c) the two-step transition
probability |h0, 0|a, niha, n|∞, 0i|2 as functions of scattering lengths a and vibrational states n. Here, N = 100.

5.4

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

The two-step transition probability distribution |h0, 0|a, niha, n|∞, 0i|2 as a functionof scattering lengths a and vibrational states n for N = 1000. . . . . . . . . . . . . . 59

5.5Transfer probability for (a) N = 100 with a∗ = 0.0859aho , and (b) N = 1000 witha∗ = 0.0332aho . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

5.6

Mean radius of the BEC in the intermediate phase versus time before quench tounitarity for N = 100 and a = 0.0859aho . . . . . . . . . . . . . . . . . . . . . . . . . 63

5.7

Transfer probability of the BEC versus scattering length a for large N . Inset showsa zoom-in profile of N = 105 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70

5.8

Transfer probabilities for N = 1000 starting from non-interacting and weakly interacting initial states. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72

Chapter 1

Introduction

It is universally acknowledged that all matter is made up of particles called atoms, which havean individual size on the order of an angstrom. These atoms may be freely moving in space independent of other atoms or external forces, or interacting with each other in the presence or absenceof some external potential. Strong interactions increase the probability of two-body collisions while
three atoms may also collide inelastically and form a molecule. However they interact (or not),they obey the laws of quantum mechanics. The distinct nature of these particles and their interactions induce some fascinating physical phenomena, among which is Bose-Einstein condensation.
Bose-Einstein condensation is a phase transition in dilute gas associated with the condensation of
atoms in the lowest quantum state which occurs under extremely low temperature conditions.
This phenomenon owes its conceptualization to Satyendra Nath Bose and Albert Einstein.
The formation of Bose-Einstein condensates (BECs) was predicted by Einstein [2, 3] based on
Bose’s [4] statistical description of light quanta in deriving Planck’s law for the black-body radiation.
Ideas from Bose and Einstein led to the identification of indistinguishable particles called bosons. A
defining feature of bosons is that they can all occupy the same quantum state. Einstein later argued
that below a certain ultracold temperature, as the kinetic energy of these particles diminishes, theywould condense to the lowest energy state. Seven decades after its theoretical prediction, BECs of
rubidium [5], sodium [6], and lithium [7] vapors were created in separate experiments performed in
JILA, MIT, and Rice University. These alkali atoms were trapped magnetically and cooled downto ∼ 100 nK temperatures. Images of the atoms’ velocity distribution with varying temperatures2
were concrete evidence of BEC. Below a certain critical temperature, a sharp peak around thezero-point of the velocity distribution was observed.
The success of the first series of BEC observations has pushed further investigations of quantum correlations in ultracold gases. When atoms in ultracold Bose gas are brought close to absolute
zero temperatures, their individual wave functions start to overlap with each other. The groundstate wave function of the non-interacting bosons is expressed as a product of identical lowest
one-body state wave functions. Note that this form of many-body wave function preserves thesymmetry property of bosons; that is, exchanging two particles leaves the many-body wave function unchanged. Further, such a many-body wave function results in a set of independent andidentical single-particle wave equations.1 The many-body problem is then reduced to a collectionof one-body problems. The non-interacting system is a standard statistical mechanics textbook
problem whose solution is well-understood [8, 9]. Nontrivial many-body effects come into play inthe presence of interactions as the wave function of a particle is affected by those of the others’.2
A basic approach to a many-body problem starts from understanding two-body interactions.
A two-body problem can generally be reduced to two effective one-body problems involving
the center of mass and the relative motions. The center of mass motion is described by its kinetic
energy (and the trapping potential if it exists). On the other hand, the relative motion, in thepresence of interaction between atoms, has an additional (interaction) potential V (r) where r is
the relative coordinate. Formally, interactions can be classified as either short- or long-ranged.
Short-ranged potentials decay faster than r−1 ,3 while long-range potentials do not decay as fastand have long tails. For short-ranged potentials, one can define the range of interaction r0 that1

If Ĥ = Ĥ1 + Ĥ2 + ... + ĤN and Ψ =QN

i=1 ψ0 (xi ) are the N -body Hamiltonian and wave function of a non-interacting system, then we get N identical set of equations Ĥi ψ(xi ) = Eψ(xi ), where E is an eigenenergy of Ĥi .
2This can be understood by incorporating an interaction potential V int in the many-body Hamiltonian Ĥ andarriving at the conclusion that one can not simply get identical and independent wave equations in the form Ĥi ψ(xi ) =Eψ(xi ) as with the non-interacting case. The Gross-Pitaevskii equation is a good illustration; the nonlinear term
describes how the other atoms affect one.
3

An exponential/Gaussian potential or a square-barrier is also considered short-ranged.

3
indicates the length scale below which the interaction takes place. For example, with a squarewell/barrier potential, the width of the well/barrier is the r0 . With a short range interaction, atom
1 can ‘see’ atom 2 only when the mean separation between the two atoms is smaller than therange of interaction. A quantum scattering event can be understood in these terms: at low-energy
scattering, a slow moving (low momentum) particle approaching a short-ranged scatterer has ade Broglie wavelength much longer than r0 . Consequently, the structure of the potential becomesinsignificant because it cannot be resolved by the incoming particle. At distances much longer thanr0 , the radial wave function of the scattering particle is, approximately, a phase-shifted sinusoidal
wave where the phase shift is a measure of how far the scattering particle’s origin is displaced(from a free particle’s origin). In the absence of an interaction potential, the phase shift is zero.
The phase shift can be expressed in terms of the low energy s-wave scattering length a whichcharacterizes the strength of interaction. The scattering length can be either positive or negative.
Positive a indicates repulsive interaction while negative a indicates attractive interaction because
the long-range wave function is pushed out or shifts to a larger r when a > 0, and is pulled in or
shifts to a smaller r when a < 0.
Early studies on BECs investigated the small a or weak interaction regime. The smallnessof a is relative to the mean interparticle spacing which is on the order of n−1/3 , with n as thenumber density of the gas. In particular, the diluteness condition n|a|3  1 is the basis forconsidering a system to be weakly interacting.4 With this condition satisfied, it is found that amean-field approximation [11, 12, 13, 14, 15, 16, 17] provides an adequate description of the weaklyinteracting system which is easily verified in a typical experimental set-up [18, 19, 20]. In thismodel, the many-body wave function, like the one for non-interacting system, is a product state ofidentical one-body wave functions. But the single-particle wave function is now a solution to theGross-Pitaevskii equation [13, 14, 15] which is a nonlinear second-order differential equation. The
4

Interestingly, Ref. [10] points out that this condition does not measure the weakness of interaction relative tothe kinetic energy. An interacting system may very well satisfy the diluteness condition, but can have an interactionenergy that exceeds the kinetic energy.

4mean-field energy, E, of the interacting Bose gas is, to first order, given by
E = 4π

~2 a
n,
m

(1.1)

where m is the mass of the boson. Beyond mean-field is the perturbative Lee-Huang-Yang (LHY)
correction where the small parameter is na3 [16, 17]. Including this correction, the energy now is~2 a128 √ 3E = 4πna .
(1.2)
n 1+ √m
15 πThis correction has also been verified experimentally [21, 22, 23, 24].
Recent experimental progress has pushed the physics of dilute BECs into the regime ofinfinite a by means of tuning magnetic fields near a Fano-Feshbach resonance [25]. These stronglyinteracting gases are also described as resonant or unitary. The latter term refers to the regimeof largest possible scattering cross-section σ, which is related to a via σ ∼ a2 ; with a maximal σ,scattering probability approaches unity. The study of resonant Bose gas poses a major challenge forboth theorists and experimentalists. Strong atomic losses are observed in the laboratory making
the resonant Bose gas difficult to measure accurately. In particular, three-body recombination
is a loss process responsible for molecule formation and subsequent destabilization, which occurs
within milliseconds. Before the catastrophe sets in, however, a resonant Bose gas exhibits a quasisteady-state behaviour within 100 µs [26], which is inferred from the time-resolved momentum
distributions. Formation of bound state dimers and trimers of atoms was also observed [27, 28,
29, 30, 31, 32]. These are rare events in dilute ultracold gases because the probability of three- or
higher body collisions is low. A typical unitary Bose gas experimental protocol starts with the gas
at a small value of scattering length, then quickly ramps or ‘quenches’ the magnetic field, which
controls the scattering length, to a Fano-Feshbach resonance within microseconds. This represents
an instantaneous projection of the many-body state at small a onto a collection of many-bodystates at a = ∞. The quench reduces the fast atom number decay from thermal heating if theprocess was done adiabatically [33].
On the theoretical side, the application of mean-field approximation leads to meaninglessresults as a diverges and na3  1 does not hold in this regime. An immediate consequence of using

5
a mean-field picture is obtaining an infinite characteristic energy scale. Because the scattering
length a becomes a physically irrelevant length scale at unitarity, length scales are expressed interms of the mean interparticle spacing, hn−1/3 i, or equivalently the size of the gas confinement.
Correspondingly, by dimensional analysis, the associated energies scale as ~2 hn2/3 i/m, a finite
quantity. Although the interaction diverges, it is physically absurd to have an infinite energy which
implies an infinite expansion of the ultracold gas. Theories on unitary Bose gas abound. They
all agree with the universal ground state energy scaling ~2 n2/3 /m. The quantity in front of thisscaling, however, remains in dispute. Field-theoretic approaches and other formalisms replace thedivergent interaction potential by an effective one that is finite, which is done via renormalizationof a or an equivalent interaction parameter [34, 35, 36, 37, 38, 1, 39, 40, 41, 42, 43].
In this work, we present an alternative approach to constructing a finite effective potential– one that does not utilize explicit renormalization. A first step is to define explicit approximate
wave function solutions for the resonant BEC. The wave functions are expressed in terms of acollective coordinate, the hyperradius, which describes the motion of the condensate as a whole,along with two-body interparticle coordinates [44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 1]. Clearly, thismethod differs from the mean-field formalism where single particle wave functions are employed.
The single-particle and two-body wave function methods, however, yield comparable results inthe weak interaction regime. Inspired by the Born-Oppenheimer formalism that is often appliedin molecular physics, we treat the hyperradius as the slow coordinate and solve the two-body
motion first. We adopt the product state of two-body wave functions – an ansatz first conceived
by Jastrow [54]. Moreover, instead of including an explicit two-body potential in the Hamiltonian,
the Bethe-Peierls boundary condition is employed with actual scattering length. This boundary
condition describes what happens when two atoms are close to or ‘in contact with’ each other.
Another boundary condition arises from interpreting the pair particle wave function as a two-bodycorrelation function – a concept developed in the lowest-order constrained variational (LOCV)
method [55, 56]. This condition is set by a critical interparticle spacing with which a pair ofatoms becomes uncorrelated. After determining the two-body motion, an effective potential in the6collective motion coordinate, which consists of the centrifugal component, the trapping potential
and the two-body interaction potential, can then be drawn for a given a. Breathing modes of thecondensate can also be extracted and studied. This potential saturates as a → ∞ which allows us
to identify the ground state of a resonant BEC.
Three-body physics [34, 57, 58, 59], to which atom losses are attributed, becomes importantin understanding unitary Bose gas as evident in the formation of of trimers. These three-bodystates are beyond the scope of our model which considers only two-body correlations. Our methodof approach, however, is able to illustrate the dearth of resonant BEC that results from a direct
projection to unitarity from a weak interaction. That is, if Ψ0 and Ψf are the wave functions
representing the initial weakly interacting and the final resonant BECs, respectively, then we findhΨ0 |Ψf i ∼ 0. A two-step process is then proposed to maximize the amount of condensate that
arrives in the resonant state. In the first step, the scattering length is jumped quickly from a low
initial value a1 ≈ 0 to a modest intermediate value a2 . Expansion takes place due to the sudden
increase in scattering length. When it reaches the size of the resonant condensate, the scattering
length is suddenly jumped from a2 to a → ∞. A properly chosen a2 yields a non-negligible fractionof atoms which are converted into a resonant BEC.5 We aim to seek the intermediate state, subject
to this a2 , yielding a maximum transfer. The theory focuses on an isotropic, harmonically trapped
Bose gas, which is one of the widely-used experimental set-ups [61, 24, 26, 31, 30]. The breathingmodes supported by the effective potential for a given a2 allow us to express the intermediate BEC
state as a wave packet which propagates in time according to the linear Schrödinger equation. The
procedure involved is reminiscent of the optical and Raman transitions [62] and vibrational wave
packet dynamics in molecular physics [63].
A notable feature of the techniques applied in this work is the resulting set of approximationsin asymptotic limits. These simple estimates may serve as a guide for interested parties seeking
to validate the efficiency of the two-step scheme. It is a goal of this work to stimulate further
5Interestingly, the opposite technique was applied in the JILA experiment, jumping to a smaller scattering length

to create a more dense condensate, to explore density effects [60].

7theoretical and experimental probing of various intriguing and unexplored aspects of the unitary
regime.
This thesis is organized as follows: Chapter 2 reviews some of the familiar and relevant
quantum mechanics concepts such as low energy scattering theory and two-body interactions. This
chapter also gives a brief discussion of the progress on the theoretical studies of BEC. Chapter 3
outlines our method of approach to understanding BEC which combines hyperspherical formalism
and the lowest order constraint variational method. Chapter 4 presents the ground-state propertiesof the weakly- to strongly-interacting BECs extracted from the hyperspherical-LOCV method.
Chapter 5 introduces a novel protocol that may improve the production of resonant BECs.

Chapter 2

Background

Early progress towards understanding a system of ultracold atomic gas exploited the dilutenature of the gas1 so that the effects of the interactions can be described in terms of the twobody scattering length alone. In other words, two interatomic potentials corresponding to thesame scattering length lead to the same properties for the condensed gas, although they may
have completely different microscopic properties. This then led to the formulation of the meanfield theory of weakly interacting bosons (na3  1) where the scattering length is a universal
parameter. Our focus will be on interacting bosons in a isotropically harmonic trap for which theGross-Pitaevskii theory most suitably applies for weak interactions.
In this chapter, we review some of the relevant aspects of low-energy two-body scattering
theory and a few short-ranged atomic interactions. We then discuss briefly the results of mean-field
theory. Lastly, we discuss some of the previous studies on unitary Bose gas when the scattering
length becomes an irrelevant parameter.

2.1

Low energy scattering
In this section, we follow the lines of discussion in Ref. [64].
Elastic scattering of two slow bosons can be described by the radial Schrödinger equation


1


~2 d2
−
+ V (r) − Erel u(r) = 0,
m dr2

(2.1)

The average interatomic spacing is much larger than the range of interatomic forces. That is, na3  1 andnr03  1.

9where r = |r1 −r2 |, V (r) is the two-body interaction potential, m is the particle mass, u(r) = rψ(r)the s-wave2 radial function, and Erel is the relative energy. In the asymptotic region r  r0 , wherer0 is the range of the potential V , the solution is given byu ∼ sin (kr + δ(k)) ,where k =

√

(2.2)

mErel /~ and δ(k) is an incident energy-dependent phase-shift in the presence ofinteraction. Without interactions, δ(k) = 0. Here, we have assumed that only the lowest orderterm, which corresponds to ` = 0 angular momentum, contributes in the partial wave expansion of
the scattering amplitude. The scattering amplitude f (k), defined as the ratio of outgoing sphericalwaves and incoming plane waves, is in this case spherically symmetric and related to the phaseshift byf (k) =

1.
k cot δ(k) − ik

(2.3)

The total cross-section σ, which is an experimentally measurable quantity, is related to f (k) via
σ = 8π|f (k)|2 for identical bosons. In the low energy limit,
lim σ`=0 (k) = 8πa2 ,

k→0

(2.4)

where the scattering length a is defined by
a = − lim f (k) = − lim
k→0

k→0

tan δ(k)
.
k

(2.5)

From this relation, one can see that for k|a|  1, δ(k) ≈ −ka. A positive (negative) a corresponds
to a wave shifted to larger (smaller) particle separation with repulsive (attractive) potentials. By
expanding the u(r) in Eq. (2.2) up to first order in k,u(r) ∼ kr cos δ + sin δ,
k→0

(2.6)

one concludes that the scattering length is the zero-crossing of the wave function as k → 0. Thus,in the limit of low-energy and region where r  r0 , the wave function takes the linear formru(r) ∼ 1 − .
k→0
a
2

(2.7)

The “s” in s-wave corresponds to the lowest angular momentum state with ` = 0. Recall “s-p-d-f” from chemistry.

10
An illustrative example is the square barrier potential given byV0 if r ≤ r0
V (r) =

0
if r > r0 ,where V0 > 0 and r0 are the height and width of the barrier. In the low momentum or energy limit
k → 0, the solution to Eq. (2.1) given this potential isA sinh (k0 r) if r ≤ r0
u(r) =B (r − a)if r > r0 ,√where A and B are normalization constants, and k0 =mV0 /~. Continuity of u(r) and u0 (r) at

r = r0 implies
a = r0 −tanh(k0 r0 ),k0

(soft sphere).

(2.8)

We note that a > 0 always. For infinitely tall barrier V0 → ∞, we recover the hard sphere solution
with the scattering length equal to the hard sphere core radius a = r0 .
If V0 < 0, thena = r0 −where k00 =tan(k00 r0 )
,.
k00

(2.9)

p
m|V0 |/~. In this case −∞ < a < +∞. The divergence of a, also known as potential

resonance or shape resonance[65], occurs when k00 = (q + 21 )π, where q is an integer or when thepotential is deep enough to support a new bound state. Two-body bound states (dimers) haveenergies −|V0 | < Eb < 0. If two colliding atoms form a bound state with an energy just below thecontinuum (unbound) threshold and |Eb |  |V0 |, then this state is known as a weakly bound stateor shallow dimer (or Feshbach molecule) and has a binding energy given by
Eb = −~2.
ma2(2.10)

One can derive this by finding the pole of the scattering amplitude in Eq. (2.3) and using Eq. (2.5).
Real experiments exploit internal atomic structure (hyperfine and Zeeman) to control ortune atomic interactions. So far, we have ignored the internal structure of the atoms in a scattering

11
process. Doing so, we are describing the process as single channel in which the Hilbert spacespanned by the states before and after collisions (entrance and exit channels) is the same. Now,
considering the hyperfine structure of the atoms like the spin introduces more degrees of freedom
and a change of internal states of the scattering atoms, which can result in multichannel scattering.
For two atoms, if the total energy is given by Etot = Erel + Eα , where Eα is the internal energy of
the individual atoms which is also the channel energy with which the atoms are prepared. Then
any channel with energy less than Etot is considered as open channel and any channel with energy
greater than Etot is called a closed channel from which the atoms cannot escape because theirenergy is too low [25, 66]. Feshbach resonance occurs when the bound state energy from a closed
channel potential is close to the (zero) energy of an unbound state from a different open scatteringchannel. The relative position of the channels can be changed by an external magnetic field, thereby
tuning the coupling. If B0 is the magnetic field at which Feshbach resonance takes place, then the
scattering length a as a function of the magnetic field B is given by [67]∆a(B) = abg 1 −,B − B0(2.11)

where abg is the background scattering length, ∆ is the width of the resonance.

2.2

Contact interactions3
From Eq. (2.9), we infer that there exist two distinct pairs of parameters {r0 , V0 } which canyield the same a. That is, we can make the width r0 as small as possible by increasing V0 to fix a.
This is also reflected if we consider the contact or zero-range interaction potential
V (r) = gδ 3 (r),(2.12)

which turns on only when two particles are in contact or the separation between them is zero, and3δ 3 (r) is the Dirac-delta function in 3-dimensional space. One can fix g = 4π
3 r0 V0 by letting r0 → 0

and V0 → ∞. However, if we go back to Eq. (2.9) and evaluate a in this set of limits, we get a = 0.
Thus, Eq. (2.12) is not a good model for contact interaction.
3

Again, we follow the lines of discussion in Ref. [64].

12Contact interaction is properly described by the Fermi pseudo-potential [68, 16, 69]
V (r)ψ(r) = gδ 3 (r)∂∂u(r)(rψ(r)) = gδ 3 (r).
∂r∂r(2.13)

Thus, if rψ (or u) is analytic or regular at r = 0, then V (r)ψ(r) = gψ(0)δ 3 (r) (or V (r)ψ(r) =3g ∂u∂r r=0 δ (r)), which is the desired effect of the contact interaction potential. The relative wave

function ψ satisfies
−

~2 2
∂
∇ ψ + gδ 3 (r) (rψ(r)) = Erel ψ.
m∂r(2.14)In spherical coordinates,∇2 ψ =

1 ∂2u− 4πu(0)δ 3 (r),r ∂r2

(2.15)

where we used the identity ∇2 1r = −4πδ 3 (r) and considered only ` = 0 so that the angular term
vanishes. Rewriting Eq. (2.14) in terms of u, we get
−

~2
m

 00

u (r)
~2 k 2 u(r)− 4πu(0)δ 3 (r) + gδ 3 (r)u0 (0) =.
r
m r(2.16)We separate the Dirac-delta terms and the regular terms, and get a second order differential wave
equation and a boundary condition:
u00 (r) = −k 2 u(r),u(0) =

mg 0
u (0).
4π~2(2.17)(2.18)

Using the definition of a in Eq. (2.5), the general solution to Eq. (2.17) is given byu(r) = A [sin(kr) − ka cos(kr)] ,(2.19)

and findg=4π~2 a.
m(2.20)The boundary condition in Eq. (2.18) can be rewritten as
1 ∂u
1
=− ,
u ∂r r=0a(2.21)

13
which is also known as the Bethe-Peierls boundary condition [70]. Note that the linear low-energy
wave result in Eq. (2.7) satisfies this boundary condition. For the zero-range interaction, theshort-range behaviour of u intercepts the r-axis at r = a. The Fermi pseudo-potential method is
instrumental in deriving the LHY and Wu corrections to the mean-field ground state energy of a
weakly interacting BEC [16, 17, 71]. The Bethe-Peierls boundary condition is equally useful andwill be applied in this work instead of using the pseudo-potential method.

Weakly-interacting BEC in a trap4

2.3

The mean-field approach has proven to provide accurate quantitative predictions for theproperties of a trapped dilute Bose gas [10]. The diluteness of atomic BEC allows one to describe the system in an independent particle picture which is similar to a system of trapped noninteracting BEC. Bogoliubov developed a single-particle-wave-function formulation by exploiting
the Fock space (second quantization) and utilizing bosonic field operators [11]. As a result, anonlinear Schrödinger equation emerges:∂i~ Ψ(r0 , t) =
∂t



~2 2
2
0
ext 0
Ψ(r0 , t),−∇ + V (r ) + g Ψ(r , t)2m

(2.22)

which is also known as the Gross-Pitaevskii (GP) equation [13, 14, 15]. Here, Ψ(r0 , t) is thecondensate wave function at position r0 and time t, V ext is the trap potential, and g = 4π~2 a/m
is the same g in Eq. (2.20). The condensate wave function is normalized to the total number ofbosons N ,Z4dr0 Ψ(r0 , t)2

References [10, 67] are the main sources for this section.

= N0 ≈ N,(2.23)

14
where N0 is the number of condensed atoms.5 The time-independent mean-field energy is given byZ
E [Ψ] =

dr0




~2
g
|∇Ψ|2 + V ext (r0 )|Ψ|2 + |Ψ|4 ,2m
2(2.24)

where the first and second terms yield the total energy of the non-interacting BEC in a trap, andthe last term is the interaction energy E int . Writing the wave function asΨ(r0 , t) = ψ(r0 ) exp (−iµt/~),(2.25)

where µ is the chemical potential, the particle density n is defined byn(r0 ) = ψ 2 (r0 ).

(2.26)

Thus, the interaction energy in Eq. (2.24) scales as E int ∼ gN hni where hni is the average density6 .
Also, the time-independent Gross-Pitaevskii Eq. (2.22) becomes

−


~2 2
2
ψ(r0 ) = µψ(r0 ).
∇ + V ext (r0 ) + g ψ(r0 )2m(2.27)For a spherically-symmetric harmonic trap, V ext = mω 2 r02 /2, where ω is the trap frequency.
p
The harmonic oscillator length is aho = ~/(mω), and the average density hni ∼ N/a3ho so that
E int ∼ N 2 |a| /a3ho . Expressing all the length and energy quantities in terms of the aho and energy~ω, the rescaled time-independent GP equation is given by

5

˜ 2 + r̃2 + 8π
−∇



Na
aho




Ψ̃2 (r̃) Ψ̃(r̃) = 2µ̃Ψ̃(r̃),(2.28)

An equivalent version of Eq. (2.22) is given byi~∂Ψ(r0 , t) =
∂t



~2 2
2
−
∇ + V ext (r0 ) + g (N − 1) Ψ(r0 , t)Ψ(r0 , t),
2mwith normalizationZ

dr0 Ψ(r0 , t)

2

= 1.

A more suggestive interpretation can be extracted from this form of GP equation: Ψ is the wave function of arepresentative atom, and the nonlinear term in the wave equation describes the effect of the rest of the atoms has on
the representative atom, which is proportional to their density (N − 1) |Ψ|2 [45]. The N − 1 factor, instead of just N
will matter if N < 102 .
RR2
4
dr0 n|ψ(r0 )|dr0 |ψ(r0 )|6
hni = R dr0 |ψ(r0 )|2 = R dr0 |ψ(r0 )|215qwhere Ψ̃ = Ψ/ N a3ho . If N a/aho  1, the first term (kinetic) in Eq. (2.28) can be ignored, so thatv2uu µ̃ − r̃2
t


Ψ̃ =
4π aNhoa(2.29)and the Thomas-Fermi (TF) approximation can be applied. Consequently, we find1µ̃T F =
2



15N a
aho

2/5
.

(2.30)

Using the thermodynamic relation µ = ∂E/∂N , the total energy per particle in the TF limit is
ET F
5= µT F .
N
7(2.31)We also note that for Ψ̃ (or the density n) in Eq. (2.29) to be real,
RT F
r̃ < R̃T F =
aho



15N a
aho

1/5
,

(2.32)

where RT F is known as the Thomas-Fermi radius. Graphically, RT F is the point where the density
profile Ψ̃2 , in the form of an inverted parabola, of the BEC goes to zero.
Local density approximation is widely used for large N . The chemical potential is estimatedby
µ = µloc n(r0 ) + V ext (r0 ),(2.33)

where µloc = gn. Including the LHY correction,
32 √ 3µloc (n) = gn 1 + √
na .
3 π


(2.34)

Solving Eq. (2.33) iteratively, we findp
µ = µT F 1 + πa3 n(0) ,


7p 3
5
πa n(0) ,E = N µT F 1 +
7
8

12/5
152/5
3
1/6 a
a n(0) =
N
.
8π
aho

(2.35)(2.36)(2.37)

Monte Carlo simulations of trapped BECs have been applied and the numerical results are foundto agree with the above GP’s with a few percentage difference [72, 73].

16

2.4

Review of previous studies
The boundlessness of the universal parameter a does not deter the pursuit of understanding

the unitary BEC. A plethora of theories have been presented. While the methodologies vary, acommon theme among the existing studies revolves around renormalizing the diverging interaction
potential or equivalently the scattering length. Many of these approaches can be mathematicallysophisticated and rigorous. Predictions for the observables, such as the chemical potential energy,
speed of sound, momentum distribution, and two-body contact7 , have been derived or extracted.
In considering a uniform unitary BEC, a general direction starts by writing the Hamiltonian
in terms of the bosonic field operators, similar to the Bogoliubov method but with the nonzero
momentum (and bound) states playing significant roles because of the expected depletion of condensate at unitarity. Song et.al. [43] constructed a square well potential for the interaction, andthe condensate (k = 0) wave amplitude in their ground state ansatz was varied to minimize the
Hamiltonian subject to some conditions related to the vacuum of the Bogoliubov quasiparticles.
They also slightly touched on BEC in a trap and found the condensate radius to be ∝ N 1/6 aho .
A number of groups derived an effective interaction potential by introducing momentum cutoffs[36, 35, 38, 40], or by expanding the coupling constant g in terms of the size of the two-bodybound state [77, 37]. Cowell et.al. [78] exploited a spatial wave function method in the form of the
pairwise Jastrow wave function and worked within the lowest order regime where only two-bodycorrelations would matter. Diederix et.al. [79, 66] presented a rigorous extension of the lowest
order constraint variational method by doing the hypernetted-chain approximation. Rossi et.al.
[80] employed a Monte Carlo simulation with a two-body square-well potential using also a Jastrow
ansatz. A universal conclusion from all these studies is that the chemical potential of a resonant
BEC is expressed as µ=a→∞#F , where F = 6π 4/3 ~2 n2/3 /(2m) is the Fermi energy, and thenumber # differs by about a factor of four in different studies (See Table 4.1).
7

Contact is a quantity that determines how the energy changes with scattering length. It describes the short-

ranged behaviour of pairs of atoms [74, 75, 76].

17
While there are a handful of sensible and conclusive theories on properties of trapped systemsof few bosons, weakly to strongly interacting, analytical studies on large-N BEC in a trap seems
sparse.8 For two atoms in a harmonic trap, the s-wave energy spectrum is given by [81]a1 Γ [1/4 − E/(2~ω)]=√,aho
2 Γ [3/4 − E/(2~ω)](2.38)

where E = ~ω (2ν + 1/2) with ν = 0, 1, 2, ... if a → ∞. For a three-body system, a Fadeev wave
function ansatz (sum of two-body wave functions) was utilized to obtain a ground state relative
energy of Erel ≈ 5.47~ω as a → ∞ [82], which agrees well with the Erel ≈ 5.52~ω result of a Monte
Carlo procedure with a Gaussian interaction potential [83]. Universal three-body bound states (ortrimers) whose energies are also dependent on the trapping frequency have been identified [82, 84].
For four bosons, the Monte Carlo technique in [83] was applied to produce Erel ≈ 8.33~ω. Thøgersen
et.al. [85, 86, 87] employed a stochastic variational method with correlated Gaussians as basis tocalculate the energy and the condensate fraction of a finite-size system (N < 100) with repulsive,
zero-range, and attractive model potentials. They found that only the attractive model potential
gave converged, independent of a, results in the large a regime. Interesting analytical properties,
such as the 2ω-spacing of the breathing mode frequencies, were derived in Ref. [88] for the N -atoms
in a trap with infinite scattering length. A combination of density-dependent renormalization ofa and mean-field model was employed in Ref. [1] to obtain the scalings of the energy per particle(∼ ~2 hn2/3 i/(2m)) and contact per particle (∼ hn1/3 i) of a unitary Bose gas. The 2ω-spacing of
the breathing mode frequencies was also verified numerically in Ref. [1].
Many of the tools and techniques, such as the hyperspherical method, the Jastrow ansatz, andthe lowest order constraint variation,9 used in some of these studies will be adapted in our formalism
to give our own description of unitary BEC. Thus, similarities and differences are expected among
the existing and our own versions.

8

At least, at this time of writing.

9

These will be discussed in more detail in the next chapter.

Chapter 3
Toolbox1

An objective of this work is to solve the Schrödinger equation for the N -particle system with
Hamiltonian
H=

N  2
X
p

1
+ mω 2 ri2
2m 2
i

i=1


+

N
X

V (rij ),

(3.1)

i<j

for a collection of identical bosons of mass m interacting via pairwise potentials V , and confined toa spherically symmetric harmonic-oscillator trap with angular frequency ω. Invoking the dilutenessof an ultracold Bose gas under typical experimental conditions, where the range of the potential
V is far smaller than the mean spacing between atoms, we may consider the atoms as mostly
independent particles, with Hamiltonian
H=

N  2
X
p


1
2 2
+ mω ri ,
2m 2
i

i=1

(3.2)

and with two-body interactions represented by the zero-range two-body boundary conditions –the Bethe-Peierls conditions2 , which describe the wave function when two atoms are touching one
another. That is, if the total wave function is denoted by Ψ, then we requirelim

1

rij →0 (rij Ψ)∂(rij Ψ)
1
=−
∂rij
a

(3.3)

for any pair distance rij , where a is the two-body scattering length. More generally, three-bodyboundary conditions are also required to describe the gas, but we do not consider them in thiswork.
1

The formalism and results presented in this chapter were published in Ref. [89].

2

See Eq. 2.21.

19

3.1

Hyperspherical formalism
Our method of attack starts by choosing a small set of relevant coordinates to the describe

the system. This includes the hyperradius and relative coordinates between pairs of atoms. The
hyperradius is an “overall” coordinate describing the size of the condensate and its collective breathing modes. The two-body relative coordinates provide convenience for implementing the boundary
condition (3.3) and accounting for correlations of atom pairs.
It should be noted that the hyperspherical formalism is not an uncommon choice for describing
an intractable many-body problem, for example, in chemistry, nuclear physics, and in the study ofatomic structure [90, 91, 92, 44, 93, 45, 94, 95, 46, 48, 49, 96, 51, 97, 98, 52, 1]. The details anddirections of how it is utilized, however quite vary. The potential harmonic expansion method, where
an explicit short-ranged two-body potential (in the form of Gaussians, for example) is expanded
in terms of the hyperspherical potential basis, seems to work well for small number (N < 100)of particles and weak interactions [99, 90, 98, 50, 52]. The K-harmonic method, where the wave
function is expressed in terms of a single harmonic basis, was explored for contact interactions fromweak [45] to resonant [1] regimes. One other route to where this hyperspherical formalism was taken
is the use of Faddeev [100, 101, 102] model, which expresses the hyperangular wave function as asum of pairwise wave functions [46, 103, 104, 59, 47, 49, 105]. This model, with simplified method ofsymmetrization as outlined in Ref. [46], is sufficient and accurate in describing three-body systems(as it was originally designed for) from weakly interacting to unitarity regimes [104, 59, 103, 82, 84],
but the accuracy3 in going above N = 3 seems to be adequate only within the small interactions
as we will show later. In all of these, the application of hyperspherical adiabatic potentials, which
we will discuss later, appears to be a natural way to proceed, taking advantage of the adiabatic
separation4 between the hyperradial and hyperangular motions.
3

Note that we are referring to the Faddeev model in Refs. [46, 47, 49, 105] where it has a unique and generalized

way of approximating the symmetrization because the exact Faddeev wave function would be numerically intractable
to construct as N becomes larger. The number of terms in the exact Faddeev wave function increases by N (N − 1)/2.
4

That is, one of the coordinates can be treated as slowly moving relative to the other.

20
3.1.1

Hyperspherical coordinates

Here, we follow the coordinate system and notation of Ref. [46], defining first the center ofmass coordinate
N

Rcm =

N

1 X
1 X
mri =
ri .
Nm
N
i=1

(3.4)

i=1

The remaining relative coordinates of the particles are conveniently described by a set of N − 1
Jacobi vectors that locate each atom k from the center of mass of the preceding k − 1 atoms:


r
N
−k
X
N −k 
1
ηk =rN −k+1 −rj  .
(3.5)
N −k+1N −kj=1These Jacobi vectors form Cartesian coordinates in a d ≡ 3(N − 1)-dimensional configuration
space. The collective coordinate is the hyperradius ρ, the radial coordinate in this space, given by
the root-mean-squared interparticle spacing of the configuration of atoms:
2

ρ =

N
−1
X
k=1

N

ηk2 =

1 X 2
rij .
N

(3.6)

i<j

Thus the hyperradius is a measure of extent of the system.
The angular coordinates on this hypersphere, collectively denoted by Ω, may be chosen in agreat many different ways [106, 44, 107]. A main point, however, is that all such angular coordinates
are bounded and therefore eigenstates of kinetic energy operators expressed in these coordinateshave discrete spectra and are characterized by a collection of as many as 3N − 4 discrete quantum
numbers. These are complemented by a description of the motion in ρ, which is also boundedwithin a finite range due to confinement by the harmonic-oscillator potential. The relative wave
function of the system can therefore be expanded in a discrete basis set, whose quantum numbers
describe the modes contributing to an energy eigenstate, or else the modes excited in a dynamical
time evolution of the system.
We adopt the hyperangular coordinates in Ref. [46]. In total, there are 3N − 4 hyperangles
in Ω. 2(N − 1) of the hyperangles are simply the angular coordinates (in real space) of the Jacobi
vectors, (θk , ϕk ). The remaining hyperangles αk are angles of radial correlations among the Jacobi

21
coordinates,
ηk
P

k
2
l=1 ηl1/2 = sin αk ,k = 2, 3, ..., N − 1.

(3.7)

In this definition α1 = π/2 is not a separate coordinate [46], i.e., if k = 1 then α1 = π/2 is fixed by
Eq. (3.7). Among these αk , we single out the angle that parametrizes a single pair,sin α ≡ sin αN −1 =

ηN −1
r12
=√ .
ρ
2ρ

(3.8)

The hyperangle α ranges from 0 to π/2. The lower bound α = 0 describes a situation whereparticles 1 and 2 precisely coincide; the upper bound α = π/2 is the point where many-body effects
most likely prevail. For our present purposes, we assume that the orbital angular momentum ineach relative coordinate is zero, whereby the angles (θk , ϕk ) are irrelevant. We are also goingto restrict our wave functions to those with an explicit dependence on r12 only – a choice thatentails symmetrizing the wave function among all pairs rij . With this assertion, the wave function
calculations will be carried out in the coordinates (ρ, α), the two “most relevant” degrees of freedom.
In this coordinate system, the volume element of the relative component is [46, 91]
ρ3N −4 dρ dΩN −1 ,(3.9)

where the solid angle element on the hypersphere is defined recursively asdΩk = sin θk dθk dϕk sin2 αk cos3k−4 αk dαk dΩk−1 .

(3.10)The most important volume element for our purposes is that for the hyperangle α = αN −1 , which
we express in the specialized notationdΩα = sin2 α cos3N −7 α dα = J(α)dα(3.11)

that defines a shorthand notation for the Jacobian J(α). We also single out the angular coordinates(θ, ϕ) = (θN −1 , ϕN −1 ) of this Jacobi vector,dΩN −1 = sin θdθ dϕ J(α)dα dΩN −2 .

(3.12)

22
3.1.2

Hamiltonian and wave function in hyperspherical coordinates

In the coordinates described in the previous section, we can write the Hamiltonian in Eq. (3.2)
as H = Hcm + Hrel , with the center of mass Hamiltonian Hcm given by
Hcm =where pcm =p2cm 1
2+ M ω 2 Rcm,
2M2(3.13)

PN

i=1 pi and M = N m are the total momentum and mass, and the relative motion

component is [94, 46]"#
Λ2N −1
1∂ 3N −4 ∂1~2ρ−+ mω 2 ρ2 .
Hrel = −
3N
−4
2
2m ρ
∂ρ
∂ρ
ρ
2

(3.14)

Hence the kinetic energy can be parametrically expressed in terms of a radial component andan angular part, given in general by the grand angular momentum Λ2N −1 . Like the surface area
element, this angular operator can be defined recursively,Λ2k = Π2k +

Λ2k−1
lk2
+
cos2 αksin2 αk(3.15)

withΠ2k = −

∂2
(3k − 6) − (3k − 2) cos(2αk ) ∂+.
2sin(2αk )∂αk∂αk(3.16)We are only considering the ` = 0 states, so that the angular momentum operators lk2 can beneglected. Moreover, only the leading term, which corresponds to two-body motion, of the recursion
relation in Eq. (3.15) is relevant5 ,Π2 ≡ Π2N −1 = −∂2(3N − 9) − (3N − 5) cos(2α) ∂.
+2∂αsin(2α)∂α(3.17)The Schrödinger equation Hrel Ψ = Erel Ψ, where Erel is the energy of the relative motion,using Eq. (3.14) represents a partial differential equation in 3(N − 1) coordinates. For convenience,
we exploit the adiabatic hyperspherical expansion method [108, 104], where we first solve the angular
5

Recall that we choose to track only the coordinates (ρ, αN −1 ), so that our wave function will be independent of

all αk except αN −1 .

23component of the Schrödinger equation for fixed ρ, and then tackle the hyperradial motion.6 To dothis, we introduce a set of basis functions defined on the hypersphere that diagonalizes the fixed-ρ
Hamiltonian at each hyperradius ρ. That is, the functions are eigenstates of hyperangular kinetic
energy:Λ2N −1 Y{λ} (ρ; Ω) = {λ} (ρ)Y{λ} (ρ; Ω).

(3.18)The set {λ} represents a set of quantum numbers, here unspecified, that serve to distinguish thevarious basis states. In the case of non-interacting particles, a = 0, these are the usual hyperspherical harmonics and are independent of ρ [106, 107]. In the present circumstance, however,
eigenfunctions of Λ2N −1 are crafted subject to the Bethe-Peierls boundary conditions, in which casethese functions depend also parametrically on ρ, as well as on a, as we will see in the next section.
We will refer to the functions Y{λ} as channel functions. Because they form a complete seton the hypersphere, it is possible to expand the relative wave function asΨ = ρ−(3N −4)/2

X

F{λ} (ρ)Y{λ} (ρ; Ω)(3.19){λ}for some set of hyperradial functions F{λ} . Using this expansion in Hrel Ψ = Erel Ψ and projecting
the resulting expression onto Y{λ0 } yield a set of coupled equations,"#!
X~2 ∂ 2~2 (3N − 4)(3N − 6) Λ2N −1
1
2 2
−
+
+
+ mω ρ − Erel F{λ} (ρ)Y{λ} (ρ; Ω) = 0.
2m ∂ρ2 2m
4ρ2
ρ2
2
{λ}

(3.20)
Here, the first term in the [ ] is a hyperradial kinetic energy and the second term is an effective
centrifugal energy that is a consequence of the hyperspherical coordinate system.
The set of coupled Eqs. (3.20) is still exact, if all terms in the expansion are kept, but this isprohibitively expensive [44]. Instead, we follow common practice and make a Born-Oppenheimerlike approximation [95, 92, 110, 111]. Namely, we assert that ρ is a “slow” coordinate in the
6

The adiabatic approximation was originally proposed by Born and Oppenheimer for calculating the structure ofa diatomic molecule. It was used in hyperspherical coordinates in studying nuclear systems by Macek[108] and Fabrede la Ripelle[109].

24
sense that we ignore the partial derivatives ∂Y{λ} /∂ρ in Eq. (3.20). When this is done, eachterm in the sum is independent of the others. The solution representing the BEC is a single
wave function specified by the quantum numbers {λ}. We adopt this approximation in what
follows; the derivative couplings between adiabatic functions can be reinstated by familiar means[95, 92, 110, 111]. It is also worth noting that in the limit of infinite scattering length, the BornOppenheimer approximation again becomes exact with the zero-range model [88].
Here we will take this procedure one step further. The particular adiabatic function of interest
to us will describe two-body correlations, and will be chosen to reduce the collective set of quantum
numbers {λ} to a single quantum number ν:Ψ = ρ−(3N −4)/2 Fν (ρ)Yν (ρ; Ω).

(3.21)The selection of this function will be described in the following section.
In the sense of having to use only a single quantum number, this approximation is not unlike
the K-harmonic approximation, in which Y is taken to be independent of all its arguments7 [45].
Such an approximation affords an easy calculation of ground state energies for small scattering
lengths. Here we extend the definition of Yν to include a dependence on both hyperradius and on a
restricted subset of hyperangles Ω that emphasize two-body correlations. This additional flexibility
will enable us to describe these correlations at any value of the scattering length.
Using a single adiabatic function, the Schrödinger equation becomes a single ordinary differential equation in ρ:



~2 d2
~2
diag
−
+
V
(ρ)
F
(ρ)
+
hν|Λ2N −1 |νiFν (ρ) = Erel Fν (ρ),ν
2m dρ22mρ2(3.22)

where
V diag (ρ) =7~2 (3N − 4)(3N − 6) 1
+ mω 2 ρ22m
4ρ2
2(3.23)

Note that in the K-harmonic approximation, ν = 0 and an explicit two-body interaction is built into the the

Schrödinger equation as a linear function of a. This method might suffice for small a and just one basis ν = 0 but
would not for larger a.

25
is the potential whose ground state supports the non-interacting (a = 0) condensate wave function. The additional term involving this matrix element hν|Λ2N −1 |νi, which is an integral over thehypersphere, represents the additional kinetic energy in the many-body wave function due to theBethe-Peierls boundary conditions. It can be viewed as a kind of “interaction” potential, since thescattering length responsible for these boundary conditions arises ultimately from the two-body interaction. Section 3.2.4 outlines how this is evaluated and an explicit representation of this quantity
in terms of the quantum number ν is given in Eq. (3.53).
As a point of reference, the interaction term vanishes in the limit of zero scattering length. The
non-interacting gas therefore has approximately the energy of the potential V diag at its minimum,
i.e., the energy of the gas isp(3N − 4)(3N − 6)
V diag (ρ0 ) =
~ω,
2

1/4
(3N − 4)(3N − 6)
ρ0 =
aho ,
4
where aho =

3.2(3.24)(3.25)

p~/mω is the characteristic trap length.

Basis functions
In this section, we construct the adiabatic basis functions Yν , focusing on the most relevant

one to describe the gas-like state of the BEC. We also construct the approximate matrix elementshν|Λ2N −1 |νi required to solve the hyperradial equation in (3.22).
3.2.1

Jastrow form and pair wave function

In describing a dilute Bose gas with two-body interactions, for our present purposes, we arecontent with a simple description that includes only atom pair distances rij , which explores only atiny fraction of the available configuration space. Specifically, consider a single pair described by
r12 . Because we are looking into ` = 0, we can ignore the direction of the pair distance vector r12 .
√The relative motion of this pair is given by the single hyperangle α through sin α = r12 / 2ρ. We

26
can therefore contemplate a hyperangular basis functionφν (ρ; α) = φν (ρ; α12 )(3.26)

that is a function of only one of the 3N − 4 hyperangles, and that may depend parametrically onρ.
Indeed, one may do the same for any pair distance rij , and define a hyperangle via sin αij =√rij / 2ρ. Each such angle is expressed starting from a different set of Jacobi coordinates. Startingfrom this nugget of a wave function, one can build a basis function that is appropriately symmetrizedwith respect to exchange of identical bosons, via
QYν = qRi<j φν (ρ; αij )

dΩ

Q

2
i<j |φν (ρ; αij )|.

(3.27)

This form of the function as a product of all two-body contribution was made famous by Jastrow’s
pioneering effort [54]. It continues to find extensive use as a form of variational trial wave function,
especially for Monte Carlo studies of many-body physics [112]. Note that at this point we deviate
from the formalism of Ref. [46], which expresses symmetrization by means of a sum of two-bodycontributions (Faddeev approach) rather than a product (Jastrow approach).
The procedure for constructing and using adiabatic basis functions therefore consists of (i)choosing a reasonable set of functions φν ; and (ii) living with the consequences of this choice. To
begin, the function φν should satisfy the Schrödinger equation for the relative motion of two atoms
when they are close to one another, that is, for small α. We define φν to be an eigenfunction of thedifferential operator Π2 ,
Π2 φν (ρ; α) = ν (ρ)φν (ρ; α),(3.28)

where Π2 is given in Eq. (3.17). The solution is determined from α = 0 (when the two particles
coincide), up to a value of α = αd to be determined below. To help visualize the consequences ofthis equation, it is sometimes useful to make the substitutionφ̃ν (ρ; α) = sin α cos(3N −7)/2 αφν (ρ; α),(3.29)

27
which satisfies the differential equation∂29N − 19 (3N − 7)(3N − 9)2− 2−+tan α φ̃ν = ν φ̃ν .
∂α2
4(3.30)

This version takes the form of an ordinary Schrödinger equation in α, with a centrifugal potential
energy term ∝ tan2 α that confines the motion of the atom pair toward α = 0, and thus prevents
the atom pair from getting too far apart. This wave function therefore automatically emphasizes
the action of this pair over the interaction of these atoms with others. As the number of particles
grows larger, the confinement is restricted to smaller values of α.
The equation we will solve is, however, Eq. (3.28) with Π2 given by Eq. (3.17). To solve it,
we make the substitution z = sin2 α, leading to


 
3
3 3N
d
ν
d2 φν
+
−
z
φν + φν = 0.
z(1 − z) 2 +
dz
2
2
2
dz
4

(3.31)

This equation has the form of the hypergeometric differential equation [113, 114, 115], yielding
two independent solutions for φν , one regular and one irregular, in terms of the hypergeometricfunctions 2 F1 :
3
3N − 52fν (α) = 2 F1 −ν,+ ν, ; sin α22
1 3N − 61gν (α) = (sin α)−1 2 F1 −ν − ,+ ν, ; sin2 α .
222Here, ν is a to-be-determined index that will in turn determine the eigenvalue,









 fν 
 fν 
 fν 
Π2= ν
= 2ν(2ν + 3N − 5).
 gν  gν  gν 


(3.32)(3.33)(3.34)A perfectly general solution to Eq. (3.28) is thenφν = Afν + Bgν .

(3.35)

In general, both constants A and B, as well as the index ν, will depend on the hyperradius and thescattering length, as we will now show.

28
3.2.2

Boundary conditions

The coefficients A and B are determined by applying boundary conditions. The first suchcondition occurs at α = 0, where the two atoms meet, and is given by the Bethe-Peierls condition(3.3). We first write this condition in hyperspherical coordinates. We have, in the limit of small α,
and for fixed hyperradius ρ,√
1 ∂(r12 φν )∂( 2ρ sin αφν )11∂(αφν )
1
√
=√
≈√
=− .
r12 φν ∂r12
∂αa2ρ sin αφν ∂( 2ρ sin α)2ρ(αφν )(3.36)

Next, expanding the hypergeometric functions near α ≈ 0 gives φν ≈ A + B/α,8 whereby
√
−

2ρ
1
∂[α(A + B/α)]
A
=
= ,
a
α(A + B/α)
∂α
B
α=0

(3.38)

which determines the ratio of the coefficients. Note that this ratio depends on ρ. It is significant
that this boundary condition is applied to a single pair of particles and is then implicitly applied
to all pairs by the form of the wave function in Eq. (3.27). This is in accord with the notion thatthe Bethe-Peierls boundary condition is local, and influences each pair independently of what theother pairs are doing.
The other boundary condition on φν is inspired by the brilliant reinterpretation of the Jastrow
wave functions by Pandharipande and Bethe [55, 56]. In this version φν is viewed as a piece of thepair-correlation function, related to the probability of finding this pair a given distance apart. To
this end, φν is required to approach unity on an appropriate length scale rd , which in hypersphericalhpicoordinates we translate into an appropriate hyperangular scale αd = sin−1 rd /(2)ρ . Beyond
this characteristic scale, the atoms are assumed to be uncorrelated and the wave function is required8

Useful asymptotic expansion of hypergeometric function:
2 F1 (a, b, c; z)

≈ 1+

z→0

a(a + 1)b(b + 1) 2
abz+z + ....
c2c(c + 1)(3.37)

29
to satisfy the following:φν (ρ; α) = 1,

α > αd ,

∂φν
= 0.
∂α αd

(3.39)

Setting the function to unity for α > αd is convenient but arbitrary.
There remains the issue of determining a reasonable value of αd . This is also done by viewing|φν |2 as a pair-correlation function. Given the location of atom 1, α encodes the distance to thenext atom, 2. If α becomes too large, then atom 2 can go explore parts of the hypersphere wherea third atom is likely to be found. At that point, the description in terms of a pair-correlation
function is not so useful, and φν should not be extended non-trivially this far.
Therefore, α should be limited to a region of the hypersphere where, on average, only one
atom will be found in addition to the fixed atom 1. This is the essence of the lowest-order constrained variational, or LOCV, approximation [116, 55]. Given the symmetrized wave function Yνin Eq. (3.27), the average number of atoms within 0 ≤ α ≤ αd of the atom presumed to lie at α = 0
is
Rα
R
Q
4π 0 d dΩα dΩN −2 i<j |φν (ρ; αij )|2RQ(N − 1),dΩN −1 i<j |φν (ρ; αij )|2(3.40)

which returns N − 1 when αd = π/2. The expression given in Eq. (3.40) can be written asRα
4π 0 d dΩα hν (α)(N − 1),R π/2
4π 0 dΩα(3.41)Q Z
RdΩN −2 i<j |φν (ρ; αij )|2Qhν (α) = 4π dΩα R.
dΩN −1 i<j |φν (ρ; αij )|2(3.42)

where hν is defined throughThis is a difficult multidimensional integral to evaluate, but it is often conveniently expanded into
powers of integrals of the (presumed small) quantities 1 − |φν (ρ; αij )|2 . The lowest-order term of
this expansion, and the approximation we will use here, then gives the approximationhν (α) = |φν (ρ; α)|2 .

(3.43)

30
Using this approximation, and setting the average number of atoms (3.41) to unity, yields thenormalization criterion
Z αd2dα sin α cos

3N −7

0

Z π/2
1
α|φν (ρ; α)| =dα sin2 α cos3N −7 αN −1 0√1
π Γ(3N/2 − 3)=.
N − 1 4 Γ(3N/2 − 3/2)
2

(3.44)

This requirement must be met self-consistently. That is, the boundary conditions (3.38) and (3.39)
determine φν for any given αd ; but αd must also be chosen so that (3.44) is satisfied. Noticethat the right-hand-side of (3.44) scales as N −5/2 as N gets large9 , whereby αd gets smaller withincreasing N . Any pair of atoms must be closer together to avoid the other atoms, when there are
more atoms.
The LOCV approach has been employed successfully when Jastrow wave functions are used[78]. Usually, one posits an unknown two-body correlation function, to be determined variationally,
minimizing some energy while varying the parameters of the trial function. This method has also
been employed, in Cartesian coordinates, to describe the energetics of a homogeous Bose gas withlarge scattering length [78, 80, 117]. Used in the present context of hyperspherical coordinates, we
dub this approach the hyperspherical LOCV (H-LOCV) method.
Putting together the boundary conditions, we have√2ρA−
B = 0,afν0 (αd )A + gν0 (αd )B = 0,

(3.46)

where fν0 (αd ) ≡ (dfν /dα)|αd , and similarly for gν0 (αd ), can be determined from the derivatives of
the hypergeometric functions. This system of equations for A and B can be solved only if
1√− 2ρ/afν0 (αd )gν0 (αd )

= 0.

det

9

√
We apply Γ(z) ∼ z z−1/2 e−z π for z → ∞ to get


−3/2
3N −6 (3N −7)/2
Γ( 3N2−6 )
3N
2
≈
≈
,

3N −3 (3N −4)/22
Γ( 3N2−3 )2(3.47)

N → ∞.

(3.45)

31
Therefore, given a (which defines the physical problem to be solved) and ρ (which defines the hyperradius at which the adiabatic function Yν is desired), zeroing the determinant (3.47) determines
a value of ν for any given αd , which is then varied to self-consistently satisfy the normalization.

3.2.3

Renormalized scattering length

The procedure outlined above generates an entire spectrum of ν values denoting various pairexcitations of the condensate. The lowest member of this spectrum, with no nodes in α, represents
a self-bound liquid-like state [118, 117], and is not what we are interested in here. Rather, forpositive scattering length, the BEC wave function in the relative coordinate r12 should contain asingle node to describe the gas-like BEC ground state.
We denote the location of this node by ac (the subscript “c” denoting the value of r12 where
the wave function crosses zero). For small scattering length, this node lies at a distance ac ≈ a
as shown in Fig. 3.1(a). However, as a gets larger this node, confined to a hyperangular range
0 ≤ α ≤ π/2, must saturate, leading to a finite ac for any finite hyperradius. The saturation
value of ac , termed a∗c , represents the effective scattering length on resonance. This saturation is
illustrated in Fig. 3.1(b), showing φν (α) for a = ∞.
Figure 3.2 tracks the value of the length ac over the entire range of positive scattering lengths.
For small a, ac grows linearly. It then rolls over, on length scales comparable to the harmonic
oscillator length, to saturate to a value a∗c when a = ∞. This auto-renormalization of the scattering
length is inherent in the H-LOCV method and independent of any local density approximation of
the gas.

3.2.4

Hyperangular kinetic energy integral hν|Λ2N −1 |νi

Having defined Yν and φν in the previous sections, we are now ready to evaluate the matrix
element hν|Λ2N −1 |νi that appears in Eq. (3.22). In general this matrix element involves complicatedmultidimensional integrals. However, far simpler, approximate versions of these integrals are often
possible by means of the cluster expansion [112]. This is the same set of ideas developed originally

32

1
0.8
0.2

0.6

0.1
0.4
0
0.2
-0.1
0
-0.2
0.03

-0.2

0.04

0.05

-0.4
0

0.5

1
(a)

1
0.8
0.6
0.2
0.4
0.2
0
0
-0.2

-0.2
0.7

-0.4
0.6

0.8

1

0.8
1.2

1.4

(b)

Figure 3.1: The angular wave function φν as a function of r12 for N = 10 atp(a) small scattering
length and on (b) resonance a = ∞. All length scales are in units of aho = ~/mω. The insets
show the zero-crossings of the curves. The crossing occurs at ≈ a for small scattering length andat some finite value√in the unitary limit. Note that φν is actually a function of the hyperangle α.
The relation r12 = 2ρ sin α is used to convert the angle α to the pair distance r12 . Here ρ has its
value at the minimum of the corresponding hyperradial potential.

33

10 0

7

10 -1

6
5

10 -2

4
3
10 -3

10 -3
10 -3

10 -2

10 -1

10 0

10 -1

10 1

10 1

10 2

10 3

10 3

Figure 3.2: The zero-crossing ac of the angular wave function φν for N = 10 as a function of thescattering length a. The inset shows the value ρmin that is used to compute ac ; ρmin is taken to be
the value of ρ where the effective potential V eff , discussed in the next section, is at its minimum.

34
in statistical mechanics to derive virial coefficients in the not-quite-ideal gas equation of state[119, 120, 121, 69].
We need to evaluate matrix elements of the grand angular momentum operator Λ2N −1 . Because our wave functions consider only one pair of atoms at a time, it should be sufficient to consider
only the leading term Π2 , acting on the pair (ij) = (12), and get the rest from symmetry.
To do so, let us for a moment return to independent-particle notation. The kinetic energy T
is a sum of single-particle operators that acts on a pairwise-symmetrized wave function:
N
X

Tk

k=1

N
Y

φν (ij).

(3.48)i<jThus,
hT i = hYν |N
X

Tk |Yν i = N hYν |T1 |Yν i,(3.49)

k=1

where Yν is defined in Eq. (3.27). Choosing a single atom, say atom 1, the operator T1 acts
identically on N − 1 of the terms in the product10 . Because N atoms do the same, the actionof T on the wave function can be written, for purposes of taking matrix elements, as (this is thesubstance of what Jastrow derives in his original paper [54])*YhT i = N (N − 1)φν (ij) [T1 φν (12)]i<j,(ij)= N (N − 1)i<j,(ij)φν (ij)i<j,(ij)6=(12)

*
Y

+
Y


φν (ij)T1 φν (12)φν (12)+ Y

φν (ij) ,(3.50)

i<j,(ij)

where T1 acts on a single pair, here chosen to be the pair (12).
√Translated into the Jacobi coordinate ηN −1 = r12 / 2, the kinetic energy acting on φν (12) isgiven by
~2 2
~2 2
~2
T1 = −∇1 = −∇12 = −
2m
2m
2m
10



1 2
∇
2 ηN −1=~2 1 2Π ,2m 2ρ2(3.51)

There is actually some subtlety in how the operator T1 acts on a Jastrow-type wave function. It should be notedthat T1 is a second-order differential operator hence cross terms are expected to appear. Inside the integral the crossterms vanish. See Appendix A for details.

35
ignoring the ρ-derivatives that are set to zero in the adiabatic approximation. Moreover, φν is aneigenfunction of Π2 in the interval 0 ≤ α ≤ αd , as described in (3.34). Beyond αd , φν (α) is constantand Π2 φν (α) vanishes (where the operator Π2 is defined in (3.17)). In other words, 2ν(2ν + 3N − 5)φν (α), 0 ≤ α ≤ αd~2 1
2
T1 φν (12) → Π φν (α) =×.
2m 2ρ2  0,α > αd(3.52)

2~
2
Thus, hT i → 2mρ
2 hν|ΛN −1 |νi, and for a single channel calculation, the expectation value of the

grand angular momentum is, therefore, given by1hν|Λ2N −1 |νi = N (N − 1)2ν(2ν + 3N − 5)
2

=

R
Rα
Q
4π 0 d dΩα dΩN −2 i<j |φν (ij)|2RQdΩN −1 i<j |φν (ij)|2N
2ν(2ν + 3N − 5),2(3.53)

where the last equality is obtained by enforcing the LOCV normalization condition that sets
Eq. (3.40) to unity.
To summarize: the adiabatic eigenfunction φν (ρ; α) that we seek has the following properties:(i) it satisfies the differential equation (3.28) in 0 ≤ α ≤ αd ; (ii) it satisfies boundary conditions(3.38) and (3.39) at a suitable αd , chosen so that (iii) the normalization (3.44) is satisfied; and(iv) the wave function that results has a single node in α. The algorithm to find such a function
is not terribly complicated, inasmuch as the wave function can be written analytically in terms ofhypergeometric functions. This procedure yields a single wave function with a particular value ofthe index ν, for each value of scattering length and hyperradius. For a single channel potential, theeffective potential V eff (ρ) corresponding to this ν is given by (see Eq. (3.22))
V eff (ρ) = V diag (ρ) + V int (ρ),(3.54)

where V diag (ρ) is given in Eq. (3.23), and the “interaction” potential, using the result in Eq. (3.53),
is given by
V int =~2 N
2ν(2ν + 3N − 5).
2mρ2 2(3.55)The effect of this interaction potential can be seen in Fig. 3.3. When a = 0 (solid curve),the ρ−2 kinetic behavior dominates for small ρ while V eff takes on the proportional to ρ2 behavior

36of the trapping potential as ρ → ∞. The ground-state solution to the non-interacting case isjust the well known solution to the quantum harmonic oscillator problem with relative energy3(N − 1)~ω/2. For small positive a (dashed line), V eff rises above the non-interacting V eff (a = 0),
and the local minimum ρ0 (a), where the condensate is centered, increases, indicating an expansion
in the overall size of the condensate. This behavior is consistent with the repulsive nature of
the contact interaction characterized by a positive scattering length. For small negative scattering
length (dotted line), the opposite is true: the atoms pull in towards each other due to the attractive
contact interaction; hence the decrease in energy and the condensate size.
Such features of V eff have been illustrated previously using the K-harmonic method [45].
In the K-harmonic method, however, V int is proportional to the scattering length, whereby thismethod suffers the same limitation to small a as does mean-field theory. In the hyperspherical
LOCV method, by contrast, the effective scattering length saturates and the effective potential
remains finite even in the a → ∞ limit. This V eff is the dash-dotted curve in Fig. 3.3.

37

60
50
40
30
20
10
0
0

2

4

6

8

10

Figure 3.3: The effective potential V eff as a function of the hyperradius ρ for a = 0 (non-interactingcase), a > 0 (repulsive interaction), and a < 0 (attractive interaction). For any a > 0 and N , alocal minimum always exists. However, for a given a < 0, there exists a maximum N when thelocal minimum of V eff starts to disappear. The highest V eff curve corresponds to the a → +∞case.

Chapter 4

Ground-State Properties of the BEC in the Hyperspherical LOCVApproximation1

In this chapter, we report on ground-state properties of the BEC in the H-LOCV approximation. Recall that the method begins by separating the center of mass energy (3/2)~ω, and thensolves for the relative energy Erel using the Hamiltonian (3.14). In the results of this chapter, we
report the full condensate ground-state energy E0 = Erel + (3/2)~ω.

4.1

General features of the ground state
The energy per particle versus scattering length is shown in Fig. 4.1 for N = 10. On theleft is the energy for attractive interaction a < 0, on the right for repulsive interaction a > 0. The
E0 (a < 0) curve connects smoothly with E0 (a > 0) at a = 0 with E0 = 3N ~ω/2, then increases
smoothly until it saturates in the large-a limit. A similar behavior can be observed for any N .
As is well known, a trapped gas is mechanically stable only for negative scattering lengthof small magnitude. In a hyperspherical picture such as this one, we find that with small a theinteraction potential V int ∝ ρa3 as shown in Eq. (4.4), and a collapse instability occurs when theattractive 1/ρ3 interaction potential overcomes the repulsive 1/ρ2 centrifugal potential and theeffective potential V eff in Fig. 3.3 has no classical inner turning point. On the left of Fig. 4.1, thecollapse region occurs near a ≈ −0.1 aho where a metastable condensate would cease to exist.
Generally, this means that for a harmonically trapped Bose condensate, for any negative value
1

Results presented in this chapter were published in Ref. [89].

39

50
45
40
35
30
25
20
15
10
-10 -1-10 -2-10 -3-10 -4

10 -3

10 -1

10 1

10 3

Figure 4.1: Ground state energy E0 as a function of the scattering length a for N = 10. The left
and right panels consider negative and positive scattering lengths.

40of a only a certain critical number Nc of atoms can be contained before collapse occurs. This is
an intrinsically small-|a| phenomenon, and was quantified via variational approach using Gaussian
ansatz to minimize the energy in Eq. (2.24) [122, 10], and also in hyperspherical terms in Ref. [45].
Our H-LOCV results are in agreement with these calculations, finding that Nc ∼ 0.671 aho /|a|.

4.2

Positive scattering length
For the rest of this thesis we will focus on the positive scattering length case. To this end,

the ground-state energy of the condensate is shown versus scattering length, for gases of N = 4
and N = 10 particles, in Figs. 4.2 and 4.3. The H-LOCV result is shown as a solid line.
In Fig. 4.2, for N = 4 atoms, the energy per particle in the H-LOCV method saturates ata finite value in the resonant limit, just over 2.5 ~ω per particle. This case, with N = 4, can becompared directly to an accurate numerical solution to the full four-particle Schrödinger equation,and is therefore an important case to check.
The full numerical calculation incorporates a two-body potential of the form√v0 exp[−(rjk / 2r0 )2 ],(4.1)

with range r0 = 0.025 aho and depth v0 adjusted to achieve the desired two-body scattering length.
This model also incorporates a repulsive three-body Gaussian potential to eliminate deep-lying
bound states of the system. Within this model, an energy spectrum is calculated using a basis setof correlated Gaussian functions. Further details are provided in Ref. [83].
The energy spectrum includes a great many states that represent bound cluster states, theanalogs of Efimov states for the trapped system. They are characterized by, among other things, a
dependence on the three-body potential. By contrast, the nearly universal state corresponding to
the gaslike BEC ground state is identified by its near independence from the three-body potential
and its vanishingly small amplitude at small hyperradii. The energy of the corresponding stateis shown as a dashed line in Fig. 4.2. The comparison shows that the H-LOCV method gets thesaturation value of the energy per particle approximately right, at least for N = 4 particles.

41
Also shown in Fig. 4.2 (dotted line) is the result of an alternative hyperspherical method thatsymmetrizes the wave function using a sum of two-body terms (Faddeev method [49, 46, 47, 123])
rather than a product as we use here in the H-LOCV method. The Faddeev method is accurate
for three atoms in a trap [82, 59], and the comparison in Fig. 4.2 suggests that it is viable for four
particles as well.
A difference occurs, however, in Fig. 4.3 for N = 10 atoms. Here a full numerical solution to the Schrödinger equation is not available, but we can compare the H-LOCV and Faddeev
hyperspherical methods side by side. The solid line represents the results from the H-LOCV approximation. Also shown are the results from the mean-field Gross-Pitaevskii (GP, dashed line)[124] and hyperspherical-Faddeev (dotted line) models [49, 46, 47, 123]. If we zoom in to thea/aho  1 domain (inset), we find that all models agree well in the weakly interacting regime.
In addition to these, other independent computations such as the K-harmonic method [45] anddiffusion Monte-Carlo [73] have predicted the weakly interacting system accurately.

3

2.5

2

Faddeev
LOCV
full-SE

1.5
10 -3

10 -2

10 -1

10 0

10 1

10 2

10 3

Figure 4.2: Ground state energy per particle E0 /N as a function of the scattering length a for
N = 4.

In the large-scattering-length limit, Fig. 4.3 shows the saturation of the energy per particle

42

5
1.58

4.5

1.56

4

1.54

3.5

1.52

3

1.5
10 -3

GP
Faddeev
LOCV

10 -2

2.5
2
1.5
10 -3

10 -2

10 -1

10 0

10 1

10 2

10 3

Figure 4.3: Ground state energy per particle E0 /N as a function of the scattering length a for
N = 10. The inset shows a blow-up of the weakly interacting regime.

43
to an asymptotic value, which grows with N (as we will show below, the energy per particle
scales as N 1/3 as anticipated by the Thomas-Fermi approximation [1]). It is well-known that the
GP approximation diverges in this limit. In addition, Monte Carlo calculations with hard spheres
cannot approach this limit, although Monte Carlo calculations with renormalized scattering lengths
can circumvent this issue [80].
Of particular interest is the hyperspherical Faddeev method (dotted line). Using the formalism of Ref. [49], one finds that the energy per particle saturates to a far lower value than for the
H-LOCV approach. This situation gets worse as the number of particles N increases. Note thatthe H-LOCV ν that goes into the V int plays a different role from the pure hyperspherical-Faddeev
ν [49, 46]; also, the latter does not place any restriction on the α domain. Using the formalism ofRef. [49, 46], we find that the asymptotic hyperspherical-Faddeev ν at large N approaches ν → 2
in the large-a limit. In this model, therefore, the energy per particle is a decreasing function of N ,and does not adequately describe the system in this limit.
One distinct benefit of the H-LOCV approach is the ease with which it extends to the large-N
limit. In this formalism, N is just a parameter in the differential equations and their boundary
conditions. For example, Fig. 4.4(a) shows the energy per particle for N = 1000 particles. Obtaining
these results is not computationally harder than the same result for N = 10. Also shown, in
Fig. 4.4(b), is the two-body contact, C2 , a thermodynamic quantity given by [125]C2 = 8πma2 ∂Em ∂E= −8π 2.
~2 ∂a
~ ∂(1/a)(4.2)

This quantity determines how the energy changes as the scattering length changes and has a dimension of inverse length. Physically, it describes the short-range behavior of pairs of atoms (explicitly
included in the H-LOCV method), or alternatively, the tail of the momentum distribution. In somecases, the intensive contact density C2 , which has a dimension of (length)−4 , is used.
Figure 4.4 shows the two-body contact per particle C2 /N for N = 1000. It increases slowly at
small a and then rises dramatically within the intermediate regime until it hits a maximum before
saturating at unitarity. A similar behavior is also observed for different N , and in the renormalized44
Thomas-Fermi (TF) method [1] albeit in different units.

4.3

Limiting cases
In the H-LOCV model, once the value of ν is determined, properties of the condensate follow

by solving the hyperradial equation for Fν . Even simpler, to a good approximation the relative
energy is given by the minimum value of the effective potential V eff (ρ). In two limiting cases,|a|  aho and a  aho , ν can be approximated analytically in the large-N limit, hence so can thecondensate’s energy and contact. These details are discussed in Appendices B and C. Here, we
explore the analytical results that follow, and compare these results to others in the literature.

4.3.1

Small scattering length

When |a| is much smaller than the trap length aho , we find in Appendix C that1
ν ≈N 3 2r

3 a 3/2
N ,
πρ

if

|a|
 1.
aho

(4.3)

In this case the interaction potential Eq. (3.55) takes the perturbative form
V

int

~2 3
(ρ) =
mρ2 4

r

3 7/2 aN,πρ

for large N .

(4.4)

This potential has the familiar a/ρ3 scaling found already in the K-harmonic approximation and its
pextensions [45, 49]. The size of the perturbation can be estimated by taking ρ ≈ ρ0 ≈ 3N/2 aho
at the minimum of the non-interacting effective potential (see Eq. (3.25)). Since V diag (ρ0 ) is of theorder N , then V int can be considered perturbative if |a|/ρ  N −5/2 . Then the total energy of the
weakly interacting system can be approximated by
E0 ≈

3N
1
a
~ω + √ N 2
~ω.
2
aho
2π

(4.5)

A similar perturbative energy result emerges from the K-harmonic [45] and hyperspherical-Faddeev[49, 105] methods. See inset of Fig. 4.4(a).

45

25
20

GP
Faddeev
K-harmonic
LOCV
3.5
3

15

2.5

10

2

5

1.5
10 -3

10 -3

10 -2

10 -1

10 0

10 1

10 -2

10 2

10 3

1

10

(a)

90
80
70
60
50
40
30
20
10
0
-3
10

10

-2

10

-1

10

0

10

2

(b)

Figure 4.4: (a) Ground state energy per particle, in units of ~ω. The inset shows the various
perturbative versions (small a region). (b) Two-body contact per particle, in units of a−1
ho , as afunction of the scattering length a for N = 1000.

46
4.3.2

Infinite scattering length

Of arguably greater interest is the resonant limit of infinite two-body scattering length, which
is far from perturbative in many models. However, the H-LOCV method is capable of producing
analytical results in this limit as well. The details are worked out in Appendix D.
In this limit and at large N , the index ν has the form"  #2x0 27 1/6
ν∞ = √
N 2/3 ≈ 2.122 N 2/3 ,
6 2π

(4.6)

where x0 ≈ 2.798 is the root to a transcendental equation defined in Appendix D. This asymptotic
value of ν∞ is independent of hyperradius in the a  ρ limit, whereby it is easy to find the minimumof the effective potential
V

eff

~2
(ρ) ≈
2mρ2



9N 2
4



1
~2
3x20
+ mω 2 ρ2 +N 8/3 .
2
2mρ2 (16π)1/3

The minimum is located at hyperradius"
  #1/2
x0 27 1/6
ρ∞ = √
N 2/3 aho ≈ 1.588 N 2/3 aho .
2 2π

(4.7)

(4.8)

From this, the condensate ground state energy is presented in Appendix D up to order 1/a (andignoring the center of mass energy that is small in the large-N limit)p1  ρ∞ ~ωE0 ≈ 3N 2 ν∞ 1 − β
2
a
"
#



1/12
2+1a
16π
27 1/6xho≈x0 N 4/3 1 −N −1/6 0 7/2
~ω
N →∞
16π
27
a
x0

aho 
≈ 2.52N 4/3 1 − 0.254N −1/6
~ω,
a

(4.9)

where the function β is defined in Eq. (D.17). A consequence of this expansion is that we have ananalytic expression for the contact of the resonant gas,
C2 /N ≈ 16.1N 1/6 /aho .

(4.10)The pair wave function used to construct the ground-state wave function at unitarity has a√node at a∗c = 2ρ∞ αc where αc is defined in Eq. (D.3), or
1 π 16π 1/12 −1/6∗
ac ≈ √
Naho ≈ 0.989 N −1/6 aho .
(4.11)
x0 2 27

47
This quantity serves as the effective scattering length and is about 30% greater than the renormalized scattering length, a∗T F , found in Ref. [1]:a∗T F =

2.182,
2(6π )1/3 hn1/3 i

hn1/3 i ≈ 0.4282

N 1/6 2.
aho

(4.12)

Further, in this limit the breathing mode frequency of the condensate is easily derived. It issimply the oscillation frequency in the hyperradial potential, given bys
1 d2 V eff.
ωb =m dρ2 ρ=ρmin
On resonance, this quantity is given bys 
ωb = ω

6

3N 2 ν∞
2



ρ∞
aho

(4.13)

−4+ 1 = 2ω,(4.14)a result already worked out long ago based on symmetry considerations [88, 126].
The ground state of a resonant Bose gas has been considered previously for a homogeneousgas, as well as, more recently, for a trapped gas. Listed in Table 4.1 are the resulting energies fromthese previous calculations.3 Because the H-LOCV method is intrinsically tied to a trapped gas, it
is hard to compare directly to other calculations. However, the result of Ref. [1] provides a link,by first calculating fixed-density quantities, then translating them into trapped values by means of
the local-density approximation.
In the homogeneous case, the energy per particle in the resonant limit a → ∞ is expected to
be a multiple of the characteristic (Fermi) energy ~2 n2/3 /2m associated with the density n. Several
such values are reported in the table, spanning a range of about a factor of 4 for the uniform
system, and 2 for the trapped gas. In the case of Ref. [1], the renormalized scattering length
affords a different energy at each value of density, whereby the energy can be represented in termsof the mean value ~2 hn2/3 i/2m averaged over an assumed Thomas-Fermi density profile. Using
this same density profile, one can write the energy per particle in terms of the characteristic energy
scale of the trap, ~ω, whereby a direct comparison can be made with the H-LOCV method 4 . This2

See Appendix B for hn1/3 i computation.

3∂E
Note that some authors reported the chemical potential µ. We use the relation µ = ∂Nfor conversion.

4

See Appendix B for hn2/3 i computation.

48
is done in the third column of Table I. Just as the homogeneous LOCV appears to come in on thehigh side for the ground state energy, so too does the H-LOCV method for the trapped gas.
Table 4.1: Ground state energy per particle, computed by various methods. The uniform gas
energies are given in units of ~2 n2/3 /2m, while the result of Ref. [1] for a trapped system is given inunits of ~2 hn2/3 i/2m with the mean density hn2/3 i determined by averaging over a Thomas-Fermi
profile in a trap. This same averaging allows this energy to be written in terms of the trap frequency
ω in the final column.

E0 /N (uniform gas)(~2 n2/3 /2m)

Cowell (LOCV) [78]

26.66

-

Song (Condensate amplitude variation) [43]

7.29

-

Lee (Renormalization group) [37]6.0-Diederix (Hypernetted chain) [79]

7.60

-

Borzov (Resummation scheme) [36]

8.48

-

Zhou (RG) [35]

8.11

-

Yin (self-consistent Bogoliubov) [42]

7.16

-

van Heugten (renormalization group) [38]

12.94

-

Rossi (Monte-Carlo LOCV) [80]

10.63

-

E0 /N (trapped gas)(~2 hn2/3 i/2m) (N 1/3 ~ω)Ding (Renormalized K-harmonic) [1]

12.67

1.205

H-LOCV [89]

-

2.52

Likewise, there are various estimates for the contact on resonance, summarized in Table
4.2. For a homogeneous system on resonance, one reports the intensive, density-dependent contact
density C2 = γn4/3 , where n is the gas density and γ is a dimensionless constant. This quantity,
like the energy, is subject to an array of values tied to the different methods.
For the trapped gas, fewer examples of the contact have been calculated. We again turn to49
Table 4.2: Contact densities, computed by various methods. For the uniform gas, the intensive
contact densities C2 are given in units of n4/3 , while the extensive contact C2 /N , in units of hn1/3 i,
is given in Ref. [1] for a trapped gas. Averaging over a Thomas-Fermi profile in a trap allows thisC2 /N to be written in terms of the characteristic trap length aho in the final column.

C2 (uniform gas)(n4/3 )

Rossi (Monte-Carlo LOCV) [80]

9.02

-

Diederix (Hypernetted chain) [79]

10.3

-

Sykes [34]12-van Heugten (renormalization group) [38]32-

Yin (self-consistent Bogoliubov) [42]158-

Smith (Virial theorem) [125]20-

C2 /N (trapped gas)

(hn1/3 i)(N 1/6 /aho )Ding (Renormalized K-harmonic) [1]

11.8

5.05

H-LOCV[89]

-

16.1

the work of Ref. [1] to make the link. For a trapped gas, one can describe an extensive contact C2 ,related to the contact density by C2 = γN hn1/3 i, using the local-density approximation (LDA). In
this approximation, the trap results of Ref. [1] give γ = 11.8 (see Table 4.2.). Further, using the
average value of n1/3 in the Thomas-Fermi approximation, Eq. (4.12), this result can be translated
into natural harmonic-oscillator units5 , yielding C2 = 5.05N 1/6 /aho . Compared to this, our value,
also in the table, is 16.1N 1/6 /aho . Like the energy, the H-LOCV seems to overestimate the contact
on resonance. It does, however, correctly identify the N 1/6 number dependence of the contact fora trapped gas.

5

See Appendix B for hn1/3 i computation.

Chapter 5
Two-Step Production of Resonant Bose-Einstein Condensates1

5.1

Motivation
A current standard procedure for producing resonant BEC prepares the gas at an initiallyweakly interacting state, then quenches the magnetic field near a Feshbach resonance. This can be
regarded as an instantaneous projection of the initial small-a many-body state Ψ0 onto a collectionof many-body states Ψf at a → ∞. With the H-LOCV method, one can write down approximate
representations of these initial and final target states. By working out the probability of the
projection |hΨ0 |Ψf i|2 , one comes to a conclusion that the sudden projection may, indeed, not be
the optimal way to generate a resonant BEC. We demonstrate the details in this chapter andpropose an alternative scheme for preparing a resonant BEC. The theory formulated here is built
on the results of the previous chapters. It may be useful to review some of the relevant and familiar
details before proceeding to the main point.
In Chapters 3 and 4, we illustrated how one can obtain an approximate ground-state BEC
solution using the hyperspherical method for any given scattering length, that is from a = 0 to
a → ±∞. A key aspect of the process is the exploitation of a coordinate representation of BEC
that is expressed in terms of an effective potential energy analogous to Born-Oppenheimer (B.-O.)
curves in molecular physics. Under the B.-O. approximation, the hyperradius ρ is treated as theslow coordinate. That is, at each value of ρ, the Schrödinger equation, is solved in the hyperangular2~
coordinates Ω to yield a set of ρ-dependent eigenenergies of the operator 2mρΛ2N −1 , along with the1

Results presented in this chapter were published in Ref. [127].

51
corresponding eigenfunctions Y{λ} (ρ; Ω). Taken as a function of ρ, these energies constitute a set ofadiabatic potential energy curves V{λ} (ρ) for the motion in ρ. Here {λ} denotes a set of quantum
numbers required to specify the wave function in all coordiantes of Ω. A coupled set of differential
equations is obtained if we expand the wave function Ψ in adiabatic hyperangular basis,Ψ = ρ−(3N −4)/2

X

F{λ} (ρ)Y{λ} (ρ; Ω),(3.19 revisited){λ}for some set of radial expansion functions F{λ} . The eigenfunctions Y{λ} (ρ; Ω) are eigenfunctionsof the 3N − 4-dimensional partial differential operator Λ2N −1 . In practice, we find approximate
eigenfunctions by invoking the Jastrow approximation, which factors the total wave function into
a product of pair wave functions of the form
QYν = qRi<j φν (ρ; αij )dΩQi<j φν (ρ; αij ).

(3.27 revisited)2

Within this approximation, each two-body function φν (α), which is given by Eq. (3.35) (alongwith Eqs. (3.32) and (3.33)), satisfies an ordinary partial differential equation (3.31), subject to theBethe-Peierls, pair correlation, and normalization boundary conditions in Eqs. (3.36), (3.39) and(3.44). By using this ansatz, each φν solves an ordinary differential equation and requires only asingle quantum number ν, which replaces the set of 3N − 4 quantum numbers {λ}.
Within the B.-O. approximation the hyperradial motion is independent within each channel
ν. This has been justified by experience in the case of the relatively small scattering lengths thatare found in nonresonant BEC experiments [45, 49]. Moreover, the B.-O. approximation has beenshown to be exact on resonance, at least when considering only two-body interactions, as we dohere [126, 88]. We will employ this approximation throughout, writing the wave function asΨ = ρ−(3N −4)/2 Fν (ρ)Yν (ρ; α)(3.21 revisited)

for the channel with the lowest value of ν, representing the ground state of relative excitation. The
ground state of Vν (ρ) represents the BEC ground state, and excitations in ρ represent breathingmodes of the condensate.

52
Using a single adiabatic function, the Schrödinger equation becomes a single ordinary differential equation in ρ:



~2 d2
diag
−
+V
(ρ) + Vν (ρ) Fν (ρ) = Erel Fν (ρ),
2m dρ2(3.22 revisited)where~2 (3N − 4)(3N − 6) 1
+ mω 2 ρ2 ,
2m
4ρ2
2
2
2
~ N
~
hν|Λ2N −1 |νi =2ν(2ν + 3N − 5).
Vν (ρ) =22mρ
2mρ2 2

V diag (ρ) =(3.23 revisited)(3.55 revisited)

V diag is the diagonal potential whose ground state supports the non-interacting condensate wave
function. Vν represents the contribution due to interactions, namely, the eigenenergy of Λ2N −1 asdetermined above. It is a function of both the hyperradius and the scattering length. Calculation of
the matrix element hν|Λ2N −1 |νi, which involves integration over the entire hypersphere, is shown in
Chapter 3 where a lowest-order constraint variational approximation has been applied to obtain ameaningful outcome even for a = ∞. In general, Vν must be determined numerically, but Chapter 4derives useful approximations to this potential in the very small and very large-a limits, which wewill employ in the following. Thus, for any scattering length a, we find a B.-O. potential
V a (ρ) = V diag (ρ) + Vν (ρ).

(5.1)

We will denote the associated hyperangular wave function by Φa (ρ; Ω) to make explicit the scattering length for which this function was calculated. Here the notation ν is suppressed, since it
is understood we are considering the lowest value of ν. Vibrational states in V a (ρ) constitutethe radial wave functions Fna (ρ), each vibration n describing a breathing mode excited above theground state condensate with n = 0. The states relevant to our model are, therefore, defined by
the scattering length a and the number of breathing quanta,|a, ni = ρ−(3N −4)/2 Fna (ρ)Φa (ρ; Ω).

(5.2)

53

1500

1000

500

0
0

10

20

30

40

Figure 5.1: The scale of the problem. Each curve represents an effective potential energy surface
for a BEC with a = 0 (bottom) and a = ∞ (top), in our hyperspherical representation. A BEC
having a = 0 (Gaussian centered at ρ = 12.2 aho ) has essentially no overlap with a resonant BEC
having a = ∞ (Gaussian centered at ρ = 33.8 aho ).

54
Figure 5.1 shows the B.-O. PES’s for the non-interacting (a = 0) and resonant (a = ∞)cases for a gas of N = 100 atoms. With a = 0, V 0 (ρ) = V diag (ρ). This potential, the lowest curveon the left, is exact. Its minimum value occurs at
ρ0where aho =

≈

N 3

q

3N
2 aho ,

(5.3)

p~/(mω). The topmost curve on the right is an approximate potential for the resonant

limit. In the large N limit, this potential is given by2
V

∞

~2
(ρ) ≈
N 3 2mρ2



9N 2+ 3c0 N 8/3
4

1+ mω 2 ρ2 ,
2

(5.4)

where c0 ≈ 2.122. For realistic values of N > 102 , the centrifugal term with 9N 2 /4 can be safely
neglected. Then the minimum of V ∞ is located near3
ρ∞ ≈ (3c0 )1/4 N 2/3 aho .
N 3

(5.5)

Near their minima, we approximate these potentials as harmonic oscillators:
3N
1~ω + m(2ω)2 (ρ − ρ0 )2 ,2
2
1
V ∞ (ρ) ≈ (3c0 )1/2 N 4/3 + m(2ω)2 (ρ − ρ∞ )2 .
2V 0 (ρ) ≈

(5.6)
(5.7)

In both cases, the excitation frequency of the radial breathing modes considered is exactly twice
the trap frequency, ωb = 2ω. For non-interacting bosons, the energies are well-known and are givenby [106]
3N − 3
En,K = ~ω 2n + K +,2

n = 0, 1, 2, ...,

(5.8)

and K = 0, 1, 2, ... is the quantum number associated with the hyperangular component. For
the resonant gas, the 2ω frequency was derived in Eq. (4.14), and was anticipated by symmetry
considerations in Refs. [126, 88]. Without considering three-body or higher order correlations, thesereferences also emphasize that the B.-O. approximation is exact in the a = ∞ limit. Corrections
beyond the B.-O. approximation arise because the adiabatic wave functions Φ change from one
2

See Eq. (4.7).

3

See Eq. (4.8).

55
value of ρ to the next. But this change is only effective if ρ changes significantly on the scale of a,i.e., the corrections are of order ρ/a and vanish in the infinite scattering length limit. Therefore,
if the atoms could be prepared in the state F ∞ Φ∞ that we describe, this state would be stableagainst non-adiabatic transitions to whatever other states there are that could lead to heating, loss,
and so on. This stability is reduced if we were to include explicit three-body correlations in the
wave function.
From the harmonic oscillator nature of the potential curves in Eqs. (5.6) and (5.7), theexpected ground-state hyperradial wave functions are Gaussians centered at the minima and with√root-mean-squared width of aho / 2:
0

F (ρ) =

2a2ho π∞1/4

F (ρ) =

2

1/4

a2ho πexp[−(ρ − ρ0 )2 /a2ho ],(5.9)

exp[−(ρ − ρ∞ )2 /a2ho ],(5.10)The unnormalized Gaussian functions F 0 and F ∞ for N = 100 are illustrated as Gaussian-shaped
humps at the bottom of the a = 0 and a = ∞ potential curves, respectively, in Fig. 5.1. From thispicture, we see that the centers are far away from each other such that quenching the gas suddenly
from a = 0 to a = ∞ will yield a low transfer probability. That is, the probability of the atoms
landing in the resonant BEC state F ∞ , upon a direct quench, is2Z2|h0, 0, |∞, 0i| =

∞

0

∞

dρF (ρ)F (ρ)
dρF (ρ)F (ρ)

Z

2
0

∞

dΩΦ (Ω)Φ (Ω)
2

Z
≤

0



≈ exp −1.3N 4/3 ,(5.11)

which is negligible for large N . For this reason, it appears that, while the sudden quench to a = ∞produces an interesting, non-equilibrium gas of strongly interacting bosons, it is unlikely to generate
the desired resonant ground state BEC. In this chapter, we present a two-step protocol to improve
the yield in resonant ground state BEC.

56

1400
1200
1000
800
600
400
200
10

15

20

25

30

35

Figure 5.2: The two-step scheme from non-interaction to small a then to resonance.

57

5.2

The two-step scheme

5.2.1

Franck-Condon factors

The miniscule overlap between the non-interacting and resonant wave functions F 0 and F ∞suggests that direct projection from a = 0 to a = ∞ will not yield a resonant ground state BEC.
As an alternative, we propose a two-step scheme, in which we identify an intermediate scattering
length a and its B.-O. curve V a (ρ). Such a potential curve is shown as the intermediate curvein Fig. 5.2. A good candidate for V a (ρ) is one that supports a set of vibrational excitations nsuch that the Franck-Condon factors |h0, 0|a, ni|2 and |ha, n|∞, 0i|2 describe a maximum transfer
probability to the resonant BEC. This idea is illustrated in Fig. 5.2 for N = 100 atoms.
Numerically calculated Franck-Condon (FC) factors are shown as color-map plots in Fig. 5.3
for N = 100; the x-axis is the scattering length, y-axis the vibrational state n, and the color
indicates the transition probability. In general, for the first step from the non-interacting state
to the intermediate state, the optimum transition occurs when a is small and for low n states,
decreasing quickly with increasing a and n as shown in Fig. 5.3(a). For the second step fromthe intermediate state to the final state, the transition is optimum when a and n are larger, anddiminishes slowly with decreasing a and increasing n as in Fig. 5.3(b). These two steps cannot be
individually at their maxima under the same conditions. However, the best overall yield, given by
the product of the FC factors, occurs when a is still small relative to the oscillator length and forhigher vibrational states. This is true for any large values of N . Further, the two-step transition
probabilities seem to decrease as a function of N . See the transition probability for N = 1000 in
Fig. 5.4.

58

0.025

150
0.02

0.015

100

0.01

50

0.005

0.06

0.08

0.1

0.12

(a)

0.05

150
0.04

0.03

100

0.02

50

0.01

0.06

0.08

0.1

0.12

(b)
10 -4
7

150

6
5
4

100

3
2

50
1

0.06

0.08

0.1

0.12

(c)

Figure 5.3: Franck-Condon factors from the (a) non-interacting to intermediate states |h0, 0|a, ni|2 ,(b) intermediate to resonant states |ha, n|∞, 0i|2 , and (c) the two-step transition probability|h0, 0|a, niha, n|∞, 0i|2 as functions of scattering lengths a and vibrational states n. Here, N = 100.

59

10 -5

4200
2

4000
1.5

3800
1

3600
0.5

3400

3200
0.032

0.034

0.036

Figure 5.4: The two-step transition probability distribution |h0, 0|a, niha, n|∞, 0i|2 as a functionof scattering lengths a and vibrational states n for N = 1000.

60
5.2.2

The optimum intermediate state

Since the intermediate state will have a small value of a/aho , we can use a perturbative
approximate expression for V a . In the limits of perturbative a  aho and large N , this is given by4
V a (ρ) ≈ V 0 (ρ) +
N 3

~2
a
d0 N 7/2 3 ,m
ρ

(5.12)

pwhere d0 = (3/4) 3/π ≈ 0.733. Because of the oscillatory nature of the intermediate radial wave
functions Fna (ρ), the FC factors are largest when the inner and outer turning points of V a (denotedby ρ1n and ρ2n , respectively,) coincide with the minima ρ0 and ρ∞ of the initial and final states,
a scenario suggested in Figure 5.2. We will make this idea more precise in what follows, but thisobservation enables us to approximately determine the optimum intermediate scattering length, a∗ ,that maximizes the product of FC factors.
Assume the intermediate state of potential V a has energy Ea . Then its inner classical turning
point is determined by V a (ρ1n ) = En . For small ρ, V a is well approximated by the interaction
term in Eq. (5.12) alone, whereby this criterion becomes~2
d0 N 7/2 ρa3 .
1nN 3 mV a (ρ1n ) = Ena ≈(5.13)

Likewise, near the outer classical turning point of V a , this potential is well approximated by the
harmonic oscillator potential, and so1
V a (ρ2n ) = Ena ≈ mω 2 ρ22n .
2(5.14)

Setting ρ1n = ρ0 and ρ2n = ρ∞ and using (5.3) and (5.5), we can solve for the optimal values ofscattering length and intermediate energy:11(ρ0 )3 (ρ∞ )2 N −7/2 4aho
1 ρ∞ 2
∗
E ≈
~ω.
2 aho
a∗ ≈

N 3 2d0(5.15)(5.16)

Using these approximations for N = 100, the results are a∗ = 0.145aho and E ∗ = 571.2~ω,and are comparable to the values a∗ = 0.0859aho and E ∗ = 598.9~ω, determined by numerically4

See Eq. (4.4).

61
maximizing the FC factors. Expressions (5.15) and (5.16) become better estimates for larger N .
For N = 1000, they yield a∗ = 0.0316aho and E ∗ = 1.26×104 ~ω, whereas the numerically optimized
values are a∗ = 0.0332aho and E ∗ = 1.28 × 104 ~ω.

5.2.3

Wave packet dynamics

While this static picture provides an overall motivation for the two-step procedure, it does notdescribe the dynamics involved. Roughly, upon the initial projection from a = 0 to the intermediate
value a∗ , a wave packet is formed at ρ1n . In approximately one half of the trap period, this wave
packet propagates to ρ2n , giving the condensate its maximum radial extent and preparing it for
projection onto the resonant BEC state. With a given intermediate potential V a , we describe thetime dynamics of the BEC by expressing the initial state after the first step as a wave packet
expanded in the basis of the vibrational states of V a :a|Ψ (t)i =∞Xa|a, niha, n|Ψ (t = 0)ie

n=0

−iEn t/~=∞X|a, niha, n|0, 0ie−iEn t/~ ,(5.17)

n=0

where at time t = 0, Ψa is at the ground state of the non-interacting potential with total energy
E ≈ 3N ~ω/2. The probability of projecting the wave packet onto the desired resonant BEC groundstate is given by
P (t) = |h∞, 0|Ψa (t)i|2 ,(5.18)

where
h∞, 0|Ψa (t)i =∞X

h∞, 0|a, niha, n|0, 0ie−iEn t/~ .

(5.19)

n=0

We first extract numerically the most appropriate choice for the intermediate scattering lengtha∗ by maximizing the product of the FC factors. We then compute the transition probability
at different times with the H-LOCV unitary BEC model found in Chapter 4 for N = 100 and
N = 1000. Figures 5.5(a) and 5.5(b) show that the first maximum transition occurring at aroundtm ≈ π/(2ω). It takes about half a period, T /2, for the BEC to expand to resonance startingfrom the left side of the V a ; the breathing mode frequency is close to 2ω, thus the dwell time is

62

0.5
0.4
0.3
0.2
0.1
0
0

2

3

2

3

(a)

0.5
0.4
0.3
0.2
0.1
0
0
(b)

Figure 5.5: Transfer probability for (a) N = 100 with a∗ = 0.0859aho , and (b) N = 1000 witha∗ = 0.0332aho .

63
tm ≈ T /2 = π/ωb = π/(2ω). Most importantly, the transfer probability is significant: it is 48% for
N = 100 and 37% for N = 1000. This is a far better yield in the resonant state that the direct
projection result in Eq. (5.11).
Figure 5.6 shows how the size of the BEC with N = 100 atoms, expressed in terms of themean hyperradius hρi, is changing over time. It starts with ρ = ρ0 , the size of the non-interacting
gas, and reaches ρ = ρ∞ , the size of the resonant BEC, at t ≈ tm . The peaks of P (t) and hρidecrease slowly over time as the wave packet gradually dephases. It is, therefore, worthwhile toinstigate the second projection, to resonance, at time t = T /2.
35

30

25

20

15

10
0

2

3

Figure 5.6: Mean radius of the BEC in the intermediate phase versus time before quench tounitarity for N = 100 and a = 0.0859aho .

5.3

Large N limit
In calculating the P (t) numerically, we notice that P (tm ) decreases with N . Determining

how P (tm ) scales with N is extremely useful. In this section, we outline a method to get a good
estimate for this scaling.

64
5.3.1

Overlap between hyperradial wave functions using the reflection formula

To determine the Franck-Condon factors |h0, 0|a, ni|2 and |ha, n|∞, 0i|2 , one has to evaluate
the overlap integrals
hFna |F 0 iρ =Z ∞dρ Fna F 0 ,(5.20)dρ F ∞ Fna .

(5.21)

0hF∞|Fna iρ =Z ∞0

Leading contribution to the Franck-Condon factors comes from the overlap of wave functions at theclassical turning points, where the wave functions Fna are sharply peaked. In between the turning
points, the wave functions are highly oscillating. Yet we can consider that the projections of Fna to
F 0 and F a are still localized to the turning points since the last wave functions are also localized(or close to zero where Fna is wildly oscillating). The idea that the Franck-Condon factors can be
estimated from properties of the potential near the turning points goes back to the early days ofquantum mechanics [128, 129]. It is widely used in theories of optical and Raman transitions inmolecules, and recently to photoassociation of cold atoms as well [130, 131, 132, 133, 134]. Outof these types of molecular spectroscopy studies, the reflection formula was developed [133, 135],
which we will adapt.
We first express Fna in terms of the energy-normalized wave function FE throughhFna |Fna0 i =Z ∞which leads to Fna =

0

 dEn
dEn
dρ Fna Fna0 = δ n − n0 =
δ (En − En0 ) =hFE |FE 0 i,dn
dn

(5.22)

pdEn /dnFE . Casting FE into phase-amplitude form, after Milne [136, 137],r
FE (k, ρ) ≈2m
ζ (k) sin [β(k, ρ)] ,π~2(5.23)

where the amplitude ζ and phase β satisfy



d2
1
2
+
k
(ρ,
E)
ζ − 3 = 0,
dρ2
ζ
dζ
1
− 2 = 0,dρ β(5.24)(5.25)

65
with the wave vectorrk(ρ) =

2m(E − V (ρ)).
~2(5.26)The rapid oscillations of FE in Eq. (5.23) will have negligible effect on the integrals in Eqs. (5.20)
and (5.21), where Fna is expressed in terms of FE , except when ρ is near a turning point whichis also a point of stationary phase. Away from a turning point, it is sufficient to use the WKB
approximations for the amplitude and phase:1ζ (k) = p,k(ρ, E)Z ρ πβ (k, ρ) =dρ0 k ρ0 , E + .
4
ρt(5.27)(5.28)

Near a turning point ρt , to avoid the divergence in the WKB amplitude, we expand the Milne phase
to second orderβ ≈ b0 + b1 (ρ − ρt ) +b2(ρ − ρt )2 + ...
2(5.29)

π
,
4
∂β
b1 =
= k(ρt , E) = 0,∂ρ ρ=ρtb0 =b2 =(5.30)(5.31)

∂2β
∂k
m
∂V
=
= − 2 ζ 2 [k(ρt )].
2∂ρ ρ=ρt
∂ρ ρt
~
∂ρ ρt

(5.32)

Now, with
r
Fna =

dEn
dn

r

2m
ζ (k) sin [β(k, ρ)] ,π~2(5.33)

the integrand Fna F 0 is sharply localized around ρ1n , the classical inner turning point. Thus,
r
r


Z ∞
Z ∞
dEn 2m
b2
a
0
0
a
0
2
hFn |F iρ ≈ F (ρ1n )dρ Fn (ρ) = F (ρ1n )ζ (k(ρ1n ))dρ sin b0 + (ρ − ρ1n ) .
dn
π~2
2
0
0
(5.34)
To evaluate the last integral, we use the formulaZ ∞2

dx cos x
0



Z ∞
=
0

(5.35)

1.
|∂V /∂ρ|ρ1n(5.36)



1
=
2

r

π
.
2

dx sin x

2

Finally, we arrive atrhFna |F 0 iρ ≈ F 0 (ρ1n )

dEn
dn

s

66
The other overlap factor in Eq. (5.21) can be approximated in a similar fashion; it is given by
s
r
dE
1
n
hF ∞ |Fna iρ ≈ (−1)n F ∞ (ρ2n ),(5.37)dn|∂V /∂ρ|ρ2nwhere the (−1)n accounts for the sign of the rightmost amplitude around the outer turning pointρ2n of the vibrational state if we set the leftmost amplitude around ρ1n to be always positive as
expressed in Eq. (5.36).

5.3.2

Overlap between LOCV hyperangular wave functions

To give a complete picture of the overlap between wave functions, the angular overlapshΦa |Φ0 iΩ and hΦ∞ |Φa iΩ should also be considered. Real calculation involves 3N − 4 dimensional
integrals since this is the size of the hyperangular space. However, here, we only consider the one
hyperangle, α, that describes the two-body interactions, and the large N case.
We start with a symmetrized Jastrow-type basis,Q
i<j φν (ρ; αij ),Yν = qRQdΩ i<j φν (ρ; αij )2(3.27 revisited)

where αij is parametrically related to the coordinate distance between to particles, rij through rij =√2ρ sin αij ; the function φν satisfies the Bethe-Peierls boundary condition. The other boundary
condition is set by treating |φν |2 as a pair correlation function such that if two atoms are more than√
distance rd = 2ρ sin αd apart, then they become uncorrelated or |φ(αij ≥ αd )|2 = 1. Therefore,
within a region bounded by αd , there is on the average only one other atom (out of N − 1) which
can be seen by a fixed atom, or
RαR
Q
4π 0 d dΩα dΩN −2 i<j |φν (ρ; αij )|2
1
R
Q
=
,
2
N −1
dΩN −1 i<j |φν (ρ; αij )|(3.40 revisited)

where dΩ = dΩN −1 = 4πdΩα dΩN −2 , and dΩα = sin2 α cos3N −7 α dα. If αd = π/2, then the rightside of (3.40 revisited) should be one. The full form of the pair correlation function g2 can be
written as
Z π/2
g2 (α) =

4π

dΩα
0

!R

dΩN −2RdΩN −1Q

Qi<j|φν (ρ; αij )|22
i<j |φν (ρ; αij )|,(3.42 revisited)

67
which is hard to evaluate. To lowest order, however, it is approximated to be g2 (α) = |φν (α)|2 . The
details of the whole procedure outlined above is found in Chapter 3. Given ρ and the scattering
length a, one can then find αd and φv . The angle αd becomes extremely small as N increases. Henceφ(ρ; αij ) is unity in large region of αij - this is an approximation that leads to g2 (α) = |φν (α)|2 .
In the following derivations, we will also treat all the pair wave functions φ(ρ, αi0 j 0 ) equivalent
to unity, except one pair namely, φ(ρ, α12 ) = φ(ρ, α). So,
a

Z π/2

0

hΦ |Φ iΩ ≈ N0 Na

dα α2 φa (ρ1n ; α)φ0 (ρ1n ; α)(5.38)

0
∞

Z π/2

a

hΦ |Φ iΩ ≈ N∞ Na

dα α2 φ∞ (ρ2n ; α)φa (ρ2n ; α)(5.39)

0where the N ’s are some normalization constants so that hΦ0 |Φ0 iΩ = 1, hΦa |Φa iΩ = 1, andhΦ∞ |Φ∞ iΩ = 1, andφ0 (ρ; α) = 1(5.40)a 1aφ (ρ; α) ≈ A 1 − √
,
2ρ α

√
cos 6N ν∞ α∞φ (ρ; α) = Bα

if α < αa(5.41)

if α < α∞ ,(5.42)

v∞ = c0 N 2/3(5.43)

The wave functions φa and φ∞ identically approach unity for α > αa and α > α∞ , which are givenby π 1/6

N −5/66 1/62π
α∞ =N −5/6 .
27
αa ≈(5.44)(5.45)

Note that αa and α∞ are extremely small for large N so that the integrals in Eqs. (5.38) and (5.39)
are integrated over large part of the α-space where φa and φ∞ are unity. The constants A and B

68
are determined from the continuity boundary condition at αa and α∞ :a 1A≈1+ √
,
2ρ αa
 
α∞
1 2π 1/6 −5/6
B=
=
N
,
c1
c1 27
p

c1 = cos
6N ν∞ α∞ ≈ −0.942.

(5.46)(5.47)(5.48)We then find
r

24
,
π3
#
r "
√
 
2 2a 24 a 2
24
1 + 3 αa + 3
αa + ... ,
Na ≈
π3
π ρ
π
ρ
r 

24
12
−5/2
N∞ ≈
1 − 3 γN+ ... ,π3π
 1/3
2π
c2
1
1
√
γ=
+ 2 − ≈ 0.199727
32c2c1 6c01
p
6N ν∞ α∞ ≈ 0.336.
c2 = sin
N0 =(5.49)(5.50)(5.51)(5.52)(5.53)

Finally, after a series of algebraic steps and careful bookkeeping of N -scaling of the relevant parameters, we find2
hφ |φ (ρ1n )iα ≈ 1 − 3
π
0

a

a
ρ1n

2
αa ,

(5.54)

hφ∞ |φa (ρ2n )iα ≈ 1 − 0.151N −5/2 ,(5.55)

which are our approximations for hΦa |Φ0 iΩ and hΦ∞ |Φa iΩ , respectively. For large N , these quantities are both essentially equal to one.
We now have useful representations of the components that go into ha, n|0, 0i and h∞, 0|a, ni.
Using Eqs. (5.36), (5.37), (5.54), and (5.55), the overlap integrals in Eq. (5.19) are approximatedto berha, n|0, 0i = hFna |F 0 iρ hΦa |Φ0 iΩ ≈ F 0 (ρ1n )

dEn
dn

s

1
|∂V a /∂ρ|ρ1nrh∞, 0|a, ni = hF ∞ |Fna iρ hΦ∞ |Φa iΩ ≈ (−1)n F ∞ (ρ2n )

dEn
dn

s

,

(5.56)
1|∂V a /∂ρ|,
ρ2n(5.57)

69
where dn/dEn is the density of vibrational states in the intermediate potential. For the hyperangular parts of the wave function we approximate hΦa |Φ0 iΩ ≈ 1 and hΦ∞ |Φa iΩ ≈ 1 since N is large
and a/ρ1n is small in Eqs. (5.54) and (5.55).

5.3.3

Transition amplitude and probability

Next we convert the discrete sum in Eq. (5.19) into a continuum integral over the energy andevaluate it at t = tm ≈ π/(2ω) around which the maximum transfer occurs.
In terms of the evaluated Franck-Condon factors in Eqs. (5.56) and (5.57), we write thetransition amplitude asah∞, 0|Ψ (tm )i ≈ e

iE0a tm∞XdEn(−1) F (ρ1n )F (ρ2n )
dn
n

0

s

∞

n=0

1
|∂V /∂ρ|ρ1ns1eiωn tm ,|∂V /∂ρ|ρ2n(5.58)

with E0a as the zero-point energy of the PES V a ,5 ωn ≈ (2 + ∆n )nω, where ∆n < 1 (∆n  1 forsmall a). Thus,(−1)n eiωn tm ≈ ei(nπ+ωn tm ) = ei2nπ = 1.

(5.59)

Also, using Eqs. (5.13) and (5.14),

1/3
∂V
m
≈ −3
En4/3 ,
∂ρ ρ1n~2 d0 N 7/2 ap∂V≈ 2mω 2 En .
∂ρ ρ2nConverting the discrete sum into an integral over energy,

P

n →

R

(5.60)(5.61)
dE, and using the form of F 0

and F ∞ in Eqs. (5.9) and (5.10), and noting that the resulting integrand is strongly peaked at√
E ∗ ≈ 3c0 N 4/3 ~ω/2 ≈ 1.26N 4/3 ~ω (see Eqs. (5.16) and (5.5)), we get
2(2d0 )1/6
a 1/6 −23/36 1|h∞, 0|Ψ (tm )i| ≈ √ √N×~ω3π( 3c0 )11/12 aho




Z ∞
(ρ1 (E) − ρ0 )2(ρ2 (E) − ρ∞ )2
dE exp −
exp −
,
a2ho
a2ho
0
a

5

a

(5.62)The phase factor eiE0 tm is actually irrelevant in the calculations, but we include it here for completeness.

702with ρ1 ≈ ( ~m d0 N 7/2 a)1/3 E −1/3 and ρ2 ≈p
2E/(mω 2 ) from Eqs. (5.13) and (5.14). Now, F ∞ (ρ2 (E))

is a peaky function of E. We can then use the saddle point approximation6 to solve the integral in
Eq. (5.62):
Z ∞
0





(ρ1 (E) − ρ0 )2(ρ2 (E) − ρ∞ )2
dE exp −
exp −
a2ho
a2ho
 


√ ρ∞

= ~ω πexp −aho

2d0 N 7/2 a4ho ρa2
∞

1/3

2 
− ρ0

a2ho



.


(5.63)

0.4
0.15
0.1

0.3

0.05
0
1.44

1.46

1.48

1.5
10 -3

0.2

0.1

0
10 -4

10 -3

10 -2

10 -1

Figure 5.7: Transfer probability of the BEC versus scattering length a for large N . Inset shows azoom-in profile of N = 105 .
Finally, expressing ρ0 and ρ∞ in terms of N , the resulting transition amplitude is

!2 

1/6

1/3

1/3 r
1/6
2(2d0 )
2daa3N 
√ N 1/36√ 0|h∞, 0|Ψa (tm )i| ≈exp −
N 13/18
−
.
5/24
aho
aho
2
3c0
(3c0 )
3(5.64)
where c0 and d0 are defined in Eqs. (5.4) and (5.12). Plots of P (tm ), calculated in this way, fordifferent N are shown in Fig. 5.7 as a function of the intermediate scattering length a. We see thatthe estimated maximum transfer for N = 103 is ∼ 33%, which is close to what the exact calculation
6

See Appendix E.

71gives, namely, 37%. The inset in Fig. 5.7 shows the sensitivity of the transition probability to theintermediate a for N = 105 . The intermediate a should at least be within 0.4% from the optimum
to get at least half of the maximum transfer. By maximizing Eq. (5.64) with respect to a, theoptimum scattering length a∗ is found to be 3/2 √33c0 −2/3a =N
aho ≈ 3.16N −2/3 aho ,2(2d0 )∗(5.65)

which matches the a∗ obtained in Eq. (5.15) using Eqs. (5.3) and (5.5). And the maximum transfer
isa2

max |h∞, 0|Ψ (tm )i|2 1/41
8−1/12≈ 1.028N −1/6 .
N
≈
1/8
3
(3c0 )(5.66)

To put this into context, for 85 Rb in a trap with frequency ω = 2π × 10 Hz, the oscillator
length is aho = 6.51 × 104 a0 . Starting with N = 105 non-interacting atoms in the trap, the twostep process would be optimized for a scattering length of a∗ ≈ 95.4 a0 with a theoretical transfer
probability of 15%.
The yield into the final state actually goes down as the number of atoms increases. Qualitatively, this is because the range of hyperradius from ρ0 to ρ∞ increases as N grows. From Eq. (5.3),
ρ0 grows as ∼ N 1/2 , whereas from Eq. (5.5) ρ∞ grows as ∼ N 2/3 . According to this scaling, thepoint of origin ρ0 and the point of termination ρ∞ for the wave packet grow farther apart as N
increases. From this, a physical picture emerges. The initial state, of a certain width, starts at
ρ0 and propagates to larger hyperradius, widening as it does so, as is the generic nature of wave
packets. By the time it has arrived at ρ∞ , it is wider than the width of the target state, which had
the same width as the original state. For larger N , the propagation distance is farther. Thus thewave packet broadens more during propagation, and its overlap with the target wave function is
reduced.
Finally, it is worth considering the effect of starting from a nonzero initial scattering length.
Some numerical experimentation finds that this would not produce a large effect. As an example,
consider N = 1000 85 Rb atoms in a spherically symmetric trap with frequency ω = 2π × 10 Hz. See

72
Fig. 5.8. Within the model where the initial scattering length a = 0, numerical optimization of thetransfer probability yields a 37% probability, passing through an intermediate state, via scattering
length 2200a0 . By contrast, starting with a more realistic scattering length for a stable, mean-field
BEC, for example a = 142a0 , raises the final transfer probability only to 39%, while changing
the intermediate scattering length to 2700a0 . This slight improvement is expected as the initialstate with nonzero scattering length is closer to the target state; that is the distance between theρa1 6=0 and ρ∞ is smaller than the distance between ρa1 =0 and ρ∞ . Thus less broadening of the wave
packet. Also the propagation time to the target state is a bit shorter. However, for greater numbersof atoms, the intermediate scattering length a∗ is reduced. Since the initial scattering length mustbe smaller than the intermediate scattering length, regarding the initial scattering length as small
becomes increasingly justified. Hence the yield in the resonant BEC is well approximated by thea = 0 initial state considered.
0.5

0.4

0.3

0.2

0.1

0
1.4

1.45

1.5

1.55

Figure 5.8: Transfer probabilities for N = 1000 starting from non-interacting and weakly interacting initial states.

Chapter 6

Conclusions and Prospects

Despite its simplicity, the H-LOCV method presented in this thesis provides a reasonable
qualitative description of the mechanically stable Bose gas on resonance. Notably, the method
affords analytical estimates of essential quantities such as energy per particle and contact, when
the scattering length is infinite. It must be remembered that the approximation used here is only
the lowest-order version of a hyperspherical theory, since it incorporates only two-body correlations.
Various improvements can be made, including:
1) Extension to excited states. We have so far incorporated only a single adiabatic channel
function, consistent with our immediate goal of approximating a ground state. Yet there exists a
whole spectrum of states corresponding to different φν (ρ; α). We can contemplate states in which
one or more particles are placed in excited states, corresponding to excitations of the BEC; or in
the nodeless state below the condensate state, standing for bound molecular pairs. We can also
contemplate placing all the pairs in this nodeless state, to approximate the liquid-like configurationof Refs. [117, 118]. In any case, having a spectrum of approximate energy eigenstates is a place tobegin looking at the dynamics of a BEC quenched to resonance, or to any value of a. Along with
this, non-adiabatic couplings and their effect can be evaluated.
2) Extension of the Hamiltonian. Thus far only two-body interactions have been contemplated, leading to pairwise Bethe-Peierls boundary conditions and universal behavior. In more
realistic treatments, additional three-body interactions are required and can also be incorporated.
In such a case, the hyperspherical basis set can be extended to incorporate triplets of atoms, just as

74pairs were used here. This involves adding an additional Jacobi coordinate and in principle several
new hyperangles. The machinery for this extension is well-known, yet incorporating it into the
H-LOCV formalism requires careful attention.
3) Extension of the Jastrow method. A key approximation in the H-LOCV method has been
to evaluate important integrals by approximating two-body correlation functions as in (3.43). Thisrather severe approximation can be ameliorated, for example by a perturbation expansion known
as the hypernetted chain approximation [112, 66].
In Chapter 5, we have presented a protocol designed to implant a nontrivial fraction of thetrapped atoms into a resonant BEC. It remains to be understood what the consequences of thispreparation step will be. It is not clear, for example, how the non-equilibrium gas produced in
the two-step method will begin to come to equilibrium, and whether this process is different fromthe case of a direct quench to resonance. It is equally unclear at present how three-body losses
would differ in the resonant BEC thus produced than in a gas of equivalent density. A useful initial
experiment might be to prepare the resonant BEC as proposed here, and compare its dynamics tothat of a gas of equal initial density as the resonant BEC, but jumped suddenly to resonance.
This experiment would unfortunately be clouded by another issue. Consider, for example,
that starting from a non-interacting BEC of N = 104 atoms, our protocol is expected to transfer
only one fifth of them to the resonant BEC. What becomes of the rest? They are presumably
projected onto other quantum mechanical states of the system, each of which has its own dynamics
and three-body loss rates. To address this, it is necessary to formulate a reliable theory of excitedstates, in our case in the hyperangular degrees of freedom.
The two-step scheme that we have proposed here allows one to vary an intermediate scattering
length. Another possible protocol to explore is to have a varying trap frequency. If the target state
has a → ∞ and is expected to be in a given trap frequency ωf , then one can prepare the initialstate with some small-a and in a smaller frequency ωi (ωi < ωf ). Within the H-LOCV model, if theinitial state has a = 0, then for some appropriate ωi < ωf , then the maximum transfer is found to75
be1 ≈ 1.54N −1/6 , which is slightly better than Eq. (5.66). This may improve further if the initial
a 6= 0. One can also explore the consequence of this second protocol outside the H-LOCV model.

1This is computed by Eli Halperin.

Bibliography[1] Y. Ding and C. H. Greene. Renormalized contact interaction in degenerate unitary Bose
gases. Phys. Rev. A, 95:053602, 2017.
[2] A. Einstein. Sitzber. Kgl. Preuss. Akad. Wiss., page 261, 1924.
[3] A. Einstein. Sitzber. Kgl. Preuss. Akad. Wiss., page 3, 1925.
[4] S. N. Bose. Plancks Gesetz und Lichtquantenhypothese. Z. Physik, 26:178–181, 1924.
[5] 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–201, 1995.
[6] K. B. Davis, M.-O. Mewes, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn,and W. Ketterle. Bose-Einstein condensation in a gas of sodium atoms. Phys. Rev. Lett.,
75:3969–3973, 1995.
[7] C. C. Bradley, C. A. Sackett, J. J. Tollett, and R. G. Hulet. Evidence of Bose-Einsteincondensation in an atomic gas with attractive interactions. Phys. Rev. Lett., 78:1687–1690,
1995.
[8] L. D. Landau and E. M. Lifshitz. Statistical Physics. Pergamon Press, 1980.
[9] R. K. Pathria and P. D. Beale. Statistical Mechanics. Elsevier, 3 edition, 2011.
[10] F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari. Theory of Bose-Einstein condensation in trapped gases. Rev. Mod. Phys., 71:463–512, 2017.
[11] N. Bogoliubov. On the theory of superfluidity. J. Phys. (Moscow), 11:23–32, 1947.
[12] E. P. Gross. Classical theory of boson wave fields. Ann. Phys., 4:57–74, 1958.
[13] E. P. Gross. Structure of a quantized vertex in boson systems. Nuovo Cimento, 20:454–477,
1961.
[14] E. P. Gross. Hydrodynamics of a superfluid condensate. J. Math. Phys., 4:195–207, 1963.
[15] L. P. Pitaevskii. Vortex lines in an imperfect Bose gas. Sov. Phys. JETP, 13:451–454, 1961.
[16] K. Huang and C. N. Yang. Quantum-mechanical many-body problem with the hard-sphereinteraction. Phys. Rev., 105:767–775, 1957.

77[17] T. D. Lee, K. Huang, and C. N. Yang. Eigenvalues and eigenfunctions of a Bose system of
hard spheres and its low-temperature properties. Phys. Rev., 106:1135–1145, 1957.
[18] M. J. Holland, D. S. Jin, M. L. Chiofalo, and J. Cooper. Emergence of interaction effects in
Bose-Einstein condensation. Phys. Rev. Lett., 78:3801–3805, 1997.
[19] M. R. Matthews, D. S. Hall, D. S. Jin, J. R. Ensher, C. E. Wieman, E. A. Cornell, F. Dalfovo, C. Minniti, and S. Stringari. Dynamical response of a Bose-Einstein condensate to a
discontinuous change in internal state. Phys. Rev. Lett., 81:243–247, 1998.
[20] L. V. Hau, B. D. Busch, Z. Dutton C. Liu, M. M. Burns, and J. A. Golovchenko. Nearresonant spatial images of confined Bose-Einstein condensates in a 4-Dee magnetic bottle.
Phys. Rev. A, 58:R54–R57, 1998.
[21] A. Altmeyer, S. Riedl, C. Kohstall, M. J. Wright, R. Geursen, M. Bartenstein, C. Chin,
J. Hecker Denschlag, and R. Grimm. Precision measurements of collective oscillations in the
BEC-BCS crossover. Phys. Rev. Lett., 98:040401, 2007.
[22] S. B. Papp, J. M. Pino, R. J. Wild, S. Ronen, C. E. Wieman, D. S. Jin, and E. A. Cornell.
Bragg spectroscopy of a strongly interacting 85 Rb Bose-Einstein condensate. Phys. Rev. Lett.,
101:135301, 2008.
[23] N. Navon, S. Piatecki, K. Günter, B. Rem, T. C. Nguyen, F. Chevy, W. Krauth, and C. Salomon. Dynamics and thermodynamics of the low-temperature strongly interacting Bose gas.
Phys. Rev. Lett., 107:135301, 2011.
[24] R. J. Wild, P. Makotyn, J. M. Pino, E. A. Cornell, and D. S. Jin. Measurements of Tan’scontact in an atomic Bose-Einstein condensate. Phys. Rev. Lett., 108:145305, 2012.
[25] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga. Feshbach resonances in ultracold gases.
Rev. Mod. Phys., 82:1225–1286, 2010.
[26] P. Makotyn, C. E. Klauss, D. L. Goldberger, E. A. Cornell, and D. S. Jin. Universal dynamics
of a degenerate unitary Bose gas. Nat. Phys., 10:116–119, 2014.
[27] B. S. Rem, A. T. Grier, I. Ferrier-Barbut, U. Eismann, T. Langen, N. Navon, L. Khaykovich,
F. Werner, D. S. Petrov, F. Chevy, and C. Salomon. Lifetime of the Bose gas with resonantinteractions. Phys. Rev. Lett., 110:163202, 2013.
[28] R. J. Fletcher, A. L. Gaunt, N. Navon, R. P. Smith, and Z. Hadzibabic. Stability of a unitary
Bose gas. Phys. Rev. Lett., 111:125303, 2013.
[29] C. Eigen, J. A., P. Glidden, R. Lopes, N. Navon, Z. Hadzibabic, and R. P. Smith. Universal
scaling laws in the dynamics of a homogeneous unitary Bose gas. Phys. Rev. Lett., 119:250404,
2017.
[30] R. J. Fletcher, R. Lopes, J. Man, N. Navon, R. P. Smith, M. W. Zwierlein, and Z. Hadzibabic.
Two- and three-body contacts in the unitary Bose gas. Science, 355:377–380, 2017.
[31] C. E. Klauss, X. Xie, C. Lopez-Abadia, J. P. D’ Incao, Z. Hadzibabic, D. S. Jin, and E. A.
Cornell. Observation of Efimov molecules created from a resonantly interacting Bose gas.
Phys. Rev. Lett., 119:143401, 2017.

78[32] C. Eigen, J. A. P. Glidden, R. Lopes, E. A. Cornell, R. P. Smith, and Z. Hadzibabic. Universal
prethermal dynamics of Bose gases quenched to unitarity. Nature, 563:221–224, 2018.
[33] R. P. Smith, S. Beattie, S. Moulder, R. L. D. Campbell, and Zoran Hadzibabic. Condensation
dynamics in a quantum-quenched Bose gas. Phys. Rev. Lett., 109:105301, 2012.
[34] A. G. Sykes, J. P. Corson, J. P. D’ Incao, A. P. Koller, C. H. Greene, A. M. Rey, K. R. A.
Hazzard, and J. L. Bohn. Quenching to unitarity: Quantum dynamics in a three-dimensional
Bose gas. Phys. Rev. A, 89:021601(R), 2014.
[35] F. Zhou and M. S. Mashayekhi. Bose gases near resonance: Renormalized interactions in acondensate. Ann. Phys., 328:83–102, 2013.
[36] D. Borzov, M. S. Mashayekhi, S. Zhang, J. L. Song, and F. Zhou. Three-dimensional Bosegas near a Feshbach resonance. Phys. Rev. A, 85:023620, 2012.
[37] Y. L. Lee and Y. W. Lee. Universality and stability for a dilute Bose gas with a Feshbachresonance. Phys. Rev. A, 81:063613, 2010.
[38] H. T. C. Stoof and J. J. R. M. van Heugten. Resummation of infrared divergencies in thetheory of atomic Bose gases. J. Low Temp. Phys., 174:159–183, 2014.
[39] J. P. Corson and J. L. Bohn. Ballistic quench-induced correlation waves in ultracold gases.
Phys. Rev. A, 94:023604, 2016.
[40] J. P. Corson and J. L. Bohn. Bound-state signatures in quenched Bose-Einstein condensates.
Phys. Rev. A, 91:013616, 2015.
[41] X. Yin and L. Radzihovsky. Postquench dynamics and prethermalization in a resonant Bosegas. Phys. Rev. A, 93:033653, 2016.
[42] X. Yin and L. Radzihovsky. Quench dynamics of a strongly interacting resonant Bose gas.
Phys. Rev. A, 88:063611, 2013.
[43] J. L. Song and F. Zhou. Ground state properties of cold bosonic atoms at large scattering
lengths. Phys. Rev. Lett., 103:025302, 2009.
[44] J. Avery. Hyperspherical Harmonics: Applications in Quantum Theory. Kluwer Academic
Publishers, Dordrecht, 1989.
[45] J. L. Bohn, B. D. Esry, and C. H. Greene. Effective potentials for dilute Bose-Einsteincondensates. Phys. Rev. A, 58:584–597, 1998.
[46] O. Sørensen, D. V. Fedorov, and A. S. Jensen. Two-body correlations in N-body boson
systems. Phys. Rev. A, 66:032507, 2002.
[47] O. Sørensen, D. V. Fedorov, and A. S. Jensen. Correlated N-boson systems for arbitrary
scattering length. Phys. Rev. A, 68:063618, 2003.
[48] T. K. Das and B. Chakrabarti. Potential harmonics expansion method for trapped interacting
bosons: Inclusion of two-body correlation. Phys. Rev. A, 70:063601, 2004.

79[49] T. Sogo, O. Sørensen, D. V. Fedorov, and A. S. Jensen. The zero-range approximation applied
to the n-boson problem. J. Phys. B: At. Mol. Opt. Phys., 38:10511075, 2005.
[50] T. K. Das, S. Canuto, A. Kundu, and B. Chakrabarti. Behavior of a Bose-Einstein condensate
containing a large number of atoms interacting through a finite-range interatomic interaction.
Phys. Rev. A, 75:042705, 2007.
[51] B. Chakrabarti and T. K. Das. Shape-independent approximation for Bose-Einstein condensates interacting through a van der Waals potential. Phys. Rev. A, 78:063608, 2008.
[52] M. L. Lekala, B. Chakrabarti, G. J. Rampho, T. K. Das, S. A. Sofianos, and R. M. Adam.
Behavior of trapped ultracold dilute Bose gases at large scattering length near a Feshbachresonance. Phys. Rev. A, 89:023624, 2014.
[53] T. K. Das. Hyperspherical Harmonics Expansion Techniques: Application to Problems in
Physics. Springer, 2015.
[54] R. Jastrow. Many-body problem with strong forces. Phys. Rev., 98:1479–1484, 1955.
[55] V. R. Pandharipande and H. A. Bethe. Variational method for dense systems. Phys. Rev. C,
98:1312–1328, 1973.
[56] V. R. Pandharipande and K. E. Schmidt. Variational calculations of simple Bose systems.
Phys. Rev. A, 15:2486–2495, 1977.
[57] J. P. D’ Incao, J. Wang, and V. E. Colussi. Efimov physics in quenched unitary Bose gases.
Phys. Rev. Lett., 121:023401, 2018.
[58] P. F. Bedaque, E. Braaten, and H.-W. Hammer. Three-body recombination in Bose gaseswith large scattering length. Phys. Rev. Lett., 85:908–911, 2000.
[59] E. Braaten and H.-W. Hammer. Universality in few-body systems with large scattering
length. Phys. Rep., 428:259–390, 2006.
[60] C. E. Klauss. Resonantly interacting degenerate Bose gas oddities. PhD thesis, University ofColorado at Boulder, 2017.
[61] J. Pino. Strongly Interacting Bose–Einstein Condensates: Probes and Techniques. PhD
thesis, University of Colorado at Boulder, 2010.
[62] K. Bergmann, H. Theuer, and B. W. Shorei. Coherent population transfer among quantum
states of atoms and molecules. Rev. Mod. Phys., 70:1003–1025, 1998.
[63] B. M. Garraway and K.-A. Suominen. Wave packet dynamics in molecules. Contemp. Phys.,
43:97–114, 2002.
[64] J. Dalibard. Collisional dynamics of ultra-cold atomic gases. In M. Inguscio, S. Stringari,and C.E. Wieman, editors, Bose - Einstein Condensation in Gases, pages 321–349. IOS Press,
Amsterdam, 1999.
[65] A. de Shalit and H. Feshbach. Nuclear Physics: Nuclear structure. John Wiley & Sons Inc,New York, 1974.

80[66] T. C. F. van Heijst. An HNC approach to strongly interacting ultracold Bose gases. Master’s
thesis, Utrecht University, 2011.
[67] L. Pitaevskii and S. Stringari. Bose-Einstein Condensation. Oxford Science Publications,
2003.
[68] E. Fermi. Motion of neutrons in hydrogenous substances. Ricerca Scientifica, 7:13, 1936.
[69] K. Huang. Statistical Mechanics. Wiley, New York, 1987.
[70] H. Bethe and R. Peierls. Quantum theory of the diplon. Proceedings of the Royal Society of
London Series A, 148:146–156, 1935.
[71] T. T. Wu. Ground state of a Bose system of hard spheres. Phys. Rev., 115:1390–1404, 1959.
[72] J. L. DuBois and H. R. Glyde. Bose-Einstein condensation in trapped bosons: A variational
Monte Carlo analysis. Phys. Rev. A, 63:023602, 2001.
[73] D. Blume and C. H. Greene. Quantum corrections to the ground-state energy of a trapped
Bose-Einstein condensate: A diffusion Monte Carlo calculation. Phys. Rev. A, 63:063601,
2001.
[74] S. Tan. Energetics of a strongly correlated Fermi gas. Ann. Phys., 323:2952–2970, 2008.
[75] S. Tan. Large momentum part of a strongly correlated Fermi gas. Ann. Phys., 323:2971–2986,
2008.
[76] S. Tan. Generalized virial theorem and pressure relation for a strongly correlated Fermi gas.
Ann. Phys., 323:2987–2990, 2008.
[77] S.-J. Jiang, W.-M. Liu, G. W. Semenoff, and Fei Zhou. Universal Bose gases near resonance:
A rigorous solution. Phys. Rev. A, 89:033614, 2014.
[78] S. Cowell, H. Heiselberg, I. E. Mazets, J. Morales, V. R. Pandharipande, and C. J. Pethick.
Cold Bose gases with large scattering lengths. Phys. Rev. Lett., 88:210403, 2002.
[79] J. M. Diederix, T. C. F. van Heijst, and H. T. C. Stoof. Ground state of a resonantlyinteracting Bose gas. Phys. Rev. A, 84:033618, 2011.
[80] M. Rossi, L. Salasnich, F. Ancilotto, and Flavio Toigo. Monte Carlo simulations of the unitary
Bose gas. Phys. Rev. A, 89:041602(R), 2014.
[81] T. Busch, B.-G. Englert, K. Rzazewski, and M. Wilkens. Two cold atoms in a harmonic trap.
Found. Phys., 28:549–559, 1998.
[82] S. Jonsell, H. Heiselberg, and C. J. Pethick. Universal behavior of the energy of trapped
few-boson systems with large scattering length. Phys. Rev. Lett., 89:250401, 2002.
[83] D. Blume, M. W. C. Sze, and J. L. Bohn. Harmonically trapped four-boson system. Phys.
Rev. A, 97:033621, 2018.
[84] F. Werner and Y. Castin. Unitary quantum three-body problem in a harmonic trap. Phys.
Rev. Lett., 97:150401, 2006.

81[85] M. Thøgersen, D. V. Fedorov, and A. S. Jensen. Trapped Bose gases with large positive
scattering length. EPL, 79:40002, 2007.
[86] D. V. Fedorov, A. S. Jensen, and M. Thøgersen. Bose-Einstein condensates and Efimov statesin trapped many-boson systems. Few Body Syst., 43:69–74, 2008.
[87] M. Thøgersen. Universality in Ultra-Cold Few- and Many-Boson Systems. PhD thesis, Århus
University, Denmark, 2009.
[88] F. Werner and Y. Castin. Unitary gas in an isotropic harmonic trap: Symmetry properties
and applications. Phys. Rev. A, 74:053604, 2006.
[89] M. W. C. Sze, A. G. Sykes, D. Blume, and J. L. Bohn. Hyperspherical lowest-orderconstrained-variational approximation to resonant Bose-Einstein condensates. Phys. Rev.
A, 97:033608, 2018.
[90] M. Fabre de la Ripelle. The potential harmonic expansion method. Ann. Phys., 147:281–320,
1983.
[91] V. Aquilanti, S. Cavalli, and G. Grossi. Hyperspherical coordinates for molecular dynamics
by the method of trees and the mapping of potential energy surfaces for triatomic systems.
J. Chem. Phys., 85:1362–1375, 1986.
[92] J. G. Frey and B. J. Howard. The calculation of the ground state energy of weakly boundvan der waals trimers using the method of hyperspherical harmonics i. the Born-Oppenheimer
and adiabatic approximations. Chemical Physics, 99:416–426, 1985.
[93] V. P. Brito, H. T. Coelho, and T. K. Das. N -body systems in the hyperspherical adiabatic
approach. Revista Brasileira de Fsica, 21:309–326, 1991.
[94] N. Barnea. Hyperspherical functions with arbitrary permutational symmetry: Reverse construction. Phys. Rev. A, 59:1135–1146, 1999.
[95] D. Blume, C. H. Greene, and B. D. Esry. Comparative study of He3 , Ne3 , and Ar3 andusing hyperspherical coordinates. J. Chem. Phys., 113:2145–2158, 2000.
[96] S. T. Rittenhouse, M. J. Cavagnero, J. von Stecher, and C. H. Greene. Hypersphericaldescription of the degenerate Fermi gas: s-wave interactions. Phys. Rev. A, 74:053624, 2006.
[97] M. Gattobigio, A. Kievsky, M. Viviani, and P. Barletta. Harmonic hyperspherical basis for
identical particles without permutational symmetry. Phys. Rev. A, 79:032513, 2009.
[98] N. K. Timofeyuk. Convergence of the hyperspherical-harmonics expansion with increasing
number of particles for bosonic systems. Phys. Rev. A, 86:032507, 2012.
[99] T. K. Das, H. T. Coelho, and M. Fabre de la Ripelle. Uncoupled adiabatic approximation
for the hyperspherical harmonic approach. Phys. Rev. C, 26:2281–2287, 1982.
[100] L. D. Faddeev. Scattering theory for a three-particle system. Sov. J. Exptl. Theoret. Phys.,
12:1014–1019, 1961.
[101] O. A. Yakubovsky. On the integral equations in the theory of N particle scattering. Sov. J.
Nucl. Phys., 5:937–945, 1967.

82[102] S. P. Merkur’ev and S. L. Yakovlev. Quantum N-body scattering theory in configuration
space. Theor. Math. Phys., 56:673–682, 1983.
[103] E. Braaten and H.-W. Hammer. Efimov physics in cold atoms. Ann. Phys., 322:120–163,
2007.
[104] E. Nielsen, D.V. Fedorov, A.S. Jensen, and E. Garrido. The three-body problem with shortrange interactions. Phys. Rep., 347:373–459, 2001.
[105] T. Sogo, O. Sørensen, A. S. Jensen, and D. V. Fedorov. Semi-analytic solution to the n-boson
problem with zero-range interactions. Europhys. Lett., 69:732738, 2005.
[106] Y. F. Smirnov and K. V. Shitikova. Method of K harmonics and the shell model. Sov. J.
Particles Nucl. (Engl. Transl.); (United States), 8:344–370, 1977.
[107] J. Avery, W. Bian, J. Loeser, and F. Antonsen. Fourier transform approach to potential
harmonics. Int. J. Quant. Chem., 63:5–14, 1997.
[108] J. H. Macek 1 831-43. Properties of autoionizing states of He. J. Phys. B: At. Mol. Phys.,
1:831–843, 1968.
[109] M. Fabre de la Ripelle. C. R. Acad. Sc. Paris, 274:104, 1972.
[110] J. S. Avery. Hyperspherical Harmonics and Generalized Sturmians. Kluwer Academic Publishers, New York, 2002.
[111] A. F. Starace and G. L. Webster. Atomic hydrogen in a uniform magnetic field: Low-lying
energy levels for fields below 109 G. Phys. Rev. A, 19:1629–1640, 1979.
[112] A. Polls and F. Mazzanti. Microscopic description of quantum liquids. In A. Fabrocini, S. Fantoni, and E. Krotscheck, editors, Introduction to Modern Methods of Quantum Many-Body
Theory and Their Applications, chapter 2, pages 49–119. World Scientific, 2002.
[113] I. S. Gradshteyn and I. M. Ryzhik. In A. Jeffrey, editor, Table of Integrals, Series, andProducts, page 1072. Academic Press, 5th edition, 1994.
[114] In S. Iyanaga and Y. Kawada, editors, Encyclopedic Dictionary of Mathematics, page 1480.
MIT Press, Cambridge, MA, 1980.
[115] D. Zwillinger. Handbook of Differential Equations. Academic Press, Boston, MA, 3rd edition,
1997.
[116] V. R. Pandharipande. Hyperonic matter. Nucl. Phys. A, 178:123–144, 1971.
[117] R. M. Kalas and D. Blume. Equation of state of cold atoms: A lowest-order-constrained variational study of systems with large non-s-wave scattering lengths. Phys. Rev. A, 76:013617,
2007.
[118] B. Gao. Zero-temperature phases of many-atom Bose systems. Phys. Rev. Lett., 95:240403,
2005.
[119] B. Kahn and G. E. Uhlenbeck. On the theory of condensation. Physica, 5:399–416, 1938.

83[120] W. Band. Condensation phenomena in a clustering Bose-Einstein gas. Phys. Rev., 79:871–
876, 1950.
[121] R. J. Riddell and G. E. Uhlenbeck. On the theory of the virial development of the equationof state of monoatomic gases. J. Chem. Phys., 21:2056–2064, 1953.
[122] G. Baym and C. J. Pethick. Ground-state properties of magnetically trapped Bose-condensed
rubidium gas. Phys. Rev. Lett., 76:6–9, 1996.
[123] O. Sørensen, D. V. Fedorov, and A. S. Jensen. Structure of boson systems beyond the mean
field. J. Phys. B: At. Mol. Opt. Phys., 37:93–116, 2004.
[124] C. J. Pethick and H. Smith. Bose-Einstein Condensation in Dilute Gases. Cambridge University Press, Cambridge, 2002.
[125] D. H. Smith, E. Braaten, D. Kang, and L. Platter. Two-body and three-body contacts for
identical bosons near unitarity. Phys. Rev. Lett., 112:110402, 2014.
[126] Y. Castin. Exact scaling transform for a unitary quantum gas in a time dependent harmonic
potential. C. R. Phys., 5:407–410, 2004.
[127] M. W. C. Sze and J. L. Bohn. Two-step production of resonant Bose-Einstein condensates.
Phys. Rev. A, 99:033606, 2019.
[128] E. U. Condon. Nuclear motions associated with electron transitions in diatomic molecules.
Phys. Rev., 32:858–872, 1928.
[129] J. G. Winans and E. C. G. Stueckelberg. The origin of the continuous spectrum of thehydrogen molecule. Proc. Natl. Acad. Sci. USA, 14:867–871, 1928.
[130] K.-A. Suominen. Theories for cold atomic collisions in light fields. J. Phys. B: At. Mol. Opt.
Phys., 29:5981–6007, 1996.
[131] J. Weiner, V. S. Bagnato, S. Zilio, and P. S. Julienne. Experiments and theory in cold andultracold collisions. Rev. Mod. Phys., 71:1–85, 1999.
[132] J. L. Bohn and P. S. Julienne. Semianalytic theory of laser-assisted resonant cold collisions.
Phys. Rev. A, 60:414–425, 1999.
[133] P. S. Julienne. Cold binary atomic collisions in a light field. J. Res. Natl. Inst. Stand. Technol.,
101:487–503, 1996.
[134] C. Boisseau, E. Audouard, J. Vigué, and P. S. Julienne. Reflection approximation in photoassociation spectroscopy. Phys. Rev. A, 62:052705, 2000.
[135] A. Jabloński. General theory of pressure broadening of spectral lines. Phys. Rev., 68:78–93,
1945.
[136] W. E. Milne. The numerical determination of characteristic numbers. Phys. Rev., 35:863–867,
1930.
[137] F. Robicheaux, U. Fano, M. Cavagnero, and D. A. Harmin. Generalized WKB and Milne
solutions to one-dimensional wave equations. Phys. Rev., 35:3619–3630, 1987.

84[138] C. M. Savage, N. P. Robins, and J. J. Hope. Bose–Einstein condensate collapse: A comparison
between theory and experiment. Phys. Rev. A, 67:014304, 2003.
[139] S. E. Pollack, D. Dries, M. Junker, Y. P. Chen, T. A. Corcovilos, and R. G. Hulet. Extreme
tunability of interactions in a 7 Li Bose–Einstein condensate. Phys. Rev. Lett., 102:090402,
2009.
[140] P. A. Altin, G. R. Dennis, G. D. McDonald, D. Döring, J. E. Debs, J. D. Close, C. M. Savage,and N. P. Robins. Collapse and three-body loss in a 85 Rb Bose–Einstein condensate. Phys.
Rev. A, 84:033632, 2011.
[141] H. W. Jackson and E. Feenberg. Perturbation method for low states of a many-particle boson
system. Ann. Phys., 15:266–295, 1961.

Appendix A

Jastrow-Jackson-Feenberg-LOCV approach to homogeneous Bose gas

Here we follow the lines of Ref. [112] and fill in some of the steps that lead to the derivationof ground-state energy of homogeneous BEC using the Jastrow-LOCV formalism in Ref. [78].
Let ΨJ = Πi<j f2 (rij ) be the Jastrow ansatz in Cartesian coordinates for the N -body wave
function. Invoking indistinguishability of particles and symmetry of the wave function, the expectation value of the kinetic energy isN
X~2 ΨJ ∇2i ΨJ~2 ΨJ ∇2i ΨJhT i = −= −N.
2m hΨJ |ΨJ i
2m hΨJ |ΨJ i

(A.1)

i=1

Using the Jackson-Feenberg identity[141]1 , one can show thatZ

1dR Ψ ∇21 Ψ =
∗

Z





dR Ψ∗ ∇21 Ψ + ∇21 Ψ∗ Ψ − 2 (∇1 Ψ∗ ) (∇1 Ψ) ,

4
Z



1
=
dR ∇21 Ψ Ψ − (∇1 Ψ) (∇1 Ψ)
2

where

R

dR =

∇1 ΨJ =

R

dr1 dr2 ...drN , and the wave function Ψ is assumed to be real. Also,NX∇1 f2 (r1i )
i=2

(A.2)

f2 (r1i )ΨJ ,∇21 ΨJ = ∇1 (∇1 ΨJ )
N
X
∇21 f2 (r1i ) f2 (r1i ) − (∇1 f2 (r1i ))2=+
2 (r )
f
1i
2
i=2

(A.3)


! N
X ∇1 f2 (r1j ) ΨJ .
f2 (r1i )f2 (r1j )NX∇1 f2 (r1i )
i=2

j=2

(A.4)
1

This can be easily shown by applying integration by parts where the boundary terms vanish.

86
Using these two expressions, we get∇21 ΨJΨJ − (∇1 ΨJ ) (∇1 ΨJ ) = Ψ2J

NX∇21 f2 (r1i ) f2 (r1i ) − (∇1 f2 (r1i ))2,2 (r )
f
1i
2
i=2

(A.5)

and in terms of the identity in Eq. (A.2), Eq. (A.1) can be written as~2 1
hT i = −N2m 2RdR Ψ2J= −N (N − 1)

PN (∇21 f2 (r1i ))f2 (r1i )−(∇1 f2 (r1i ))2
f22 (r1i )

i=2

(A.6)

hΨJ |ΨJ i
22R
2 (∇1 f2 (r12 ))f2 (r12 )−(∇1 f2 (r12 ))
2
dR
Ψ
J
~ 1
f 2 (r )212hΨJ |ΨJ i
R(∇21 f2 (r12 ))f2 (r12 )−(∇1 f2 (r12 ))2 R
dr3 ...drN Ψ2J~2 1 dr1 dr2
f22 (r12 )R= −N (N − 1)2m 2
dr1 ...drN Ψ2JZ∇2 f2 (r12 ) f2 (r12 ) − (∇1 f2 (r12 ))2~2 1
dr1 dr2 g2 (r12 ) 1
= −n2
2m 2
f22 (r12 )Z∇2 f2 (r) f2 (r) − (∇f2 (r))2
~2 1
= −nN
,
dr g2 (r)2m 2
f22 (r)
2m 2

(A.7)

(A.8)
(A.9)
(A.10)

where we used the definition of the pair correlation function g2 in Cartesian coordinates [112]:R
N (N − 1) dr3 ...drN |Ψ|2R
g2 (r12 ) =,n2
dr1 ...drN |Ψ|2

(A.11)

with n as the number density. In the last line, we expressed r12 = r for neatness.
In the LOCV method, g2 ≈ f22 . Thus,Z
2 f (r) f (r) − (∇f (r))2
∇
hT i
~2 1
2
2
2
= −n
dr f22 (r)2N
2m 2
f2 (r)
Z2
2∇ f2 (r)~dr f22 (r)
,
= −n
2m
f2 (r)(A.12)(A.13)

where we wrote the second term of the integral in Eq. (A.12) asZZ(∇f2 (r))22
dr f2 (r)= 4π
f22 (r)Z∂f2
dr r (∇f2 (r)) = 4π dr r
∂r



Z
∞
∂f2
∂
∂f2
= 4π r2 f2− dr f
r2
,
∂r 0
∂r
∂r
2

2

2



∂f2
∂r


(A.14)(A.15)and noted that1 ∂∇ f (r) = 2
r ∂r
2


r

2 ∂f2

∂r


.

(A.16)

87
The boundary terms in Eq. (A.15) vanish because of the LOCV boundary conditions f2 (r > rd ) = 1
and f20 (r ≥ rd ) = 0, where rd is some healing distance, which is chosen such that on average, there
is only one other particle within a sphere of radius rd around an arbitrary particle in the gas:
Z rd
4πn

dr r2 f22 (r) = 1.

(A.17)

0

The above mathematical steps are one way to show Jastrow’s statement [54], “The cross terms in
the kinetic energy vanish in the integration.”
Now, the pair particle wave function f2 satisfies the Schrödinger equation [78]



~2 d2
−
+ V (r) − λ2 rf = 0,m dr2

(A.18)

where λ2 is the pair energy. Note that this expression is similar to Eq. (2.1). Also, instead of having
an explicit two-body interaction potential V , the Bethe-Peierls boundary condition is used so that
−

~2 d2
~2 2
(rf
)
=
−
r∇ f = λ2 rf,m dr2m(A.19)

and
ELOCV
hT i
n
=
= 4π
N
N
2
which is the model used in Ref. [78].

Z rd0

dr r2 f 2

λ2 f
= 2πnλ2f

Z rd
0

dr r2 f 2 ,(A.20)

Appendix B
Mean Density Calculations at a → ∞Here we compute hn2/3 i and hn1/3 i using the Thomas-Fermi wave function with renormalized
scattering length at a → +∞.
For a → +∞, the Thomas-Fermi wave function is given by[1] #3/42 1/3 R2 − r 02
3
6π
TF
ψT F (r0 ) =16πN 2/3 ζ(+∞)a4ho√ !1/6 
256 2
ζ(+∞) 1/4 1/6
N ,RT F = aho
9
π
"

(B.1)

(B.2)

with ζ(+∞) = 2.182. The number density is defined as n(r0 ) = N |ψ|2 . So,
hn

2/3

Z
i=



0

dr ψT F N |ψT F |"= 4πN 2/3
=

2

2/3

Z
ψT F =

dr0 N 2/3 |ψT F |10/3(B.3)#5/2 Z
1/3
RT F
5/2
3 6π 2
dr0 r02 RT2 F − r02
42/3
16πN ζ(+∞)aho
0

5 32/31
N 1/3
p.
8 (2π)5/6 ζ(+∞) a2ho

(B.4)
(B.5)

Similarly, for hn1/3 i
hn

1/3

Z
i=



0

dr ψT F N |ψT F |

2

1/3

ψT F = 4πN

1/3

Z RT F0

8/3

dr0 r02 ψT F

#2 Z

2 1/3
RT F
2
3
6π
dr0 r02 RT2 F − r02= 4πN 1/3
4
2/3
16πN ζ(+∞)aho0 5/12 1/6
2
1
N 1/6
128 2
=
.
35π π
91/3 ζ(+∞)1/4 aho

(B.6)

"

(B.7)
(B.8)

Appendix C

Weak Interaction in the Large-N Limit

First, let us consider the noninteracting case. From Eqs. (3.55) and (3.46), the ground statesolution for the non-interacting system, a = 0, gives ν = 0, B = 0, andφν (αd ) = A = 1,

(C.1)

where we used Eq. (3.35) and the definition fν=0 (α) = 1 from Eq. (3.32). The LOCV boundary
condition (3.44) becomes
4
1≈ √
π

r
 3/2Z αd
6 5/2 αd3
3
5/2
2 2
N
dα α A = 3
N
,
2
π
3
0

(C.2)

where we used sin α ≈ α1 andΓ( 3N2−3 )Γ( 3N2−6 )

≈




3N −3 (3N −4)/2
3N 3/22,≈

3N −6 (3N −7)/2
2
2

N → ∞..

(3.45)

Thus,
αd ≈ (π/6)1/6 N −5/6 ,1

if a = 0,(C.3)

Note that we can assume that α  1 where 0 ≤ α ≤ αd since αd is expected to be a small quantity in the large

N -limit.

90
which is an extremely small quantity when N → ∞. We assume the same scaling for αd when
0 < a  aho . Now, for small nonzero a, B 6= 0, boundary condition (3.39) becomes2∂φν≈ Afν0 (αd ) + Bgν0 (αd ),
∂α αd


4
3N − 50fν (α) √ ≈− ν+ ν α,
32
3N να1,
0=

(C.5)

α1

gν0 (α) √13N − 61− 2 + 4 −ν −
≈
+ν ,
α
2
2
3N να1,

(C.6)

α1

which givesr
3 1 3/2A≈ −N .
B N →∞
2π ν(C.7)

With the Bethe-Peierls boundary condition (3.38), we get
1
ν≈
2

r

3 a 3/2
N ,
πρ

if

|a|
 1.
aho

(C.8)

Note that ν is not an absolutely small quantity as it also depends on N .

2Here, we used the following asymptotic series representation of 2 F1 :2 F1 (a, b, c; z)

≈ 1+

z→0

ab
z + ...
c(C.4)

Appendix D

Strong Interaction in the Large-N Limit

As a → +∞, we assume from Eq. (3.38) that A ≈ 0 and√

6N ν∞ α2
,
α

(D.1)

 X

∞
√
(−1)k
3N 1 2,;α=(3N να2 )k = cos 6N να2 .
F
−ν,
2 1
2 2
(2k − 1)!!k!

(D.2)

φν (α)≈

N >ν∞  12

Bα

−1


2 F1

3N 1 2
−ν∞ ,
, ;α
2 2


=B

cos

where we used

k=0

Note that α may be small but the product√6N ν∞ α need not be. Also, this wave function is zero

when α = αc :
αc =

π
1
√
.
2 6N ν∞(D.3)

Boundary conditions (3.39) and (3.44) yield the relationsαq dcos 6N ν∞ αd2qq2
0 = 1 + 6N ν∞ αd tan 6N ν∞ αd2

B=

(D.4)

(D.5)


 q

 3/2
2
5/2
3
sin
2
6N
ν
α
∞ d 
N αd
4
3
1
q
q
1≈ √
.
1 +
π 2
cos2 6N ν∞ αd2 2
2 6N ν∞ αd2(D.6)

From Eqs. (D.5) and (D.6), we get

αd =

2π
27

1/6

N −5/6 ≈ 0.7843 N −5/6 ,if a → +∞.

(D.7)

Equation (D.5) has the form
1 + x0 tan x0 = 0,

(D.8)

92
q
with solutions x0 ≈ 2.798, 6.121, .... If x0 = 6N ν∞ αd2 , then"  #2x0 27 1/6
ν∞ = √
N 2/3 ≈ 2.122 N 2/3
2π6

if x0 = 2.798.

(D.9)

Now, if 0 < aho /a  1, then

3N 1 23N 3 2−1, ; α + Bα 2 F1 −ν,, ;αφν (α) ≈ A 2 F1 −ν,
2 2
2 2
N >ν 12√√sin 6N να2
cos 6N να2
≈A √
+B
.
α
6N να2(D.10)

Boundary condition (3.39) yields
q
q
2 tan
1
+
6N
να
6N ναd2 √
A
d
q
= q
6N ν.
B
6N ναd2 − tan 6N ναd2

(D.11)

With the Bethe-Peierls boundary condition (3.38), we get the relation
1 + x tan x= −, wherex − tan xqx = 6N ναd2 ,
= √

ρ 1
.
3N ν a(D.12)(D.13)(D.14)

Let x = x0 − ∆. If  = 0, then we recover Eq. (D.5) and ∆ = 0. Suppose 0 <   1 so that ∆ isalso a small varying quantity. Then, using
tan x ≈

tan x0 − ∆,1 + ∆ tan x0

Eq. (D.12) gives∆≈

x20 + 1.
x20

(D.15)

Expressing x and  back in terms of ν and a,√

6N ναd = x0 −x20 + 1 ρ 1√.
x20
3N ν a(D.16)

Note that if the second term vanishes, then ν = ν∞ as defined in Eq. (D.9). Let ν = ν∞ (1 − β),where β is a function of 1/a. Then Eq. (D.16) leads to




p
1
x20 + 1ρ
1
1√
6N ν∞ 1 − β αd ≈ x0 −
1+ β
2
2
x20
3N ν∞ a⇒x2 + 1ρ
1
β≈2 0 3 √.
x0
3N ν∞ a(D.17)

93

Appendix E

Saddle Point Approximation

The saddle point approximation or the method of steepest descent is useful in evaluating
integrals of the form
Z x2
I=

dxf (x)eAg(x) ,

(E.1)

x1

where f and g are some real functions of x and A > 0 is some parameter. The integral is dominatedby the tallest of the sharp narrow peaks of the integrand. Each peak is located at a maximum ofg(x). Let x0 be the location of the biggest maximum of g between x1 and x2 , thens
2πAg(x0 )
I ≈ f (x0 )e,−Ag 00 (x0 )

(E.2)

2

d g
where g 00 (x) = dx
2.

We apply the above formula to evaluate
"
#
"
#
√
Z ∞
(C 1/3 E −1/3 − ρ0 )2(D−1/2 2E − ρ∞ )2I2 =dE exp −
exp −
,
(E.3)
a2ho
a2ho
0
√
h
i
−1/2 2E−ρ )22∞
where C = ~m d0 N 7/2 a, D = mω 2 . Here, exp − (Dis the more peaky function. Thus,2aho

we have A = 1,"#
(C 1/3 E −1/3 − ρ0 )2
f (E) = exp −,a2ho√(D−1/2 2E − ρ∞ )2
g(E) = −
,
a2ho
r

1
2 −1/2  −1/2 √0g (E) = − 2ED
2E − ρ∞ ,
aho D
r 

1
2
1 1
1 −3/2  −1/2 √00√
g (E) = − 2− E
D2E − ρ∞ ,aho D
2D E 2

(E.4)
(E.5)
(E.6)
(E.7)

95
and maximum of the integrand occurs at Em = Dρ2∞ /2. Hence,
 


√ ρ∞

I2 ≈ ~ω π
exp −aho

2d0 N 7/2 a4ho ρa2
∞
a2ho

1/3

2 
− ρ0



.


(E.8)