Hyperspherical Approach to Quantal Three-body Theory
by
Jia Wang
B.A., Tsinghua University, 2004
M.Phil., Chinese University of Hong Kong, 2006

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
2012This thesis entitled:
Hyperspherical Approach to Quantal Three-body Theory
written by Jia Wang
has been approved for the Department of Physics

Chris H. Greene

Prof. John Bohn

Prof. Eric Cornell

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
Wang, Jia (Ph.D., Physics)
Hyperspherical Approach to Quantal Three-body Theory
Thesis directed by Prof. Chris H. Greene

Hyperspherical coordinates provide a systematic way of describing three-body systems. Solving three-body Schrödinger equations in an adiabatic hyperspherical representation is the focus of
this thesis. An essentially exact solution can be found numerically by including nonadiabatic couplings using either a slow variable discretization or a traditional adiabatic method. Two different
types of three-body systems are investigated: (1) rovibrational states of the triatomic hydrogen ion
H+
3 and (2) ultracold collisions of three identical bosons.

Dedication

To my love.

Acknowledgements

First and foremost, I’d like to thank my advisor Chris H. Greene, for his patience, motivation,enthusiasm, and immense knowledge. Working with Chris is a great experience, not only becauseof his unparalleled knowledge of atomic and molecular physics, but also because he has always been
supportive of my career. He is also a nice and easygoing person, and will always be a role modelin my life.
Besides Chris, working with other people in Greene’s group is also a very nice experience.
Jose D’Incao has always been a great teacher and big brother for me. Learning and discussing
physics with him is always pleasant. I also owe many thanks to Yujun Wang. He has been very
helpful in various projects during the last year. I also want to thank all the other group members:
Michal Tarana, Chen Zhang, John Papaioannou, Nathan Morrison, Stephen Ragole, Victor Colussi,
Minghui Xu and Javier von Stecher for all the discussions about physics during our fun groupmeetings and vast spare time. Previous group members Daniel Haxton, Zach Walters, Nirav Mehta
and Seth Rittenhouse have also helped me a lot during my PhD years.
I would also like to thank the people outside the Greene group that I have gotten the pleasureof working with directly. Viatcheslav Kokoouline has given me tremendous help in the H+
3 project.
I would like to thank Richard Saykally and his group for sharing their beautiful experimental resultswith us. From the collaboration with Brett Esry, I learned how to think critically and how to bescientifically rigorous. Having the opportunity to collaborate with Eric Cornell is also a valuable
experience. I am always amazed by his great intuition of telling the underlying physics out ofcomplex phenomena.

vi
I would also like to thank my Master advisor Chi Kwong Law and Ming Chung Chu in the
Chinese University of Hong Kong, for bringing me into the beautiful world of atomic and molecularphysics and giving me a strict training in scientific research.
I also owe many thanks to Peter Ruprecht and James McKown for all their technical support,
and to Pam Leland for her administrative support.
I would like to thank all of my friends for supporting me and encouraging me with their best
wishes. A special thanks goes to Yiqing Xie for her endless support during the last five years. Last
but not least, I would like to thank my family for their unconditional support throughout my life.
I have received funding support from the National Science Foundation, the Department ofEnergy, and the Air Force Office of Scientific Research for this work.

vii

Contents

Chapter
1 Introduction

1

2 Adiabatic Hyperspherical Approach

8

2.1

Hyperspherical coordinate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

8

2.2

Adiabatic hyperspherical representation . . . . . . . . . . . . . . . . . . . . . . . . . 10

2.3

Coupling matrices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.4

Slow variable discretization (SVD) method. . . . . . . . . . . . . . . . . . . . . . . 17

2.4.1

Bound-state calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

2.4.2

Scattering calculations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18

3 Rovibrational states of H3 + and quantum-defect analysis of H3 Rydberg states

22

3.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22

3.2

Rovibrational states of H+
3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.3

3.4

3.2.1

Adiabatic representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25

3.2.2

Accuracy of the rovibrational energies of H+. . . . . . . . . . . . . . . . . . 28
3

3.2.3

Rovibrational-frame transformation . . . . . . . . . . . . . . . . . . . . . . . 28

p-wave energy levels of H3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
3.3.1

Body-frame quantum defects for p-waves . . . . . . . . . . . . . . . . . . . . 31

3.3.2

3p1 energy levels of H3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

Higher angular-momentum states . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33

viii
3.5

Recombination-pumped triatomic hydrogen infrared lasers . . . . . . . . . . . . . . . 37

3.6

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

4 Numerical study of three-body recombination

44

4.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

4.2

Adiabatic hyperspherical representation . . . . . . . . . . . . . . . . . . . . . . . . . 47

4.3

Three-body recombination rates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48

4.4

Dominant recombination pathways . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.5

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61

5 Origin of the Three-body Parameter Universality in Efimov Physics

62

5.1

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63

5.2

Van der Waals interaction and classical suppression . . . . . . . . . . . . . . . . . . . 64
5.2.1

Van der Waals interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65

5.2.2

Classical (WKB) Suppression . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

5.3

Three-body adiabatic hyperspherical representation . . . . . . . . . . . . . . . . . . . 73

5.4

Three-body parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76

5.5

5.4.1

Universality of the three-body parameter . . . . . . . . . . . . . . . . . . . . 80

5.4.2

Single channel approximation . . . . . . . . . . . . . . . . . . . . . . . . . . . 85

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89

6 Efimov physics on the positive side

92

6.1

Stückelberg minima. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93

6.2

Three-body recombination resonances associated with d-wave interactions . . . . . . 93

6.3

Three-body state associated with the d-wave dimmer . . . . . . . . . . . . . . . . . . 101

6.4

Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103

ix
Bibliography

106

Appendix
A Symmetry of R-matrix

112

B R-matrix propagation method with the traditional adiabatic method.

114

C Permutation symmetry of the basis functions

116

D Body-frame quantum defect matrix elements

119

E Effective adiabatic potentials

123

x

TablesTable
3.1

Comparison of several calculated rovibrational energy levels of H+
3 with experimental
results [1]. Only states with N + ≤ 3 for the {0, 00 } and {0, 11 } bands are shown here. 29

3.2

A comparison between several of our calculated 3p1 H3 energy levels with empiricallyfitted experimental energy levels [2, 3, 4, 5]. . . . . . . . . . . . . . . . . . . . . . . . 34

3.3

Comparison between several of our calculated 3d energy levels of H3 with experimentallydetermined energy levels[2, 3, 4, 5]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 36

3.4

Possible assignment of laser lines observed in this work. We calculate the theory
lines by using models described in this chapter. The experiment lines are chosenfrom the H2 O – He laser lines which are also observed in H2 O – He or H2 O – Ne
experiments.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

C.1 Permutation symmetry for basis functions in the different degrees of freedom. . . . . 117

xiFiguresFigure
3.1(Color online) Lowest 60 adiabatic potential curves U (R) of H+
3 with total angularmomentum N + = 1, odd parity and gI = 1. The dashed horizontal line shows thelowest rovibrational ground state of this system. . . . . . . . . . . . . . . . . . . . . 27

3.2

Comparison of experimental results with calculated nd → n0 p transitions of H3 . The
experimental laser strength is in arbitrary units with linear scale, while the present
theoretical B-coefficients are in the units of 1022 (m/Js2 ) on a logarithmic scale. The
calculated line positions have a precision of about 13 cm−1

4.1

. . . . . . . . . . . . . . 40

(Color online) Adiabatic potential curves U (R) at short distances R. Red solid lines
represent the three-body continuum channels, i.e., the initial channels for recombination processes, and black dashed lines represent the final recombination channels.

4.2

49

(Color online) Same as Fig. 4.1, but for the adiabatic potential curves U (R) at large
distances R. This figure contrasts with Fig. 4.1 in the characteristically smoothbehavior at large distances. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50

xii
4.3

(f ←i)(Color online) Partial recombination rate K3

as a function of the collision energy

E for the positive-scattering length case. The solid, dashed, dot, dash-dot and shortdashed lines indicate partial recombination rates to different recombination channels:
f = 1, 2, 3, 4, 5, respectively. The thick black lines, thinner red lines and thinnest
blue lines indicate recombination rates from different incident channels: i = 1, 2, 3.
The three green dash-dot-dot lines are proportional to E 0 , E 1 and E 4 as indicatedin the figure. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
4.4

(Color online) Same as Fig. 4.3 but for the negative scattering length case. . . . . . 54

4.5

(Color online) Branching ratio of the calculated recombination rates r3(f ←i)

as func-

tions of the collision energy E for the positive-scattering length case. The line-styles
solid, dashed, dot, dash-dot and short-dashed indicate recombination rate to different
recombination channels: f = 1, 2, 3, 4, 5, respectively. The thick black lines, thinner
red lines and thinnest blue lines indicate recombination rates from three different
incident channels: i = 1, 2, 3. The inset shows the branching ratio between the deepbound states only and it excludes the shallowest bound state (see the text for further
details). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.6

(Color online) Same as Fig. 4.5 but for the negative scattering length case. . . . . . 56

5.1

Lennard-Jones potential with the 2nd, 3rd and 4th s-wave bound state exactly at
the threshold (denoted by λ = λ∗2 , λ∗3 , λ∗4 ). The insert shows the corresponding zeroenergy wave function. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

5.2

Effective ranges at the pole as a function of the potential depths Dn∗ for vsch . . . . . 70

xiii
5.3

The red solid curve represents the scattering length, a/r0 , while the green dashedcurve represents the parameter ξpin (k → 0). Both quantities are plotted as functionsof the depth D of the two-body interaction model vsch [Eq. (5.2)], whose values
for which |a| = ∞ are indicated in the figure as Dn∗ , where n is the number ofs-wave states. The black circles, open circles, and open squares are the values ofξpin at |a| → ∞, a = 5r0 , and a = −5r0 , respectively. Their trends documents the
suppression of the ξpin as the number of bound states increases. The results for ξpin also
show higher efficiency of the suppression inside the well for |a|/r0  1. The blackdash-dotted curve shows the suppression factor 1/kL (r0 ), confirming the classicalorigin of the suppression mechanism.

5.4. . . . . . . . . . . . . . . . . . . . . . . . . . 72(a) Full energy landscape for the three-body potentials at a = ∞ for our vλa model
potential. (b) Effective diabatic potentials Wν relevant for Efimov physics for vλa withan increasingly large number of bound states (λ∗n is the value of λ that produces a =∞ and n s-wave bound states). The Wν converge to a universal potential displaying
the repulsive barrier at R ≈ 2rvdW that prevents particles access to short distances.
(c)–(e) demonstrate the suppression of the wave function inside the potential wellthrough the channel functions Φν (R; θ, ϕ) for R fixed near the minima of the Efimov
potentials in (b). (c) shows the mapping of the geometrical configurations onto thehyperangles θ and ϕ. (d) and (e) show the channel functions, where the “distance”
from the origin determines |Φν |1/2 , for two distinct cases: in (d) when there is a
substantial probability of finding two particles inside the potential well (defined by
the region containing the gray disks) and in (e) with a reduced probability — seealso our discussion in Fig. 5.5. In (d) and (e), we used the potentials vsch and vλa ,respectively, both with n = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74

xiv
5.5

Density plot of the three-body probability density |Φν (R; θ, ϕ)|2 sin 2θ that determines the three particle configuration [see Fig. 5.4 (c)] in the θ-ϕ hyperangular plane
for a fixed R (sin 2θ is the volume element). (a)–(d) show the results for an R near
the minima of the Efimov potentials in Fig. 5.4 (b) for the first four scattering length
poles of the vλa model, as indicated. (a) shows that there is a negligible probabilityof finding the particles at distances smaller than rvdW (outer dashed circle) and, ofcourse, inside the 1/r12 repulsive barrier (inner solid circle). For higher poles, i.e., as
the strength of the hard-core part of vλa potential decreases, the potential becomesdeeper and penetration into the region r < rvdW is now classically allowed. Nevertheless, (b)–(e) show that inside-the-well suppression still efficiently suppresses theprobability of finding particle pairs at distances r < rvdW , found to be in the range
2–4%. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 75

5.6

The Efimov potential obtained from the different two-body potential models used
here. The reasonably good agreement between the results obtained using models
hs[obtained by replacing the
supporting many bound states (vsch , vλa and vλb ) and vvdW

deep potential well with a hard wall but having only one (zero-energy) bound state]
supports our conclusion that the inside-the-well suppression of the wave functionis the main physical mechanism behind the universality of the three-body effective
potentials. The differences between these potentials are seen to cause differences ofa few percent in the three-body parameter. . . . . . . . . . . . . . . . . . . . . . . . 77
5.7

Fitting the Efimov resonance using a fano lineshape [Eq. (5.16)] for a system using
the two-body model potential vλa with λ = λ∗2 . The blue circles, red triangles,
and black squares are the numerically calculated |1 − Sii |2 for the three deeper
atom-dimer channels (a g-wave dimer, a d-wave dimer, and a deeper s-wave dimer,
correspondingly.) The curves are fitting results from using Eq. (5.16). . . . . . . . . 79

5.8

Efimov scenario calculated for vλb with λ = λ∗2 . Both the ground and first excitedEfimov spectra as a function of sign(a)(a/rvdW )1/4 are shown. . . . . . . . . . . . . . 81

xv
5.9

The Efimov resonance corresponding to the Efimov ground state. The black squares
are numerically calculated results, and the solid curve is the fitting result from using
Eq. (5.17).

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 82

5.10 Values for the three-body parameter (a) κ∗ and (b) a−
3b as functions of the number
n of two-body s-wave bound states for each of the potential model studied here. (c)133 Cs [6] (red: ×, +, , and ∗), 39 K [7] (magenta:
Experimental values for a−
3b forM), 7 Li [8] (blue: •) and [9, 10] (green:  and ◦), 6 Li [11, 12] (cyan: N and O) and[13, 14] (brown: H and ♦), and 85 Rb [15] (black: ). The gray region specifies ahsresults. The inset of (a) showsband where there is a ±15% deviation from the vvdW

the suppression parameter ξpin [Eq. (5.11)] which can be roughly understood as thedegree of sensitivity to nonuniversal corrections. Since ξpin is always finite — even in
the large n limit — nonuniversal effects associated with the details of the short-range
interactions can still play an important role. One example is the large deviation inκ∗ found for the vsch (n = 6) model, caused by a weakly bound g-wave state. Forn > 10 we expect κ∗ and a−
3b to lie within the range of 15% established for n ≤ 10. . 84
5.11 (a) This figure compares the energies (as characterized by the three-body parameterκ∗ ) obtained from a single channel approximation with our full calculations. The
three-body parameter κ∗ is shown for the vλa model in the single channel approximation (open triangles) as well as for our full numerical results (open circles). The
single channel approximation can be improved by imposing a simple change in the
adiabatic potentials near the barrier, as is shown in (b). There we smoothly connecths
the potential for vλa (red solid line) to the barrier obtained for vvdW
(black solidline), resulting in the potential labeled by v̄λa (green solid line). This new potential
is actually more repulsive and has energies [filled circles in (a)] that are much closer
to our full numerical calculations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86

xvi
5.12 This figure shows the three-body potentials obtained using the vλa (λ = λ∗10 ) model
supporting a total of 100 bound states. Roughly speaking, the potential of Eq. (5.18) [16](black solid line) can be seen as a diabatic potential since it passes near one of theseries of avoided crossings. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 88
5.13 Comparison between the effective potential obtained from Ref. [16] and the potentialfrom Eq. (5.20).
6.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90

The three-body recombination rate K3 as a function of scattering length a. The
black curve with square symbols shows the results from a Lennard-Jones potential
with two s-wave bound states; the red curve with circles illustrates the results forthree s-wave bound states. A Stückelberg minimum appears at about a = 28.09rvdW
for both cases; the minimum is indicated by a vertical dashed line. Two enhancement
features are also shown for both cases near the small scattering length of a ≈ 0.956
rvdW . The green dashed line is proportional to a4 , the overall scaling of K3 . . . . . . 94

6.2

The total and partial three-body recombination rate K3 as a function of scattering
length a. The black curve with square symbols shows the total recombination rate,
and the red curve with circle symbols shows the partial recombination rate for theshallowest s-wave dimer channel. The other curves shows the partial recombination
rate for deeper atom-dimer channels. The Stückelberg minimum only shows up in
the shallowest s-wave dimmer channel, but not the deeper atom-dimer channels. . . 95

6.3

The two-body s-wave scattering length a∗l when a d-wave (l = 2) dimer (black curvewith square symbols) or an i-wave (l = 6) dimer becomes bound (red curve withcircle symbols), shown as functions of the number of two-body s-wave bound states.

97

xvii
6.4

Adiabatic potential curves for a ≈ 0.977. The 2+1 channels, which corresponding toa dimer plus an atom at very large distance, are labeled by a combination of a letterand a number. The letter denotes the angular momentum quantum number of the
dimer, and the number labels the channels for the same dimer angular momentum
from low-to-high dimer energies. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 98

6.5

Total and partial three-body recombination rate as a function of two-body scattering
length a, shown on a linear scale, near a d-wave dimer that is just becoming bound.
The inset shows the same graph with a logarithmic scale for the y axis. The solid
vertical line shows the s-wave scattering length when the d-wave dimer becomesbound. The two dashed vertical lines show the two peaks of recombination rateenhancement, denoted A and B respectively. . . . . . . . . . . . . . . . . . . . . . . . 100

6.6(I) Wνν (R) as a function of hyperradius R at the K3 resonance peak A. (II) W (R)νν
as a function of hyperradius R around the K3 enhancement peak B. The inset shows
the three points where the potential curves correspond to. . . . . . . . . . . . . . . . 102

6.7

Energy of the three-body bound state associated with a d-wave dimer as a functionof scattering length a. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 104

6.8

The enhancements for the total three-body recombination rate at about a = 0.995rvdW
for Lennard-Jones potential with 2 and 3 s-wave bound states. . . . . . . . . . . . . 105

xviiiE.1 Comparison of Uν (R) (green dashed curves) and Wν (R) (red solid curves) demonstrating the importance of nonadiabatic effects introduced by Qνν (R). (a) and (b)show the results for the potential models vλa (λ = λ∗2 ) and vλa (λ = λ∗10 ) [Eq. (5.3)],
respectively, and in (c) we show the results obtained for vsch (D = D6∗ ) [Eq. (5.2)]. (a)and (b) also show the effect of the diabatization scheme used in order to prepare someof our figures in the main text (dash-dotted curves). The goal of the diabatization isto eliminate the sharp features resulting from Qνν (R) which should not contribute
substantially to the three-body observables. The case shown in (c), however, does
not allow us to easily trace the diabatic version of the potentials Wν . In this case,however, the “less” sharp features have a larger contribution due to the crossing witha three-body channel describing a collision between a g-wave molecular state and afree atom, giving rise to the anomalous n = 6 point for the vsch model in Fig. 4 ofthe main text. Although such cases are relatively infrequent in our calculations (andoccur mostly for the vsch model), they do illustrate nonuniversal effects that canaffect the three-body parameter. Nevertheless, the three-body observables obtainedfor cases like the one shown in (c) are still within the 15% range we claimed for theuniversality of the three-body parameter. . . . . . . . . . . . . . . . . . . . . . . . . 124

Chapter 1

Introduction

The studies of three-body systems have a long and interesting history. In 1687, Sir Isaac
Newton started studying the gravitational problem of three-body systems and presented someresults in his famous “Principia” [17]. One of the most important applications at that time wasfor studies of the motion of the Moon under the gravitational influence of the Earth and the Sun
as lunar theory [18]. During work to improve the accuracy of the lunar theory Henri Poincaré’s
research in the late nineteenth century led to the beginning of chaos theory. Poincaré extended the
problem from the Earth-Moon-Sun system to general three-body systems with mutual gravitational
interactions, as the “three-body problem”. Because of this study, he won the prize competition
in honor of the 60th birthday of King Oscar II of Sweden in 1889. Poincaré discovered that sucha system, under certain conditions, can exhibit chaotic behavior that is highly sensitive to initialconditions but impossible to predict in the long term. Chaotic behavior is described by the nonlinear
differential equation governing the dynamics of a classical three-body system. However, such chaotic
behavior would not exist in a classical sense for a quantal three-body system whose dynamic isgoverned by the linear Schrödinger equation. Nevertheless, the study of how chaotic classicaldynamics can be described in terms of quantum theory became an interesting question, and leads
to the beginning of quantum chaos theory in the twentieth century [19].
The studies of quantum chaos raise an interesting question: can quantum mechanics describe
the exotic behaviors of three-body systems? Quantum mechanics has been proved to be very
successful in describing two-body systems. One example is that for two particles interacting with2each other via a short-range interaction at ultracold temperatures, one single parameter, s-wave
scattering length a, is sufficient to characterize the system. This simple description is verified byusing a to determine the properties of a sufficiently dilute homogeneous Bose-Einstein condensation(BEC) [20] and compare them with experimental observations [21, 22, 23]. However, when theexperimentally achieved densities become greater and greater, physicists have realized that thedescription of the system only in terms of two-body interactions is no longer sufficient. Quantum
calculations of a three-body system are then a natural extension of two-body analogies and serve
as a key meeting point for theoretical and experimental efforts to understand few-body physics. In
this thesis, we use the hyperspherical coordinate approach to study three-body systems.
Hyperspherical coordinates have been applied successfully in several areas of theoreticalphysics ranging from nuclear physics [24, 25, 26] to atomic structure [27, 28, 29, 30] and fundamental few-body scattering [31, 32]. In a three-body system, there are nine degrees of freedom
in total. After separating out the center-of-mass motion, the remaining six degrees of motion canbe described by three Euler angles (α, β, γ) and three internal coordinates: the hyperradius R
and two hyperangles, θ and φ. In particular, the hyperradius R represents the overall size of thesystem. Some of the deepest insights into the nature of the three-body problem have emerged
from Macek’s adiabatic hyperspherical methodology [27]. In this method, the Hamiltonian of thesystem is initially diagonalized at fixed values of the hyperradius R, and the eigenvalues yield aset of 1D-coupled adiabatic potential curves that represent the energy of the system as a functionof R. The resulting eigenfunctions can be used to develop the coupling matrices between thesepotential curves. These potentials and couplings not only can be used for almost exact numerical calculation, but also help to build an intuitive understanding of the system. The potentials
can describe available reaction pathways by indicating the threshold laws and scaling laws for thecorresponding reactions. They also give information about the bound/quasi-bound states of thesystem and the excitation and decay mechanisms of these states. Furthermore, the coupling matrix
elements allow accurate calculations of the full three-body Schrödinger equation for bound-state
problems and two-body inelastic and rearrangement collisions (A + BC → AB + C), three-body3collisions (A + B + C), and photon-assisted collision processes.
The application of the hyperspherical-coordinates approach to quantal three-body problems
dates back at least as far as the pioneering work of Llwellen Hilleth Thomas in 1935 [33]. He realized
a striking quantum effect in a three-body system: when the ratio between the potential range (r0 )and s-wave scattering length (a) becomes arbitrarily small, r0 /a → 0, the ground state energy ofthe system can “collapse” to E → −∞. This effect is known as the Thomas collapse. In thiscase, all three particles collapse into an infinitely small size with an infinitely large binding energy.
This phenomenon can be understood in the hyperspherical picture. The effective hypersphericalpotential of such a system has the form of − s20 + 1/4 / 2µ3b R2 for r0  R  a, where µ3b is
the three-body reduced mass. The parameter s0 is universal, i.e., it does not depend on the form ofinteraction as long as r0 /a → 0. For three identical bosons, s0 = 1.00624. It is well known that suchpotentials can support an infinite number of bound states and that all the nearby eigen energies are
related with a geometric scaling factor En+1 /En = exp (−2π/s0 ). Clearly, this geometric scaling
factor is also universal [exp (π/s0 ) ≈ 22.7 for three identical bosons]. This universal scaling factor
implies that when n approaches −∞, infinitely tightly bound states (called Thomas’s collapse
states) exist. In real physical systems, however, the range of interparticle interactions r0 can neverbe zero. Presumably, this fact prevents a Thomas collapse from being observed. On the other hand,
there is no similar obstacle to observing a closely related quantum effect, namely the Efimov effect,as n approaches +∞. In 1970, Efimov predicted that when a two-body bound state is exactly at
the threshold, i.e., a → ∞, there is an infinite number of bound states just below the three-bodybreak-up threshold [26]. These states, called Efimov states, also obey En+1 /En = exp (−2π/s0 )
as a result of the − s20 + 1/4 / 2µ3b R2 potential for R  r0 , which is usually called the Efimov
potential.
The fact that Efimov effect is universal implies that it can exist in any systems of threeidentical bosons interacting with each other via short-range interactions. Ultracold atomic gases
are perfect systems for studying Efimov effects experimentally, because of the extraordinary degreeof control for such systems. Using techniques such as laser cooling and subsequent evaporative

4
cooling of atomic gases, the experimentalists can now reach the nano Kelvin range with highdensities (between about 1012 and 1015 cm−3 ), and finally attain BEC [21, 22, 23]. In addition,
applications of Feshbach resonances allow physicists to control the scattering length a between
ultracold atoms and study the properties of condensates. In particular, when three free particles
collide at ultracold temperatures, they can form a two-body bound state and a free particle, whichis called a three-body recombination (A + A + A → A2 + A). This recombination process normallyreleases a large amount of kinetic energy, producing atomic losses that often limit the lifetimesof Bose-Einstein condensates [34]. Theoretical studies indicate that there is an a4 scaling of the
field-free recombination rate of three identical bosons that leads to a catastrophic loss of atoms
even if a is not quite large. Three-body recombination is a process that is also closely related with
the Efimov effects. When a is much larger than the range of two-body interaction r0 but still finite,
Efimov states can cause interference and resonant effects in three-body recombination processes
when they cross thresholds. The hyperspherical approach gives a comprehensive description of
three-body recombination and a fundamental understanding of how Efimov states affect threebody recombination [35, 36].
Hyperspherical coordinate has also been applied to study triatomic molecules. As the simplest
triatomic molecule, H+
3 is an interesting system that attracts theorists to high-accuracy quantum
calculations. H+
3 also plays an important role in astrophysics since it acts as a proton donor inchemical reactions occurring in interstellar clouds [37, 38]. Furthermore, this ion also helps tocharacterize Jupiter’s atmosphere from afar [39, 40]. H+
3 is the dominant positively charged ion
in molecular hydrogen plasmas and was first identified in 1911 by J. J. Thomson with an early
form of mass spectrometry [41]. Without a stable electronic excited state and a permanent dipole
moment, H+
3 cannot be observed by electronic spectroscopy or rotational spectroscopy. Therefore,
an infrared rotation-vibration spectrum is the only mean to observe this ion. The first observation
was carried out by T. Oka in 1980 [42]. By 2012, more than 600 low-lying rovibrational statesof H+
3 had been identified. The good agreement between the experimental spectrum and a firstprinciples calculation provided a benchmark for calculations on other polyatomic molecules such

5
as water. Application of the hyperspherical approach to study H+
3 leads to a better understandingof a quantum phenomenon that once was considered mysterious and esoteric. This phenomenon isthe dissociative recombination (DR) of H+
3e− + H +
3 →





H2 + H(1.1)

 H + H + H,
which is one of the key process to understand chemical reactions in diffuse interstellar clouds[43]. However, there once was a three-orders-of-magnitude discrepancy between the theoretical andexperimental DR rate of H+
3 . This discrepancy was finally dissolved by Kokoouline and Greene
in 2003 using the hyperspherical approach [31]. The use of hyperspherical coordinates has botha practical computational advantage and a qualitative conceptual advantage in this problem. For
instance, the theory of DR is much better understood for a diatomic target than for a polyatomic
target, so the use of an adiabatic hyperspherical representation of the nuclear positions ultimately
maps polyatomic DR theory back in terms of the more familiar diatomic DR theory. In addition,
applying the hyperspherical approach to describe the coupling between an incident electron and thevibrational or dissociative degrees of freedom of H+
3 permits a natural inclusion of the Jahn-Teller
coupling and dissolves this discrepancy [31].
In this dissertation, we apply the adiabatic hyperspherical approach to investigate two different types of three-body systems: (1) rovibrational states of the triatomic hydrogen ion, H+
3 and (2)ultracold collisions of three identical bosons. Both systems provided interesting questions and rich
physics to explore. The remainder of this thesis is organized as follows. Chapter two discusses thedetails of our numercial approach, the adiabatic hyperspherical representation. The hypersphericalcoordinates used in this thesis are first introduced. Then the numerical methods for both boundstates and scattering-state calculations are elaborated.
In chapter three, we calculate the rovibrational states of H+
3 with the tools described in thelast chapter. We model H+
3 as three protons interacting with each other on a potential surface. Our
calculation gives rovibrational energies that are in good agreement with experimental results. In
addition, using accurate rovibrational-state wave functions, we apply multichannel quantum defect6
theory to studying the Rydberg energy levels of H3 , which consists of a Rydberg electron and anH+
3 ion core. The interaction between the Rydberg electron and the ion core is described via bodyframe quantum defects. We perform the rovibrational-frame transformation with rovibrational wave
functions of H+
3 to obtain laboratory-frame quantum defects that are used to calculate both theenergy levels and the mid-infrared spectrum of the H3 Rydberg states. The mid-infrared spectrum
corresponds reasonably well with the laser lines recently observed in hydrogen/rare gas discharges[44, 45], indicating that H3 is a likely candidate for the carrier of these lasing transitions. A lasing
mechanism for the population inversion is also proposed.
In chapter four, we study another type of three-body system: three-body recombination at
ultracold temperature. In this chapter, we study three-body recombination processes numericallyfor a system of three identical bosons with a much more realistic model than used in previous
studies. Our two-body model potentials support many bound states, which is a major leap incomplexity. Our study indicates that recombination into deeply bound states can be describedby the dominance of one decay pathway, resulting from the strong coupling between different
recombination channels. Moreover, the usual Wigner threshold law must be modified for excitedincident recombination channels. Three-body recombination has also been recognized as one of themost important scattering observables in which features related to the universal Efimov physics
can be manifested, which will be studied it in the next chapter.
In chapter five, problems related to Efimov physics are studied. Efimov physics predicts that
there are an infinite number of three-body bound states for a three-identical-bosons system when
the two-body s-wave scattering length a → ∞ [26]. Tuning a from infinity to finite and negative(but still much larger than the range r0 of the interactions), these three-body states disappear
into the continuum one by one at different scattering lengths a−n . Whenever an Efimov trimer
hits the threshold, a resonance in recombination is observed. It is remarkable that the valuesof a−n for two nearby resonances are related by a scaling factor of approximate by 22.7 that isuniversal, i.e., independent of the details of the short-range interaction. Hence, we only need one
three-body parameter to determine the absolute positions of all the resonances. However, this7
three-body parameter was once believed to depend on the details of two- and three-body shortrange interactions, suggesting that this parameter would not be universal. Surprisingly, recent
experiments support the idea that the three-body parameter exhibits universality. We explore theorigin of this universality. Our study shows that the three-body parameter universality emerges
because a universally effective barrier in the three-body potentials prevents the three particles fromsimultaneously getting closer to each other. Our results show limitations to this universality, as it
is more likely to occur for neutral atoms and less likely to extend to light nuclei.
Chapter six focuses on three-body collisions on the positive side. When the scattering length islarge and positive, minimums show up in the three-body recombination rate because of a destructive
interference effect. When the Efimov states hit the atom-dimer thresholds, resonances in atomdimer relaxation can also be observed. The universality of three-body parameter represented by theminimum position is also found to be universal for cold atoms. In addition, a universal three-bodystate attached to the d-wave two-body state is found in our model. The d-wave state crosses thethreshold and becomes bound at a universal s-wave scattering length a∗d . Near a∗d , two enhancementpeaks in three-body recombination can be found. The positions of the two peaks are also quite
universal. The one larger than a∗d corresponds to the three-body state attaching to the d-wave dimer,
and the one smaller than a∗d corresponds to the competition between a new d-wave threshold andthe a4 power law.

Chapter 2

Adiabatic Hyperspherical Approach

Hyperspherical coordinates have a long and distinguished history in atomic, molecular andnuclear physics. This chapter discusses the details of the adiabatic hyperspherical approach, which
is a numerical method we found to be suitable for calculations of three-body systems. Section
2.1 introduces the hyperradial coordinates used in this thesis. Section 2.2 sets up the adiabatic
hyperspherical approach. Some particularly important quantities, called the coupling matrices, arestudied in section 2.3. Finally, section 2.4 presents applications of the adiabatic hyperspherical approach to three-body problems, including both bound-state calculations and scattering studies.

2.1

Hyperspherical coordinate
There are numerous conventions for defining hyperspherical coordinates. The modified ver-sion of Whitten-Smith’s democratic coordinate described in this section is one of the most convenient conventions for our study. After separating the center-of-mass motion, the six remaining
degrees of freedom are described by three Euler angles (α, β,γ), two hyperangles (θ and φ), and one
hyperradius R. In this convention, the hyperangles and the hyperradius can be best constructedin two stages, as shown in Ref. [46]. We first introduce the mass-scaled Jacobi coordinates,1[~r2 − ~r1 ] ,d12
m1~r1 + m2~r2= d12 ~r3 −,m1 + m2

ρ
~1 =
ρ
~2

(2.1)

9
where ~ri is the position of particle i with mass mi , and the mass-weighting factor d is given byd2ij=µ23b =

mk (mi + mj ),µ3b (mi + mj + mk )m1 m2 m3.
m1 + m2 + m3

(2.2)

Here, µ3b is the three-body reduced mass. In the case of three identical particles, d12 = d23 =√d31 = 21/2 /31/4 and µ3b = m/ 3, where m1 = m2 = m3 = m is the mass of a single particle. The
hyperradius R describing the overall size of the system can be defined as
R2 = ρ21 + ρ22 ,(2.3)

where R ∈ [0, ∞). The two hyperangles (θ, ϕ) describing the shape of the system are defined by(~ρ1 )x = R cos (θ/2 − π/4) sin (ϕ/2 + π/6) ,(~ρ1 )y = R sin (θ/2 − π/4) cos (ϕ/2 + π/6) ,(~ρ1 )z = 0,

(2.4)

(~
ρ2 )x = R cos (θ/2 − π/4) sin (ϕ/2 + π/6) ,(~ρ2 )y = −R sin (θ/2 − π/4) sin (ϕ/2 + π/6) ,(~ρ2 )z = 0.
Here, x, y, z are body frame axes that refer to the principal axes of the body frame whose orientationis specified by the Euler angles (α, β, γ): the z axis is parallel to ρ
~1 ×~
ρ2 , and the x axis is associatedwith the smallest moment of inertia. The hyperangle θ spans the ranges [0, π/2] and the hyperangle
ϕ spans the range [0, 2π), after requring the wave function to be single valued. The hyperangle ϕcan be further restricted to a smaller range if two or more particles are indistinguishable [47].
Usually, the three-body interaction V (R, θ, ϕ) can be expressed as a function of the interparticle distances V (R, θ, ϕ) = Ve (r12 , r23 , r31 ). For three identical particles, the interparticle distances
rij can be expressed in terms of hyperspherical coordinates [47, 48, 49, 50, 51], where
r12 = 3−1/4 R [1 + sin θ sin (ϕ − π/6)]1/2 ,10
r23 = 3−1/4 R [1 + sin θ sin (ϕ − 5π/6)]1/2 ,(2.5)

r31 = 3−1/4 R [1 + sin θ sin (ϕ + π/2)]1/2 .
The three-body Schrödinger equation in this hyperspherical representation can then be written as

−




1 ∂ 5 ∂~2
Λ2 (θ, ϕ)
+
V
(R,
θ,
ϕ)
Ψ = EΨ,R−
2µ3b R5 ∂R ∂R
R2

(2.6)

where Λ2 is the “grand angular-momentum operator” [52] defined asΛ2= Tθ + TϕC + Tr ,2µ3b R2

(2.7)

where
∂
∂
2
sin 2θ ,µ3b R2 sin 2θ ∂θ
∂θ


2
∂
Jz 2
TϕC =
i
− cos θ,2
µ3b R2 sin2 θ ∂ϕTθ = −

(2.8)
(2.9)

and
Tr =

Jy2
Jz2
Jx2
+
+
.
µ3b R2 (1 − sin θ) µ3b R2 (1 + sin θ) 2µ3b R2(2.10)The operators (Jx , Jy , Jz ) are the body-frame components of the total angular momentum of thesystem. One convenient transformation is to introduce a rescaled wave function ψE = R5/2 Ψ. The
Schrödinger equation for ψE then becomes,


 2


~2
∂
Λ2 (θ, ϕ) + 15/4−−
+ V (R, θ, ϕ) ψE = EψE .
2µ3b ∂R2
R2(2.11)The volume element relevant to integrals over |ψE |2 is 2dR sin 2θdθdϕdα sin βdβdγ.

2.2Adiabatic hyperspherical representation
Solving the Schödinger equation [Eq.(2.11)] directly should in principle provide accurate

results for three-body problems, but it would require extensive computational time and memory to
diagonalize the full Hamiltonian matrix. Instead, we break the problem into two steps: first solve
the hyperangular Schödinger equation in the adiabatic representation, and then later include the

11nonadiabatic coupling. The adiabatic hyperspherical representation is a numerically efficient andstable way to solve the three-body Schrödinger equation [Eq. (2.11)]. It is currently a standard
method. The first step in solving the three-body Schrödinger equation is to find the adiabatic
potentials and channel functions, which are defined as solutions of the adiabatic eigenvalue equations
Had (R, Ω) Φν (R; Ω) = Uν (R) Φν (R; Ω) ,(2.12)

whose solutions depend parametrically on R. Here, Ω denotes the Euler angles and the two hyperspherical angles. The adiabatic Hamiltonian, containing all angular dependence and interactions,is defined asHad (R, Ω) =


~2 Λ2
15~2
+
+
V
(R,
θ,
ϕ)
.
2µ3b R2 8µ3b R2(2.13)

Therefore, the adiabatic potentials and nonadiabatic couplings obtained by solving Eq. (2.12) for
fixed values of R contain all the correlations relevant to this problem. For each R, the set ofΦν (Ω; R) is orthogonal,ZdΩΦµ (R; Ω)∗ Φν (R; Ω) = δµν ,(2.14)and completeXΦτ (R; Ω) Φτ R; Ω0∗= δ Ω − Ω0 .

(2.15)

τIn practice, calculating all the channel functions is time consuming and impractical. However,numerical studies show that only a small number of channels are needed as a truncated set of thebasis to expand the whole wave function, e.g.,ψE (R, Ω) =Nc
XFνE (R) Φν (Ω; R),(2.16)

ν=1

where Ω = {α, β, γ, θ, ϕ} denotes the Euler angles and the two hyperspherical angles, Nc is thenumber of channels adopted. Insertion of Eq. (2.16) into Eq. (2.11) leads to a set of coupled
one-dimensional equations:

−




~2 d2
~2 X
d
E
E
+
U
(R)
−
E
F
(R)
−
2P
(R)
+
Q
(R)
Fµν
0 (R) = 0,
ν
νµ
νµ
νν 0
2µ3b dR22µ3b µdR(2.17)

12
where P and Q are the coupling matrices defined below, and ν 0 denotes the ν 0 −th independent
solution. (Hereafter, unless otherwise specified, we use an underline to denote the matrix form, e.g.,
P denotes a matrix with matrix element Pνµ . ) Numerically, if there are Nc coupled equations,then there are Nc independent solutions in general, before imposing any boundary condition at
R → ∞. For a scattering calculation, these solutions can be used to compute the R-matrix (R);this R is a fundamental quantity that can be subsequently used to determine the scattering matrix
S, which is the main goal of the scattering study [see Eqs. (2.21) and (2.22) below]. As usual, theR-matrix R is defined for some large, fixed radius R ashi−1,
R (R) = F (R) Fe (R)(2.18)

where matrices F and Fe are given in terms of the solutions of Eqs. (2.11) and (2.12) by:ZFνν 0 (R) =ZFeνν 0 (R) =dΩΦν (Ω; R)∗ ψν 0 (Ω, R) ,dΩΦν (Ω; R)∗∂ψν 0 (Ω, R) .
∂R

(2.19)(2.20)

Once we have the R-matrix at large distances, the physical scattering matrix S (and its close relative, the reaction matrix K) can be simply determined by applying asymptotic boundary conditions[53], i.e.,−1
K = f − f 0R g − g0R,(2.21)

S = (1 + iK) (1 − iK)−1 ,(2.22)

where f , f 0 , g and g 0 are diagonal matrixes whose matrix elements are the energy-normalized
asymptotic solutions fν , gν and their derivatives fν0 , gν0 , respectively. fν and gν are given in terms ofspherical Bessel functions: fν (R) = (2µ3b kν /π)1/2 Rjlν (kν R), gν (R) = (2µ3b kν /π)1/2 Rnlν (kν R),
where kν and lν are determined by the asymptotic behavior of the potential in Eqs. (2.28-2.29).1
In general, R is symmetric, which guarantees that K is symmetric and S is unitary (see Appendix
A).
1

For the f ’th recombination channel, lf is given by Eq. (4.5), and kf =√break-up channel, li = λi + 3/2, and ki = 2µ3b E.

r



(f )
2µ3b E − E2b . For the i’th three-body13
For a bound state calculation, we need to apply the boundary condition that the wave function
vanishes at R → ∞. The eigen-energy will be quantized, and only one solution (denoted as FνEi )survives for a given eigen energy Ei .

2.3

Coupling matrices
In Eq. (2.17), the coupling matrices describing nonadiabatic coupling are critical for a quan-titative calculation. They are defined as
ZPνµ (R) =Z
Qνµ (R) =dΩΦν (R; Ω)∗∂Φµ (R; Ω) ,∂R

(2.23)dΩΦν (R; Ω)∗∂2
Φµ (R; Ω) .
∂R2

(2.24)

In practice, only the2Pνµ(R) = −ZdΩ∂∂Φν (R; Ω)∗Φµ (R; Ω) ,∂R
∂R

(2.25)

component of Qνµ (see the Appendix B) is needed to solve the coupled equations. The relationdbetween Q and P is given by dR
P = −P 2 + Q. From the definition of the P and Q matrices, it is
2 = P2 ,easy to see that the coupling matrices have the following properties: Pνµ = −Pµν and Pνµµν
2 . We usually define a matrix as the addition of thewhich leads to Pνν = 0, and Qνν = −Pνν

coupling to the potential, with the matrix elements being~2dWνµ (R) = Uν (R) δνµ −2Pνµ (R)+ Qνµ (R) .
2µ3bdR(2.26)The diagonal terms of this matrix
2eν (R) ≡ Wνν (R) = Uν (R) − ~ P 2 (R) ,U
2µ3b νν(2.27)are called “adiabatic potentials with diagonal correction” or “effective potentials”. The effective
potentials are usually more physical than the adiabatic potentials. They describe the system better.
For example, the effective potential gives physical asymptotic behaviors of the system at large R. In
the case of three particles with only short-range interactions, the zero energy of the system can be14defined as the energy of three stationary particles that are infinitely far away from each other (thethree-body break up threshold). The effective potentials can then be classified into two categories.
The potential curves converging to asymptotic limiting values below the three-body thresholdat very large R are called recombination channels. These recombination channels have asymptotic
behavior at large R as,lν (lν + 1)(ν)
eν (R) R→∞U≈+ E2b ,
22µ3b R(2.28)(ν)

where E2b is the two-body bound-state (dimer) energies, and lν is the corresponding angularmomentum of the third particle relative to the dimer. One can easily see that these recombination
channels have a strong “dimer plus atom” character; they are sometimes called 2 + 1 channels.
The potential curves above the three-body threshold at very large R are called the three-bodybreak-up channels. Their asymptotic behavior is described byλν (λν + 4) + 15/4eν (R) R→∞,U≈
2µ3b R2(2.29)

where λν (λν + 4) is the eigenvalue of the grand angular momentum operator Λ2 . In the large R
limit, the corresponding eigenstates will be hyperspherical harmonics. Therefore, these channels
have a strong “atom plus atom plus atom” character, and are sometimes called 1 + 1 + 1 channels.
In addition, the asymptotic values of the coupling matrix elements in the large R limit alsoobey some simple power laws of R [54, 55]. Therefore, we usually calculate the values of thesematrix elements at large R then fit them with power law expansions, and extrapolate them to aeven larger distance, if desired.
A traditional method for calculating the coupling matrices is to apply a simple differencing
scheme for the derivative of Φµ (R; Ω), i.e.,Φµ (R + ∆R; Ω) − Φµ (R − ∆R; Ω)∂Φµ (R; Ω) ≈
.
∂R
2∆R

(2.30)

However, this scheme is only accurate up to the first order of ∆R. In addition, the value chosenfor ∆R in a realistic numerical calculation can sometimes be tricky. When Φµ (R; Ω) changes very
rapidly, e.g., near a sharp avoided crossing, we need to choose a very small step size ∆R. In

15contrast, when Φµ (R; Ω) changes very slowly, e.g., at very large distances R, we need to choose a
relatively larger step size ∆R, or else Φµ (R + ∆R; Ω) − Φµ (R − ∆R; Ω) would be too small, andthe accuracy would be limited by the machine precision.
One way to improve the accuracy is to apply the Hellmann–Feynman theorem. The Hellmann–
Feynman theorem can give us analytical formulas for the coupling matrices if we know the derivative∂of the adiabatic Hamiltonian ∂R
Had from the following derivation. First, taking the derivative ofboth sides of Eq. (2.12) leads to


∂
∂
∂
[Had (R, Ω) − Uν (R)]Φν (R; Ω) = −Had (R, Ω) −Uν (R) Φν (R; Ω) .
∂R
∂R
∂R

(2.31)

Next, multiplying Φµ (R; Ω)∗ on both sides of Eq. (2.31) and integrating over Ω gives
Z
Pµν =


 ∂
R
Had (R, Ω) Φν (R; Ω)dΩΦµ (R; Ω)∗ ∂R∂dΩΦµ (R; Ω)Φν (R; Ω) = −,∂R[Uµ (R) − Uν (R)]∗

(2.32)

where µ 6= ν, and∂Uν (R) =∂RZ∗dΩΦν (R; Ω)∂
Had (R, Ω) Φν (R; Ω)∂R(2.33)

after some manipulation of algebra. The Hellmann–Feynman theorem is believed to be numericallyexact (assuming that the numerical basis expansion for Ω is complete enough). The matrix elements
for P 2 can be obtained by
2
Pµν
=

Nc
X

Pµτ Pτ ν ,(2.34)

τ =12
where Nc is the number of channels. However, numerical studies show that the convergence of Pµν

with respect to number of channels is very slow, making this method impractical.
∂
We now introduce a new method to calculate ∂RΦµ (R; Ω). It is numerically exact (again,

assuming that the numerical basis expansion for Ω are complete enough). The first hint of the∂derivations of this method is that Eq. (2.31) seems plausible to be directly solved for ∂RΦµ (R; Ω)

by


∂
∂
∂
−1
Φν (R; Ω) = − [Had (R, Ω) − Uν (R)]
Had (R, Ω) −Uν (R) Φν (R; Ω) .
∂R
∂R
∂R

(2.35)

16
However, this solution is forbidden since Had (R, Ω) − Uν (R) is singular:|Had (R, Ω) − Uν (R)| = 0(2.36)

meaning that Had (R, Ω) − Uν (R) is not invertible. The singularity can also be understood fromthe fact that the equation
∂∂[Had (R, Ω) − Uν (R)] χν (R; Ω) = −Had (R, Ω) −Uν (R) Φν (R; Ω)∂R∂R(2.37)

does not have a unique solution, χν (R; Ω). In fact, any functions with the form ofχν (R; Ω) =∂Φν (R; Ω) + cΦν (R; Ω) ,∂R

(2.38)(where c is an arbitrary number) can be a solution of Eq. (2.37). The singularity of matrix
Had (R, Ω) − Uν (R) can be removed by considering the additional condition thatZdΩΦν (R, Ω)∗∂Φν (R, Ω) = 0,∂R

(2.39)

which can be derived from the normalization condition of Φν (R, Ω), as shown in Ref. [56]. Nevertheless, our numerical studies show that even without removing the singularity, applying numericalsolver packages such as “Linear Algebra PACKage” (LAPACK) [57] or PARDISO [58] to solve
Eq. (2.37) directly can still give an accurate solution χν (R; Ω) in the form of Eq. (2.38) with an
unknown c. And once we have the numerical solution χν (R; Ω), c can be calculated byZc=dΩΦν (R; Ω)∗ χν (R; Ω) .

(2.40)

Finally, the derivative of Φν (R; Ω) can be written as∂Φν (R; Ω) = χν (R; Ω) − cΦν (R; Ω) ,∂R

(2.41)

which can be inserted into Eq. (2.23) and Eq. (2.25) for the coupling matrices. The P matrices
obtained in this way are found to be numerically the same as the one calculated from the Hellmann–∂
Feynman theorem up to machine precision, proving that our ∂RΦν (R; Ω) are numerically accurate.

Therefore, this method can give numerically very accurate coupling matrices.

17

2.4

Slow variable discretization (SVD) method
The traditional method using Eq. (2.17) works well, however, only when Pνµ (R) and Qνµ (R)are smooth functions of R. In this case, the P and Q (actually P 2 ) matrices can be calculated ona sparse grid and then interpolated and/or extrapolated on a much denser grid and even larger
distances. Clearly, this scheme suffers from tremendous numerical difficulties arising from sharp
nonadiabatic avoided crossings. In that case, the SVD approach offers a much more stable andaccurate approach for solving Eq. (2.11). One key ingredient for implementing the SVD approachis the use of the discrete variable representation (DVR) [59, 60]. Our DVR basis functions πi (R)are defined by the Gauss-Lobatto quadrature points xi and weights wi [61]. This quadratureapproximates integrals of a function g (x) as
Z 1
NX
g (x)dx ∼g (xi )wi .
=−1(2.42)

i=1

After scaling the quadrature points and weights, the above equation is generalized to treat definite
integrals over an arbitrary interval R ∈ [a1 , a2 ]:
Z a2
N
X
∼
g (R) dR =g (Ri )w
ei ,a1

(2.43)

i=1

where
w
ei =

a2 + a1
a2 − a1
a2 − a1
wi , Ri =
xi +
.
2
2
2

(2.44)

Equation (2.43) is exact for polynomials whose degree is less than or equal to 2N − 1. We construct
the DVR basis functions asr
πi (R) =N

1 Y R − Rj,wei
Ri − Rj(2.45)

j6=i

which have the important property thatrπi (Rj ) =

1
δij .
w
ei

Hence, over an interval R ∈ [a1 , a2 ], the DVR approximation based on quadrature gives
Z a2
πi (R) H (R) πj (R) dR ∼= H (Ri ) δij
a1

for matrix elements of any function H (R), which is usually called the DVR approximation.

(2.46)(2.47)

18
2.4.1

Bound-state calculations

This subsection discusses the implementation of SVD in bound-state calculations. The solution ψν 0 is expanded in the radial DVR basis πn (R) and in hyperangles in terms of the adiabatic
hyperspherical channel functions as
ψEi (R, Ω) =Xcinν πn (R) Φν (Ω; Rn ).

(2.48)

ν,nIt is now possible to rewrite Eq.(2.11) under the DVR approximation Eq. (2.47) asX

Tnn0 Onν,n0 µ cin0 µ + [Uν (Rn ) − Ei ] cinν = 0,(2.49)

1 ∂2dRπn (R) −πn0 (R) dR,
2µ ∂R2(2.50)

n,µ

where
Z
Tnn0 =

and Onν,n0 µ is the overlap matrix given byZOnν,n0 µ =dΩΦν (Ω; Rn )∗ Φµ (Ω; Rn0 ) .

(2.51)

Finally, Eq.(2.49) is solved for the expansion coefficients cinν and Ei .
2.4.2

Scattering calculations

In our scattering calculations, the R-matrix propagation method is combined with the SVD
approach (following the logic of Ref. [62]) and uses the DVR basis given by Eq. (2.45). For agiven R-matrix [Eq. (2.18)] at R = a1 , one uses the R-matrix propagation method to calculate
the corresponding R-matrix at another point R = a2 , as follows. The solution ψν 0 is expanded in
the radial DVR basis πj (R) and in hyperangles in terms of the adiabatic hyperspherical channel
functions asψν 0 (R, Ω) =X

cjµ,ν 0 πj (R) Φµ (Ω; Rj ),(2.52)

jµwhere Φν (Ω; Rj ) is the ν−th hyperspherical adiabatic channel function calculated at R = Rj .
Substituting Eq. (2.52) into Eq. (2.19) yields the values of the matrix elements of Fνν 0 andFeνν 0 at the R = a1 and R = a2 boundaries in terms of the coefficients of Eq. (2.52):

19

Fνν 0 (a1 ) =X

cjν,ν 0 πj (a1 ),(2.53)

cjν,ν 0 πj (a2 ),(2.54)

1j 0
cjµ,ν 0 Oνµπj (a1 ),(2.55)

Nj 0
cjµ,ν 0 Oνµπj (a2 ),(2.56)dΩΦν (Ω; Rj )∗ Φµ (Ω; Ri ) .

(2.57)

jFνν 0 (a2 ) =XFeνν 0 (a1 ) =

X

j

jµ

Feνν 0 (a2 ) =

X
jµ

ij
where Oνµare the overlap matrix elements, and

ji
Oνµ
=

Z

Note that the determination of F and F̃ according to the above expressions only depends onderivatives of the well-behaved DVR basis [πj0 (R)] . Therefore, this approach is much better suitedfor handling the complex structure of avoided crossings present in systems.
Over an interval R ∈ [a1 , a2 ], the DVR approximation gives
Z a2
πi (R) Had (R, Ω) πj (R) dR ≈ Had (Ri , Ω) δij .

(2.58)

a1

Expansion of the Schrödinger equation in the same numerical basis functions as in Eq. (2.52) andintegration by parts yields the equation for the expansion coefficients cjµ,ν 0 (in vector notation, ~cν 0 )
as
h

i
e
H − E ~cν 0 = L~cν 0 ,(2.59)

hi−1
e −E
L~cν 0 .
~cν 0 = H

(2.60)or, equivalently,

e and L are given byHere, the matrix elements of He iν,jµ =HLiν,jµ =Z a2
1dπi (R) dπj (R)ijdR Oνµ+ Uν (Ri ) δνµ δij ,2µ3b a1dRdR
dπj (R) ij a2
1
πi (R)
Oνµ
.
2µ3b
dR
a1

(2.61)(2.62)

20
e over the range [a1 , a2 ] gives,
Diagonalizing H
e x~0 = εn δn,n0 ,x~n T Hn(2.63)and the completeness relation of ~xn ,X~xn ~xTn = 1,

(2.64)nwhere 1 is an identity matrix. Equation (2.60) is then rewritten as
i−1 XhX ~xn ~xTne −E
~xn ~xTn L~cν 0 =~cν 0 = H
L~cν 0 .
ε
−
E
n
n
n

(2.65)

Substitution of the matrix elements of L from Eq. (2.62) and insertion of the definition of Fνν 0 andFeνν 0 at a1 and a2 into Eq. (2.53) finally givesFνν 0 (a1 ) =(n)
X u(n)
ν (a1 ) uµ (a2 )nµFνν 0 (a2 ) =2µ3b (εn − E)nµX uν(n) (a1 ) u(n)µ (a1 )nµ(n)
X u(n)
ν (a2 ) uµ (a2 )2µ3b (εn − E)Feµν 0 (a2 ) −Feµν 0 (a2 ) −2µ3b (εn − E)

X uν(n) (a2 ) u(n)µ (a1 )
nµ2µ3b (εn − E)Feµν 0 (a1 ),(2.66)

Feµν 0 (a1 ),(2.67)

where,u(n)ν (R) =

X

xjν,n πj (R),

(2.68)jand xjν,n are elements of the vector ~xn .
Our next step introduces the following matrices(R11 )νµ =(n)
X u(n)
ν (a1 ) uµ (a1 )nµ(R12 )νµ =(n)
X u(n)
ν (a1 ) uµ (a2 )nµ(R21 )νµ =

2µ3b (εn − E)(n)
X u(n)
ν (a2 ) uµ (a1 )nµ(R22 )νµ =

2µ3b (εn − E)

2µ3b (εn − E)(n)
X u(n)
ν (a2 ) uµ (a2 )nµ2µ3b (εn − E),(2.69),(2.70),(2.71).

(2.72)

After some manipulation, the matrix equation is finally obtained that determines the R-matrix
propagation from a1 to a2 :R (a2 ) = R22 − R21 [R11 + R (a1 )]−1 R12 .

(2.73)

21
ji
In the SVD method, the overlap matrix Oνµrequires us to calculate the channel functionsΦν (Ω; Rj ) at every grid point Rj , which can be very memory demanding if one needs to perform
calculations over a broad range of R. At large distances, therefore, we apply the traditional adiabatic
approach combined with the R-matrix propagation method. In the traditional adiabatic method,
the P and Q matrixes can be calculated on a sparse grid, and then interpolated and/or extrapolated
on a much denser grid and larger distances. This strategy makes the calculation faster and it also
requires less memory. The main difference between the traditional adiabatic approach and the SVD
method is the use of a different three-body numerical basis. The details of this traditional approachand its connection with the SVD method are discussed in Appendix B.

Chapter 3
Rovibrational states of H3 + and quantum-defect analysis of H3 Rydberg states

In this chapter, the hyperspherical approach is used to study rovibrational states of triatomic
hydrogen ion (H+3 ). These rovibrational states have important applications in a multichannel quantum defect theory (MQDT) analysis of Rydberg energy levels of the triatomic hydrogen molecule(H3 ). In MQDT, interactions between the Rydberg electron and the ion core H+
3 are describedby quantum defects. We extract the body-frame p-wave quantum defects from highly accurate abinitio electronic potential surfaces and calculate the quantum defects of higher angular momentum states in a long-range multipole potential model. Laboratory-frame quantum defect matrices
emerge from a rovibrational-frame transformation carried out with accurate rovibrational states ofH+
3 . Finally, the laboratory-frame quantum defects are used to calculate 3p and 3d Rydberg energy
levels for the fundamental neutral triatomic molecule H3 . In addition, calculations of radiativetransitions for higher Rydberg states give explanations for a recent experiment. In this experiment, mid-infrared laser lines observed in hydrogen/rare gas discharges are assigned to three-bodyrecombination processes involving an electron, a rare gas (He or Ne) atom, and the H+
3 [44, 45]. A
mechanism for the population inversion is proposed.
Note that the material in this chapter has been published in Ref. [63] and Ref. [45].

3.1

Introduction
The triatomic hydrogen molecule (H3 ) plays an important role in astrophysics because its

cation H+
3 acts as a proton donor in chemical reactions occurring in interstellar clouds. As the

23
simplest triatomic neutral molecule, H3 also attracts fundamental interest. Ever since its emission
spectra were first observed by G. Herzberg in the 1980s [2, 3, 4, 5], H3 has been studied extensively.
Herzberg and co-workers measured infrared and visible emission spectra of H3 in discharges throughhydrogen and assigned them to Rydberg-Rydberg transitions between n = 2 and n = 3 electronic
states using empirical fits [2, 3, 4, 5]. Helm and co-workers investigated the higher Rydberg statesand ionization potentials of H3 by analyzing the photoabsorption spectrum [64, 65]. In 2003,
building on previous work of Schneider, Orel and Suzor-Weiner [66], it was shown [31, 67] thatintermediate Rydberg states of H3 play an important role in the dissociative recombination (DR)− → H → H + H or H + H + H. Prior to the study of Ref. [31, 67], theprocess, H+
3
2
3 + e

large discrepancy between the DR rate determined by experiment and previous theory had notbeen resolved. Refs. [31, 67] found that Jahn-Teller effects in H3 neglected in previous theoreticalstudies couple the electronic and nuclear degrees of freedom and generate a relatively high DR rate
via intermediate p-wave Rydberg-state pathways. A recent alternative formulation developed byJungen and Pratt provides supporting evidence for this interpretation [68]. Vervloet and Watson
improved both the experimental techniques and empirical fits and reinvestigated the low Rydberg
states that G. Herzberg had observed [69]. Here we undertake an analysis of the Rydberg stateswith ab initio theory. One of the most successful techniques in treating Rydberg states by abinitio theory is MQDT [70, 71]. Earlier studies [31, 72, 73, 74] have utilized MQDT to successfully
describe the DR process.
The application of MQDT to study molecular Rydberg energy levels treats the H3 molecule
as a Rydberg electron attached to a H+
3 ion. The interaction between the Rydberg electron and theion core is described through a smooth reaction matrix K or quantum defect matrix µ. K and µ are
simply related, e.g., for a single-channel, K = tan (πµ). We extract a body-frame reaction matrix
from ab initio electronic potential surfaces for p-wave Rydberg states and calculate the body-frame
reaction matrix for higher angular momentum states (l > 1) by using the long-range multipole
potential model. For the higher angular-momentum states, we neglect short-range interactions
due to the nonpenetrating nature of the high l states. Here l denotes the quantum number of

24
the Rydberg electron orbital angular momentum. We then construct the total laboratory-frame
reaction matrix K through a rovibrational frame transformation, obtaining
Kii0 =

X

hi| αiK̃αα0 α0 i0 .

(3.1)

α,α0

Here, Kii0 is the laboratory-frame reaction matrix element between the laboratory-frame eigenchannels |ii and |i0 i, and K̃αα0 is the body-frame reaction matrix element between body-frame
eigenchannels |αi and |α0 i. The rovibrational frame transformation is specified by the unitarytransformation Uiα = hi |αi.
The process of constructing the rovibrational transformation is similar to that describedin Ref. [31] and is based on the rovibrational wave functions of H+
3 . To calculate them, there
are two important approximations adopted in Ref. [31], the rigid rotator approximation and theadiabatic hyperspherical approximation. The nonadiabatic coupling between different adiabatic
hyperspherical channels was included in later studies by using the slow variable discretization(SVD) approach in Ref. [75, 76]. In the present study, we abandon the rigid rotator approximationand consider the Coriolis interaction. In this way, we obtain very accurate rovibrational energy
levels and wave functions of H+
3 that allow us to construct the rovibrational transformation.
After the rovibrational transformation described by Eq. (3.1) is carried out, we obtain thelaboratory-frame K matrix and calculate the eigenenergies E of the H3 molecule by solving thesecular equation [53]det |tan (πν) + K| = 0,
where ν is a diagonal matrix with elements νii = 1/

(3.2)

p
2 (Ei − E). Here Ei denotes the ith rovibra-

tional energy level.
The remainder of the chapter is organized as follows. Section 3.2 describes the detailed
calculation of rovibrational states of H+
3 and shows how to use them to construct the rovibrationaltransformation. Section 3.3 describes the calculation of the p-wave energy levels of H3 using ab
initio quantum defects. Section 3.4 discusses the long-range multipole potential model for higher
angular momentum Rydberg states, and Section 3.5 gives our conclusions.

25

Rovibrational states of H+
3

3.2

e
In the cation H+
3 , three protons interact with each other under a potential surface V (r12 , r23 , r31 ).
The potential surface was created by Refs. [77, 78, 79] and is sub-micro-hartree accurate. With
this potential surface, the three-body Schrödinger equation is solved in the adiabatic hypersphericalrepresentation.
In a previous study [31], two approximations are adopted: adiabatic hyperspherical approximation and rigid rotator approximation. The adiabatic hyperspherical approximation refers to neglecting coupling matrices in Eq. (2.11). Under this approximation, all the channels are decoupledwith each other, and can be treated as a 1-D single channel problem. Rigid rotator approximationrefers to neglecting the Coriolis interaction that couples rotational and vibrational dynamics. The
vibrational energy are first calculated for J = 0, where the grand angular momentum operator canbe simplified as,
Λ2 =

−4 ∂
∂
4 ∂2
sin 2θ
−.
sin 2θ ∂θ∂θ sin2 θ ∂ϕ2(3.3)

The rotational energy are latter included approximately as the rotational energy of a rigid rotator,
which only depends on the rotational quantum numbers and the moment of inertia.
Here, we improve both of these approximations: (1) using the SVD method to include thenonadiabatic couplings between different adiabatic channels, and (2) including the Coriolis interaction that couples rotational angular momentum with vibrational angular momentum by implementing the full “grand angular momentum operator”.

3.2.1

Adiabatic representation

The first step is to solve Eq. (2.12) numerically for the adiabatic potentials and couplings.
+

+

m gI
In order to do so, we expand Φν (Ω; R) in a set of basis ΦN, such that
jm2 K +Φν (Ω; R) =X
jm2 K +

(ν)

+

+

m gI
ajm2 K + (R) ΦN,
jm2 K +

(3.4)

26
where j, m2 , N + , K + and gI are quantum numbers labeling basis functions in different degrees of
+

+

m gI
freedom. ΦNsatisfies the permutation symmetry of a three-fermion system, e.g.,
jm2 K +
+

+

+

+

+

+

m gI
m gI
P12 ΦN= −ΦN,
jm2 K +
jm2 K +

(3.5)

and
+

+

m gI
m gI
AΦN
= ΦN
,
jm2 K +
jm2 K +

(3.6)

A = 1 − P12 − P 23 − P31 + P12 P31 + P12 P23 .

(3.7)

where

+

+

m gI
is given in Appendix C. Expanding in these basis
The explicit form of the basis functions ΦN
jm2 K +

functions, Eq. (2.12) can be solved for the adiabatic potentials Uν (R) and the channel functionsΦν (Ω; R).
For the purpose of illustration, Fig. 3.1 presents the lowest 60 hyperspherical adiabatic+
potentials of H+
3 . The total angular momentum of the system here is N = 1, with odd parity andgI = 1 (spin para state). Near 3 ≤ R ≤ 5, these hyperspherical potentials show a series of avoid
crossing, implying the existence of important nonadiabatic couplings between different channels.
These couplings are included through the SVD method. The dashed line shows the position of theground rovibrational level of the ion. This is in fact the lowest possible rovibrational state of H+
3
shown in table 3.1 below.
Once Uν (R) and Φν (Ω; R) are obtained, the expansion coefficients cinν and Ei can be solved
from Eq.(2.49). The total rovibrational wave function is therefore given by++IΠψEi (R, Ω) ≡ ΨvN +gm+ (R, Ω) =Xcinν πn (R)

nν

X

(ν)

+

+

m gI
ajm2 K + (Rn ) ΦN,
jm2 K +

(3.8)

jm2 K +

corresponding to the eigenenergy Ei . Here i is the set of good quantum numbers N + m+ v + gI Π+ ,where v + denotes the vibrational quantum numbers.

27

0.30

Energy (a.u.)

0.25

0.20

0.15

0.10

0.05

0.00
1

2

3

4

5

6

R(a.u.)
Figure 3.1: (Color online) Lowest 60 adiabatic potential curves U (R) of H+
3 with total angular momentum N + = 1, odd parity and gI = 1. The dashed horizontal line shows the lowest rovibrational
ground state of this system.

28
3.2.2

Accuracy of the rovibrational energies of H+
3

Next we compare our theoretical rovibrational energy levels Ei of H+
3 with experimental
energy levels [1]. Adopting the notation used in Ref. [1], we label the rovibrational states i byquantum numbers (N + , G){v1 , v2l2 }(l|u). v1 is the symmetric-stretch vibrational quantum number,v2 denotes the quantum number of the asymmetric-stretch mode, l2 describes the quantum numberof the vibrational angular momentum, and G ≡ |K + − l2 |. The fact that G instead of K + is a goodquantum number implies that the Coriolis interaction couples rotational and vibrational angularmomenta and makes levels with the same G nearly degenerate. However, for levels with l2 6= 0 and(N + − |l2 |) ≥ G ≥ 1, the degeneracy breaks, and we utilize u (or l) to denote the upper (or lower)
energy level; these levels with a u or an l cannot be described by rigid rotator approximations.
Table 3.1 compares the rovibrational energy levels, calculated for N + ≤ 3 states of {0, 00 }and {0, 11 } bands, with the experimental results of Lindsay and McCall [1]. The agreement isgood, with a rms difference of 0.281 cm−1 for the levels shown in the table. Higher-rovibrationalenergy-level calculations also exhibit good agreement. For our calculated energy levels up to around
9000 cm−1 with N + ≤ 4, the rms difference between our calculation and the experimental resultsof Ref. [1] is 0.657 cm−1 .

3.2.3

Rovibrational-frame transformation

Next we describe in detail how to construct the rovibrational-frame transformation using theionic rovibrational eigenstates. In the laboratory frame, the H3+ + e− system is described by theelectron orbital angular momentum l and its projection λ onto the laboratory z-axis, and by N + ,m+ , gI , and parity of the ion core. Hence, we construct the wave function of the H3+ + e− system
as a sum of products of the ionic rovibrational wave function and electronic wave function, of the
form
+

+

v gI Π
ΨN+ m+ (R, Ω) Ylλ (θe , ϕe ) ,

(3.9)

29

Table 3.1: Comparison of several calculated rovibrational energy levels of H+
3 with experimental
results [1]. Only states with N + ≤ 3 for the {0, 00 } and {0, 11 } bands are shown here.
Q.N.1(1, 1){0, 00 }(1, 0){0, 00 }(2, 2){0, 00 }(2, 1){0, 00 }(3, 3){0, 00 }(3, 2){0, 00 }(3, 1){0, 00 }(3, 0){0, 00 }(0, 1){0, 11 }(1, 2){0, 11 }(1, 1){0, 11 }Ecal 2 
cm−1
64.12886.960169.288237.335315.317427.974494.712516.8232521.183
2547.996
2609.302

Eexp 3
cm−1
64.121(00)86.960(00)169.295(04)237.356(05)315.349(04)428.018(07)494.775(07)516.873(07)2521.411(05)2548.164(11)2609.541(05)Q.N.1(2, 3){0, 11 }(2, 2){0, 11 }(2, 1){0, 11 }l(2, 1){0, 11 }u(2, 0){0, 11 }(3, 3){0, 11 }(3, 2){0, 11 }l(3, 2){0, 11 }u(3, 1){0, 11 }l(3, 0){0, 11 }(3, 1){0, 11 }u

Ecal 2 
cm−12614.0342723.708
2755.313
2790.086
2812.504
2876.566
2931.091
2992.151
3002.348
3025.663
3063.181

Eexp 3
cm−12614.279(11)2723.962(06)2755.565(04)2790.344(04)2812.857(05)2876.847(06)2931.366(05)2992.436(05)3002.905(05)3025.951(08)3063.478(05)a The quantum numbers labeling the energy levels, in the notation (N + , G){v , v l2 }(l|u) described
1 2

in the text.
b Theoretically calculated results from this work.
c Experimentally determined energies from Ref. [1].

30
where θe and ϕe are spherical angles of the electron in the laboratory frame. The radial part ofelectronic wave function is not shown here since we apply the MQDT method for that degrees
of freedom. Specifically, we want to construct the laboratory eigenchannel function |ii with a
N,l,m,Π+ ,gI )definite laboratory-frame total angular momentum N : |ii = |N + , ν + i(

. In the position

representation, it takes the form as,
X

+

+

N,m
v gI Π
CN+ ,m+ ;l,λ ΨN + m+ Ylλ ,(3.10)

λN,mwhere CN+ ,m+ ;l,λ denotes the appropriate Clebsch-Gordan coefficient.

In the body-frame, the is specified state by the projection of theelectron angular momentum
on the molecular Z-axis Λ, by the total angular momentum of the system N including the electron
contribution l and the projection m of N on the laboratory z-axis. Applying the transformation
between body-frame and lab-frameYlλ (θe , ϕe ) =

Xh

i∗

l
(α, β, γ) YlΛ θe0 , ϕ0e ,DλΛ(3.11)

Λthe expansion of the product of two Wigner functions,
+

N
l
Dm
+ K + DλΛ =

X

N,m
N,K
N
DmK
CN
+ ,K + ;l.Λ CN + m+ ;l,λ ,(3.12)

N

we have the following equation,
X

h
i∗
X

N,m
N + ,K +  N ∗N+CNYlλ (θe , ϕe ) =(−)l−Λ Cl,−Λ;N,K
DmK YlΛ θe0 , ϕ0e ,+ ,m+ ;l,λ Dm+ K +λΛ(3.13)

with some manipulation of algebra, where θe0 and ϕ0e are the spherical angles of the electron in
the body-frame. Using Eq. (3.13), we derive the transformation between the body-frame andlaboratory-frame states as,
N,l,Π+ ,gI )hα |ii = hR, Ω; Λ N + , ν + (

+

mgI Π
= Ψ̃N
,
N + v+ Λ(3.14)+

mgI Π
where Ψ̃N
is defined as,
N + v+ Λ
+

mgI Π
Ψ̃N
=
N + v+ Λ

X (v+ )
X (v)mgI.
cnν πi (R)ajm2 K + (Rn )Φ̃N
jm2 N + K + Λnν

jm2 K +

(3.15)

31
mgI
The explicit form of Φ̃Nis given in Appendix C. The rovibrational frame transformation
jm2 N + K + Λ

can be accomplished as follows:(N mgI Π+ )KN + v+ ;N +0 v+0 =

XZ
Λ,Λ0

+

+

mΠ gI ∗mΠ gI
dQdΩE Ψ̃N
KΛΛ0 (Q) Ψ̃N,
N + v+ Λ
N +0 v +0 Λ0(3.16)

where ΩE denotes the Euler angles, and Q contains the three vibrational degrees of freedom.

3.3

p-wave energy levels of H3
The body-frame reaction matrix for a p-wave electron is described by the short-range inter-action extracted from an ab initio calculation [80]. In practice, the quantum defects are smootherthan the reaction matrix elements because the latter can have poles. Hence in this work, we extract the body-frame quantum defects µ (Q) from the ab initio energy surface directly. Afterreplacing KΛΛ0 (Q) by µΛΛ0 (Q) in Eq.(3.16), we perform a rovibrational transformation to get thelaboratory-frame quantum-defect matrix. Finally, we get the laboratory K matrix by using the
eigenvalues µe from the laboratory-frame quantum-defect matrix,
K = U tan (πµe ) U T ,(3.17)

where U denotes the unitary transformation that diagonalizes the laboratory-frame quantum-defectmatrix.

3.3.1

Body-frame quantum defects for p-waves

Because of Jahn-Teller effects, the body-frame quantum-defect matrix is generally not diagonal in the electronic projections Λ, Λ0 . Similar to the body-frame K matrix proposed by Staib
and Domcke [81], it has the form




0
0
 µ00 (Q)µ (Q) = .
0µ11 (Q)µ1−1 (Q) 




0
µ−11 (Q) µ−1−1 (Q)(3.18)

32
We express each matrix element by using the vibrational symmetry coordinates Q = (Q1 , Qx , Qy )as
1
Q1 = f √ (∆r1 + ∆r2 + ∆r3 ) ,31
Qx = f √ (2∆r3 − ∆r2 − ∆r1 ) ,3
Qy = f (∆r1 − ∆r2 ) .

(3.19)(3.20)(3.21)(3.22)

where f = 2.639 255 bohr−1 is a constant and ∆ri describe displacements of the nuclei from theequilibrium geometry at which r12 = r23 = r31 = requi = 1.6504 a.u.. For example, ∆r1 = r23 −requi .
(Qx , Qy ) can be alternatively described by another pair of coordinates (ρ, φ) as Qx = ρ cos φ and
Qy = ρ sin φ. Q1 describes the symmetric stretch of the molecule, while (Qx , Qy ) or (ρ, φ) describe
bends and the asymmetric stretch. These coordinates are very useful here for our Taylor expansionof the body-frame quantum defects around an equilibrium position. We use the following forms,µ00 (Q) = µ00 (Q = 0) + a1 Q1 + a2 Q21 + a3 Q31 + a4 ρ2 ,(3.23)µ11 (Q)(3.24)= µ−1−1 (Q)= µ11 (Q = 0) + b1 Q1 + b2 Q21 + b3 Q31 + δρ2 ,andµ1−1 (Q) = µ−11 (Q) = λρ.

(3.25)

The form of our off-diagonal matrix elements µ1−1 (Q) and µ−11 (Q) differs from the usual adopted
form in Ref. [31] by a phase factor exp (±iφ). This different phase convention is due to the factthat the usual adopted form was derived in a body frame that is rotated from our body frame by
φ/2. We develop the detailed proof in Appendix D. The effective quantum numbers are calculatedby diagonalizing the quantum defect matrix. This calculation gives,
νn,π1,2 (Q) = n − [µ11 (Q) ± |µ1−1 (Q)|] ,(3.26)

33
an expression that can be used to fit the effective quantum numbers calculated ab-initio in Ref.
[80].

3.3.2

3p1 energy levels of H3

We calculate the 3p1 energy levels of H3 and compare them with empirical fits from Ref.
[69]. Quantum defect parameters in Eqs.(3.23–3.25) are extracted from the ab-initio calculation
in Ref. [80]. To fit the experiment results, we shift the quantum defects at equilibrium positionsµ00 (Q = 0) = 0.0683 and µ11 (Q = 0) = 0.4069 by a small amount, 0.0043 and 0.0021 correspondingly.
In Ref. [69], Vervloet and Watson studied the H3 emmision lines of (3s, 3p0 , 3d) → 2p0
bands and (3s, 3p0 , 3d) → 3p0 bands. They then fitted the lines with effective Hamiltonians of thefollowing form,
BN (N + 1) + (C − B) K 2−DN N 2 (N + 1)2(3.27)−DN K N (N + 1) K 2−DK K 4 + . . .
where the explicit expressions can be found in Ref. [2, 3, 4, 5] and Ref. [69]. Table 3.2 compares our
MQDT result with the experimental energy levels calculated from the fitted effective Hamiltonians.
The labels N, g, U are fitting parameters, where N is also the total angular momentum of H3 , and gis related to the quantum number G by G = |g|. Evidently our calculations are in good agreementwith the fitted and recalculated experimental results, with differences of around a few cm−1 .

3.4

Higher angular-momentum states
For higher electronic angular-momentum states with l > 1, the orbits are nonpenetrating

and the short-range interaction is negligible. The long-range multipole potential model employing
perturbation theory has successfully described the high orbital angular Rydberg states of H2 [82,83, 84]. In this work, we include the perturbations and interactions between levels of different n34

Table 3.2: A comparison between several of our calculated 3p1 H3 energy levels with empiricallyfitted experimental energy levels [2, 3, 4, 5].
Label1
N, g, U0, 1, 1
1, 0, 1
1, 1, 1
1, 2, 12, 0, −1
2, 1, −1
2, 1, 12, 2, 12, −3, 1
3, 0, 13, 1, −1
3, 2, −1
3, 2, 13, 3, 1

Ecal 2 
cm−112967.812999.1
13052.3
13066.913139.913056.1
13221.1
13234.3
13203.813450.713300.6
13160.5
13485.4
13453.2

Ef it 3cm−1
12966.863
12998.196
13050.966
13068.700
13138.608
13056.588
13219.125
13235.522
13212.055
13446.072
13300.11913165.030
13483.545
13460.934

4
Differences

cm−1
0.9
0.9
1.3
-1.8
1.3
-0.5
1.9
-1.2
-8.3
4.6
0.5
-4.5
1.9
-7.7

a The label denotes the values of N, G, U adopted in Ref. [69] to fit the experimental energy levels.
b Theoretical results calculated in this study.
c Empirical fits for experimental energies determined in Ref. [69].
d Differences theory - experiment.

35(principle quantum number) and l (angular-momentum quantum number) in a systematic fashion
by incorporating the formalism of MQDT [85]. We use this long-range model to calculate theRydberg states of H3 with l ≥ 2.
For a Rydberg electron with high orbital angular momentum (l ≥ 2 in the case of H+
3 ), theeffects of core penetration are negligible. Hence, the interaction between the Rydberg electron
and the ion core can be approximately described by two effects. First, the interaction potential
between the Rydberg electron and the molecular ion is expanded into a multipole series, where thequadrupole moment of the H+
3 core is the leading anisotropic term. Second, the induced dipole
moment of the ion core interacts with the Rydberg electron by a potential characterized by thepolarizability of the H+
3 core. All higher angular momenta and higher-order polarizabilities are
neglected here, as well as the anisotropic portion of the polarizability interaction.
In this approximation, the Hamiltonian is given in atomic units by
1
1
H = − ∇2 − + Vef f + Hcore ,2r(3.28)

where Hcore is the rovibrational energy of the H3+ core. Vef f includes quadrupole and polarizability
interactions:Vef f = Vquad + Vpol = −



Q2
γ
α
P2 cos θe0 − 4 − 4 P2 cos θe0 .
3
r
2r
3r

(3.29)

where Q2 , α, and γ are respectively the quadrupole moment, isotropic polarizability and thecylindrically-symmetric anisotropic polarizability. Other components of the quadrupole moment
tensor vanish for the undistorted equilateral triangle configuration. For the vibrational groundstate, Q2 , α, and γ are taken from table III of Ref. [86, 87]. The polarizability and quadrupole
interactions are much smaller than the Coulomb interaction and hence will be treated in perturbation theory. We also find that the quantum defect is small (of the order of 0.01), and the coupling
between vibrational ground states of H+
3 to excited vibrational states are negligible. Hence, in
the rovibrational transformation, we only include the vibrational ground state. The body-frame
reaction matrix thus can be written asZ0KΛΛ ≈ −π drfnl (r) hYlΛ | Vef f |YlΛ0 i fnl (r),(3.30)

36

Table 3.3: Comparison between several of our calculated 3d energy levels of H3 with experimentallydetermined energy levels[2, 3, 4, 5].
Label1N +, K +, N
2, 1, 02, 2, 01, 0, 13, 0, 1
1, 1, 12, 1, 13, 1, 12, 2, 13, 2, 13, 3, 1
1, 0, 23, 0, 2

Ecal 2 cm−117399.14
17058.41
17284.81
17742.32
17005.99
17403.89
17698.40
17107.24
17540.96
17204.46
17011.36
17643.36

Ef it 3cm−117415.86
17039.61
17296.57
17741.29
16991.72
17412.83
17700.43
17094.12
17557.32
17188.48
17001.08
17655.58

4
Differences−1
cm-15.8918.80-11.76
1.03
14.28
-8.94
-2.02
13.13
-16.35
15.98
10.27-12.21a The label denotes the quantum numbers N + , K + , N .
b Theoretical results calculated in this work.
c Empirical fits to the experimental energies determined in Ref. [69].
d Differences between theoretical and experimental results.

37
where fnl is the regular Coulomb function with l = 1 as the angular momentum quantum number
and n as the principal quantum number. As the quantum defect for d-wave electrons are small, we
can use integers for n in calculating the radial functions. Here, r is the electronic radial coordinate.
Again, we perform a rovibrational transformation (with only the vibrational ground states) toobtain the laboratory-frame K matrix and finally, calculate the energy levels using the standard
determinantal equation of MQDT.
Table 3.3 compares our theoretical calculations with the experimentally-determined 3d energy
levels [69]. The agreement somewhat poorer than the p-wave case.

3.5

Recombination-pumped triatomic hydrogen infrared lasers
The nice comparisons between our calculated low-lying Rydberg energy levels and the exper-

iment results give us confidence of our method, and encourage us to apply it for higher Rydberg
states. In a recent experiment, mid-infrared laser lines observed in hydrogen/rare gas discharges
are assigned to three-body recombination processes involving an electron, a rare gas (He or Ne)
atom, and the triatomic hydrogen ion (H+3 ). Calculations of radiative transitions between higher
Rydberg states of neutral (H3 ) Rydberg states support this interpretation, and link it to recent
results for hydrogenic/rare gas afterglow plasmas. A mechanism for the population inversion is also
proposed in this section.
In this experiment [44], laser emission lines were produced from an ultrahigh finesse opticalsupercavity containing a supersonically expanding plasma in the spectral range of 930 – 4370 cm−1 ,with the majority occurring near 7 µm (1430 cm−1 ). The supersonic slit expansion generates aweakly ionized plasma with neutral gas density expanding as the reciprocal of distance from theslit nozzle, such that the pressure in the region of observation is in the range of 0.1 to 10 Torr(densities of 1015 – 1017 cm−3 ). Fifty-seven laser transitions were observed in both H2 – He andH2 – Ne gas mixtures. Twenty-nine laser transitions in H2 – He were also observed in H2 O – Heor H2 O – Ne discharges, which is not surprising, as these plasmas are known to produce significant
amounts of hydrogen [88].

38
It was found that laser intensity was highly dependent on partial pressure of H2 in the
expanding gas mixture, maximizing when the mixing ratio was < 1%. No lasing was observed
when pure H2 gas was used for the expansion. This most probably reflects the conversion of the+dominant H+
3 ion to the H5 complex at higher hydrogen concentrations, as characterized by Glosik

et al [89].
To test whether contaminants such as O2 or H2 O could be involved in the lasing observed in
the H2 /He plasma, three experiments were performed. First, O2 was carefully added to the gas line
through a needle valve to determine how the presence of O2 affected laser intensity, revealing thatit decreased significantly when O2 was added to the gas mixture. Second, no lasing was observed
when the H2 cylinder was closed, indicating that the lasing species was indeed generated from thecylinder gas. Third, in case the cylinder was itself contaminated with traces of condensible, thegas was flowed through a long (∼8 ft) liquid nitrogen trap. No difference was observed in thelasing intensity with or without the liquid N2 trap, indicating that the lasing species originatesfrom species produced from pure H2 and He or Ne.
We applied the long-range multipole potentialmodel discussed in the section 3.4 to the d-wave
Rydberg states of H3 . While the d-wave Rydberg states of H3 were computed with this long-range
model, the p-wave Rydberg states are computed separately using the method discussed in section
3.3, since they are mainly dominated by the short-range interaction between the outer electron
and the ion core [90, 91]. These p-wave Rydberg states either autoionize or predissociate much
more rapidly than do the d-wave Rydberg states because of their stronger Jahn-Teller couplings,
or through ordinary l-uncoupling, which causes rotational autoionization and does not rely on the
Jahn-Teller effect. The corresponding d-p H3 dipole transitions thereby satisfy one of the golden
rules of astrophysical lasers, viz. the lifetime of the lower level of a lasing transition should bemuch shorter than the lifetime of the upper level independent of the excitation mechanism;[] hence,
a population inversion can be created by the recombination process, and mid-IR lasing can occur
under proper conditions. Comparison of the measured lasing lines with the theoretical nd to np
lines suggests that metastable H3 created in the expanding plasma is indeed a likely candidate for39
the carrier of these lasing transitions. Note, however, that the candidate lasing lines calculated
in this study are primarily of the type 4d–4p, whereas the ternary recombination mechanism isbelieved to predominantly produce much higher Rydberg levels with principal quantum numbersin the range n ≈ 40 ∼ 60 or even higher. This view of the process therefore requires a cascade, i.e.,a decay of these highly excited nf levels to the 4d states, which then have a population inversion
and are able to lase while decaying to the 4p levels of H3 .
Because there are no laboratory measurements of the 4d – 4p transitions that appear tobe the most likely candidates for the lasing lines reported here, the theoretical model has beentested (in section 2.3 and 2.4) against the experimental 3d to 3p transitions of H3 that were
measured by Herzberg, showing good agreement (the rms error for 14 tabulated transitions is
below 13 cm−1 ) [92]. A comparison between the current experiment and our theoretical calculationis shown in Fig. 3.2 and table 3.4. However, there are a few residual experimental lines having nocorresponding theoretical transition. These may correspond to Rydberg levels of the ion core withother {N + , K + }, although we do not have specific candidates in mind at this point.
The conditions obtaining in the pulsed supersonic slit jet plasma used to generate the observed
IR laser action described here are similar to those employed by Glosik et al. in their studies of H+
3
recombination in flowing afterglow plasmas [89, 93]. The total gas pressure drops from the 1 – 2
atmosphere nozzle backing pressure as the reciprocal of distance from the 300 – 400 µm wide exit
slit, yielding neutral gas densities of 1014 – 1017 cm−3 in the regions probed. Typical fractional
ionizations in such plasmas are ∼10−5 . Studies of the H3 O+ ion in similar pulsed slit jet discharges
yielded number densities near 1010 cm−3 and a rotational temperature near 110 K for this ion.
Glosik et al. measured the temperature-dependent recombination rate of H+
3 with electrons inrecombination-limited He/Ar/H2 plasmas at He densities of 0.5 – 6 × 1017 cm−3 and temperaturesof 77 – 300 K. They identified three different behaviors of the measured recombination coefficientwith respect to hydrogen content. For [H2 ] < 1012 cm−3 (or <10−5 of the total gas density), it
12 <increases with increasing H2 density and the spin states of H+
3 are not equilibrated. For 10[H2 ] < 5 × 1013 cm−3 , the rate is independent of H2 density and the ion is in thermal equilibrium.

40

2500

1500

Lines

Experimental

2000

1000

500

0

B-coefficient

0.01
1200

1300

1400

1500

0.1

1

10

Wavenumber
-1

(cm )

Figure 3.2: Comparison of experimental results with calculated nd → n0 p transitions of H3 . The
experimental laser strength is in arbitrary units with linear scale, while the present theoretical
B-coefficients are in the units of 1022 (m/Js2 ) on a logarithmic scale. The calculated line positions
have a precision of about 13 cm−1

41
Table 3.4: Possible assignment of laser lines observed in this work. We calculate the theory lines
by using models described in this chapter. The experiment lines are chosen from the H2 O – Helaser lines which are also observed in H2 O – He or H2 O – Ne experiments.
Initiala

Finalb

N + ,K + ,N ,nd

K + ,N ,n0 p3,1,3,5d
3,1,1,5d

1,2,5pπ−
1,2,5pπ−

1,1,1,4d
1,1,3,4d
1,0,1,7d

1,2,4pπ+
1,2,4pπ+
0,1,6pπ

3,1,3,6d
3,1,2,6d
1,0,2,4d2,1,3,6d
1,1,1,6d
2,2,3,6d
2,1,3,7d2,1,1,7d2,2,2,4d
1,0,3,4d
3,1,4,7d3,1,3,7d
3,1,2,7d2,2,3,4d
1,1,1,4d
1,1,1,4d
1,1,1,4d
1,1,1,4d
1,0,1,4d

1,3,5pσ
1,3,5pσ0,2,4pπ
1,2,5pσ1,1,5pσ2,2,5pσ1,2,6pπ−
1,1,6pπ2,2,4pπ
0,2,4pπ
1,3,6pπ−
1,3,6pπ−
1,3,6pπ−2,2,4pπ
1,1,4pπ
1,2,4pπ−
1,1,4pπ
1,2,4pπ−
0,2,4pπ

1,1,3,4d

1,2,4pπ2,1,2,4d
1,1,2,4d2,2,0,4d
1,0,2,4d2,1,1,4d2,2,1,4d
1,1,1,4d

1,2,4pπ+
1,1,4pπ2,1,4pπ
0,1,4pπ
1,2,4pπ+2,1,4pπ
1,0,4pπλ

Theorylines(cm−1 )1053.411084.27
1218.47
1231.06
1248.761256.36
1266.441293.941302.57
1312.781316.66
1332.071337.161346.831350.49
1350.74
1368.811370.69
1371.75
1377.141387.33
1392.30
1399.421392.30
1399.421411.911409.74
1417.111421.091439.631442.931445.25
1457.61
1469.151471.771475.02

H2 O-He
lines(cm−1 )1054.631082

Differencec(cm−1 )
1.22
-2.27

B/1021(m/Js2 )2.26
4.24

1231.59
1249.371254.13

0.53
0.61
-2.23

0.10
3.71
0.08

1297.831301.311312.26
1317.49
1332.05
1336.66
1348.071349.141349.85
1367.29
1368.51
1371.181375.35
1385.871392.6
1400.55
1401.131403.421407.2

3.89
-1.26
-0.52
0.83
-0.02
-0.50
1.24
-1.35
-0.89
-1.52
-2.18
-0.57
-1.79
-1.46
0.30
1.13
8.83
4.00
-4.71

961.95
0.44
6.28
16.11
3.90
8.22
20.04
2.92
8.40
25.12
2.34
5.78
0.08
9.59
8.12
0.66
8.12
0.66
0.70

1418.84

1.73

23.93

1440.54
1444.221447.851455.451462.081472.731473.27

0.91
1.29
2.60
-2.16
-7.07
0.96
-1.75

19.17
14.62
35.43
26.57
5.42
26.57
19.68

H2 O – He
intensity
130
100
300
2000
700
110
1700
1400
650
800
500
2400
1000
2400
500
1500
2300
2400
10
300
1500
40
80
2400
60
10
2400
2300
2400
2400
500
2400
2400
2400
125
2300

a The initial states, labeled by the quantum number N + , K + , N, nd.
b The final states, labeled by the quantum number K + , N, n0 p . An additional + or − denotes theλpossible Jahn–Teller splitting of two p states: +(−) denotes the higher (lower) energy level.
c The difference between the theory lines and H O – He laser lines.
242
For [H2 ] > 1013 cm−3 , it again increases with increasing hydrogen density due to the conversion+of H+
3 to the H5 complex [89], which recombines very rapidly, but which apparently does notproduce observable laser action in our experiments. Analysis of these data indicated that threebody recombination processes dominated over two-body processes at pressures of a few Torr (1016 –1017 cm−3 ) and temperatures near 300 K, and that this ternary rate maximizes at temperatures of
130 – 170 K. Hence, we can reasonably conclude that the three-body recombination processes also
dominate the ion loss in our recombination-limited supersonic plasmas. Glosik et al. [89, 93] showthat such processes are ca. 100 times faster than that previously described by Bates and Khare[94], and they present a theoretical model for them. In this model, the resonant H+
3 + e. complex[i.e., the highly excited (n ≈40 – 60) H3 Rydberg molecule], which efficiently decays back into theion and electron via rotational autoionization, collides with a He atom during its lifetime, changing
the l-state of the complex to higher, longer-lived values. In particular, the l = 1(p) electronic states
are strongly coupled to the ion core and autoionize rapidly, whereas l ≥ 2(d) states are much more
longlived. These l-changing collisions of the H3 molecule thus produce the population inversion
required for the laser action observed in our experiments.
Glosik et al. have shown that in the low H+
3 + e. collision energy regime, rotational excitationof the vibrational ground state of the ion core with {N + , K + } = {1, 1} produces long-lived Rydberg
states most efficiently. Hence, the most likely product immediately following collision would be ahighly excited (n ≈ 40 – 60) p-wave Rydberg electron attached to a {N + , K + } = {2, 1} ion core.
A subsequent collision of the H3 neutral molecule with a rare gas atom changes the Rydberg
electron angular momentum to a longer-lived, higher value, and/or the collision of the He atomwith the H+
3 core might deexcite the ion core, which would close the rotational autoionization decayroute. Radiative cascade can then populate the 4d levels that are the primary lasing candidates intransitions to the 4p states. Hence, in the theoretical calculation, we include the transitions froma 4d electron attached to {2, 1}, {1, 1}, and {1, 0} cores to the n→p energy levels of H3 . Other
transitions of the type nf→nd and ng→nf of H3 have also been checked, but they do not appear to
be connected with the present experimental results.

43

3.6

Summary
In this chapter, we have calculated the Rydberg energy levels of H3 molecules. Using an

accurate ab-initio quantum-defect surface and ab-initio core energies of H+
3 , our theoretical results
for the p-wave Rydberg states from the present MQDT calculations are in good agreement withexperimental results from J. K. G. Watson [69]. We also study higher-momentum states by a using along-range multipole potential model in conjunction with MQDT, and find encouraging agreementwith experimental results from Ref. [69]. Finally, applying such models for higher Rydberg states
gives supportive explanation of the mid-infrared laser lines in a recent experiment [44].

Chapter 4

Numerical study of three-body recombination

This chapter focus on another type of three-body systems: recombination of three particles
at ultracold temperature. Here, three-body recombination processes are treated numerically fora system of three identical bosons. The two-body model potentials utilized support many boundstates, a major leap in complexity that produces an intricate structure of sharp nonadiabatic
avoided crossings in the three-body hyperradial adiabatic potentials at short distances. This model
thus displays the usual difficulties of more realistic systems. To overcome the numerical challenges
associated with sharp avoided crossings, the slow variable discretization (SVD) approach is adopted
in the region of small hyperradii. At larger hyperradii, where the adiabatic potentials and couplings
are smooth, the traditional adiabatic method suffices. Despite the high degree of complexity,
recombination into deeply bound states behaves regularly due to the dominance of one decay
pathway, resulting from the strong coupling between different recombination channels. Moreover,the usual Wigner threshold law must be modified for excited incident recombination channels.
Note that the material in this chapter has been published in Ref. [95].

4.1

Introduction
Three-body recombination has attracted much theoretical and experimental research interestin recent years. Recombination is the process in which three free particles collide to form a twobody state and a free particle, with the released kinetic energy being distributed between the final
collisional partners. Such reactions are common and important in chemical reactions and in atomic,

45
molecular, and nuclear physics. In ultracold degenerate Fermi gases [96, 97, 98] recombination
has been used as a process to form weakly bound diatomic states, crucial for the experimental
realization of the BEC-BCS crossover physics. In fact, it was shown in Ref. [99] that the use ofrecombination as an efficient way to produce weakly bound diatomic molecules can be extended
to systems other than fermionic gases. For colliding BECs at precisely-tuned relative velocity, theformation of molecules via 3-body recombination can also be used to form molecules efficiently
owing to a double Bose enhancement [100]. Recombination processes normally release a highamount of kinetic energy, producing atomic losses that often limit the lifetimes of Bose-Einsteincondensates (BEC) [101]. Moreover, three-body recombination has been recognized as one of themost important scattering observables in which features related to the universal Efimov physics canbe manifested [102, 35, 103, 104, 32, 105]. Near a two-body Fano-Feshbach resonance, i.e., when thes-wave scattering length a is much larger than the range r0 of the interactions, Efimov states canoccur, causing interference and resonant effects in recombination. The experimental observation ofthese features in recombination has been recently used as evidence of Efimov physics in ultracold
quantum gases [106].
From the theoretical viewpoint, quantitative calculations of recombination for the typicalalkali atoms used in experiments in ultracold gases are limited by the large number of diatomic
states existing in such systems. Most of the available calculations for recombination for realistic
systems have been confined to model systems possessing just a few bound states and/or systems withsmall scattering lengths, and even these are challenging calculations [107, 108, 47, 109]. As a result,recombination calculations relevant for ultracold gases can only be made in the universal regime|a|  r0 , by using simple potential models (with a few-bound states) or else by simply modeling the
decay into all deeply bound molecular states through a single inelastic parameter [104]. However,in ultracold gases experiments the condition of universality is typically not satisfied, making it
desirable to perform more realistic calculations involving more sophisticated two-body models with,
eventually, a larger number of deeply bound states. This paper develops a methodology within thehyperspherical adiabatic representation that permits the treatment of systems with many bound46
states.
The present study still utilizes two-body potentials models that are, however, designed tosupport many bound states, and therefore mimic three-body collisions for more realistic scenarios.
In the hyperspherical representation, the existence of many bound states leads to a complex set of
sharp nonadiabatic avoided crossings in the hyperspherical potential curves at short distances. The
large number of sharp avoided crossings creates numerical difficulties for the traditional adiabatic
representation as formulated with d/dR couplings [110]. To overcome these numerical difficulties,
one solution is to use the slow variable discretiazation (SVD) method proposed by Tolstikhin etal [75]. The SVD method has been successfully applied to three-body bound-state calculations[110, 63] and three-body collisions for the H + Ne2 system [111, 112]. Those calculations, however,
did not require numerical solution of the hyperradial equation [see Eq. (2.17) below] out to large
distances. To study ultracold collision processes such as recombination in the large scattering length
limit, it is crucial to solve the hyperradial equation out to very large distances. Since application of
SVD over the entire space is demanding in terms of memory and cpu-time, it is in fact much more
efficient to separate the domain of hyperradii into two regimes. At short distances, where many
avoided crossings appear, the SVD method is applied, while at large distances, where the adiabatic
potential curves are smooth, the traditional adiabatic method [54] is utilized.
This two-pronged strategy enables efficient calculation of the three-body recombination rate
at low collision energy, with extremely stable results for a two-body potential model supporting
up to about 10 bound states. This numerical capability of calculating recombination with manybound states permits us to study the final state distribution of the recombination rate, K3 . One
unexpectedly simple finding is that the branching ratio of recombination into a particular final (f )channel, defined as(f ←i)

r3(f ←i)K= P 3 (f ←i)
K3

(4.1)

f

(f ←i)

is the same for different initial (i) three-body collision channels. In the above equation, K3

is the partial recombination from the initial three-body channel i to a particular final channel f .

47
The threshold laws for the partial recombination rates have also been considered, i.e., recombination processes occurring from excited three-body continua. These partial rates are observed todeviate from the usual Wigner threshold law. Specifically, the energy dependence of the partial
recombination rates display a much weaker suppression than the usual Wigner analysis [113, 103]for excited continuum channels. These numerical results can be interpreted as the manifestationof a strong coupling between three-body continuum channels. This is further quantified through aperturbation series expansion of the scattering matrix that reveals the three-body recombination
pathways at low collision energies.
This article is organized as follows. Section 4.2 discusses the numerical details in this study.
Section 4.3 shows the numerical results for the three-body recombination rates , and Section 4.4
presents our analysis of the recombination pathways. Section 4.5 summarizes and concludes.

4.2

Adiabatic hyperspherical representation
The system studied here consists of three identical bosons with masses m1 = m2 = m3 = m

with total angular momentum J = 0. The present study uses a two-body potential model in formofv (rij ) = Dsech

2



rij
r0


,

(4.2)

where the coefficient D is negative, and its magnitude is chosen to be large enough to support 8
to 10 two-body bound states (4 to 5 s-wave bound states). The total interaction potential can beapproximated accurately by a pairwise sum of two-body interaction:
V (R, θ, ϕ) = v (r12 ) + v (r23 ) + v (r31 ) .

(4.3)

The first correction to this expression, the Axilrod-Teller term [114], is included latter and foundunimportant. Using the numerical basis discussed in Ref. [47], Eq. (2.12) is solved for the adiabatic
potential and channel functions.
Figures 4.1 and 4.2 show the typical adiabatic potential curves for the parameter m = 7.2963×
103 a.u., D = −5.500 × 10−5 a.u. and r0 = 15 a.u. Figure 4.1 exhibits several sharp nonadiabatic

48
avoided crossings at short distances. Although in our representation sharp crossings are associatedwith vanishingly small transition probabilities, these avoided crossings can cause several numericaldifficulties when solving for the hyperradial motion in the traditional adiabatic approach. As
mentioned above, such difficulties are overcome by implementing the SVD method described inRef. [75]. Figure 4.2 shows, however, that the adiabatic potentials at long distances are smoothand, therefore, are more suitable for traditional approaches. The R-matrix propagation method
was discussed in subsection 4.2.4 and Appendix B. With this method, R-matrix at very large
distance can be obtained. The scattering matrix S is then given by Eq. (2.22).
The three-body recombination rate is therefore given by [47](f ←i)

K3=

192π 2 (2J + 1)
|Sf i |2 ,
4
µ3b k

(4.4)

where Sf i is the appropriate S-matrix elements, J is the total angular momentum of the system,√
and k = 2µ3b E gives the hyperradial wave numbers in the incident channels.

4.3

Three-body recombination rates
The present numerical study focuses on systems of three identical bosons with total angularmomentum J = 0, with parameters adjusted to represent the 4 He system (m1 = m2 = m3 =7.2963 × 103 a.u. and r0 = 15 a.u., the same as in Fig. 4.1). The two-body potential depth D [seeEq. (4.2)] is adjusted to tune the scattering length a and explore both the positive- and negativescattering length cases while supporting 8-10 bound states. The recombination rate near the unitary
regime (k |a| ≈ 1), is explored next for two sets of typical parameters: D = −5.500 × 10−5 a.u. forthe positive-scattering length case, a = 1020.36 a.u.; D = −5.467×10−5 a.u. for negative-scattering
length case, a = −1096.07 a.u. For both cases a is much larger than r0 (|a|/r0 ≈ 70) and therefore
such calculations are solidly within the universal regime [104, 103].
The black dashed lines in Figs. 4.1 and 4.2 denote the recombination channels, i.e., theeν (R) ≡final state channels of the recombination process. The effective hyperradial potentials U

49

-5

4.0x10

f=1
f=2

f=4
f=3

-5

f=5

2.0x10U(R)(a.u.)

f=6

f=7

f=8

f=9

0.0

-5

-2.0x10-5-4.0x10-5-6.0x10

0

20

40

60

80

100

R (a.u.)Figure 4.1: (Color online) Adiabatic potential curves U (R) at short distances R. Red solid lines
represent the three-body continuum channels, i.e., the initial channels for recombination processes,
and black dashed lines represent the final recombination channels.

50

-10

2.0x10-10

1.0x10U(R)(a.u.)

0.0

i=2

i=3

i=4

i=5

i=1
-10

-1.0x10-10-2.0x10

5000

10000

15000

20000

25000

30000

R (a.u.)Figure 4.2: (Color online) Same as Fig. 4.1, but for the adiabatic potential curves U (R) at large
distances R. This figure contrasts with Fig. 4.1 in the characteristically smooth behavior at large
distances.

51
Uν (R) − Qνν / (2µ3b ) for these channels have asymptotic behavior given bylf (lf + 1)
(f )ef (R) R→∞U≈+ E2b ,
22µ3b R

(4.5)

(f )

where E2b is the two-body bound state (dimer) energies, and lf is the corresponding angularmomentum of the third particle relative to the dimer. The subscript f labels these recombination
channels in ascending order, i.e., from high-to-low two-body bound state energies. In Figs. 1 and
2, red solid lines denote the three-body break-up channels (or entrance channels) whose asymptotic
is described byλi (λi + 4) + 15/4
ei (R) R→∞U≈,
2µ3b R2(4.6)

where λi (λi + 4) is the eigenvalue of the grand angular momentum operator Λ2 (here, λi =0, 4, 6, 8..., where λi = 2 is absent for symmetry considerations). The subscript i labels threebody break-up channels in ascending order, i.e., from low-to-high eigenvalues λi .
ei determines the Wigner threshold
As Ref. [113, 103] points out, the asymptotic form of U
laws for recombination, i.e., the low energy behavior of the recombination rate. A simple extensionof the results of Ref. [113, 103] yields the Wigner threshold laws for all three-body channels as(f ←i)K3∝ E λi(4.7)

Our numerical results, however, show that Eq. (4.7) only holds for the lowest entrance channel(i = 1) while it fails to describe the threshold laws for higher incident channels (i > 1). This isapparent from Figs. 4.3 and 4.4, which show our numerical calculations for recombination withpositive and negative values of the scattering lengths, respectively.
In fact, Figs. 4.3 and 4.4 illustrate that for the partial recombination rate from the lowest
three-body incidence channel (λ1 = 0), the threshold behavior does follow the Wigner threshold(f ←1)

law prediction, K3∝ E 0 . However, for higher incident channels λ2 = 4, λ3 = 6, the thresholdenergy exponent is independent of λi and recombination rates are only proportional to E 1 . Therefore the low-energy suppression for higher three-body break-up channels is much weaker than what
Wigner’s threshold law would predict (see Eq. (4.7)). Note that we have used different line-styles

52
to indicate the recombination rate for different incident channels, and use different color and thickness of lines for different final channels. The solid, dashed, dot, dash-dot and short dashed lines
indicate recombination rate to different recombination channels: f = 1, 2, 3, 4, 5. The thick blacklines, thinner red lines, and thinnest blue lines indicate recombination rates from different incident
channels: i = 1, 2, 3.
Another important property that has emerged from our numerical calculations is that thebranching ratio [Eq. (4.1)] for the three-body recombination rates into different final channels are
the same for the few lowest initial channel in the low collision energy limit (see Figs. 4.5 and
4.6). For instance, for the three different initial channels shown in Fig. 4.5, a case with positive
scattering length, the branching ratio into the highest bound state is about 0.35 for each of thethree lowest incident channels throughout the energy range E . 10µK. Similar results are seen forthe branching ratios at negative scattering lengths, as is documented by Fig. 4.6. Note, however,that the branching ratios for positive scattering lengths are not the same for E > 10µK (see Fig.
4.5) while they remain the same for negative scattering length (see Fig. 4.6). This is a resultof destructive interference effects that reduces the recombination probability for the most weakly
bound recombination channel for positive scattering length, and it is related to the universal Efimov
physics [104, 103]. In fact, such interference effect is only significant in the shallowest final channel.
Hence, if the f = 1 channel is excluded from the summation in the denominator on the right handside of Eq. (4.1), the calculated branching ratio between the deep bound states should be the samefor the whole energy range considered, as the inset of Fig. 4.5 shows.
Both the branching ratio properties uncovered in the present numerical exploration and thedeviations from the recombination Wigner threshold laws can be understood using the analyticalmodel developed in the next section. As we will see, these results are driven simply by the strong
long-range coupling between the three-body incident channels; this analysis gives further insight
into the pathways controlling three-body recombination.

53

-23

10

0

E
-24

-25

10

4

E

1

E

-26

10

3

(

Kf i

)

6

(cm /s)

10

-27

10

-28

10

1

10

2

3

10

10

E( K)
(f ←i)

Figure 4.3: (Color online) Partial recombination rate K3
as a function of the collision energy
E for the positive-scattering length case. The solid, dashed, dot, dash-dot and short dashed lines
indicate partial recombination rates to different recombination channels: f = 1, 2, 3, 4, 5, respectively. The thick black lines, thinner red lines and thinnest blue lines indicate recombination rates
from different incident channels: i = 1, 2, 3. The three green dash-dot-dot lines are proportional toE 0 , E 1 and E 4 as indicated in the figure.

54

-21

10

0

E
-22

10

-23

4

E
-24

10

-25

10

1

E

3

(

Kf i

)

6

(cm /s)

10

-26

10

-27

10

-28

10

0

10

1

10

E( K)

2

10

Figure 4.4: (Color online) Same as Fig. 4.3 but for the negative scattering length case.

55

0

10

-1

r

3

(f

i)

10

4x10
-2

10

3x10

2x10

10
-3

10

-1

-1

-1

-1

10

1

10
1

10

2

10

3

2

10

E( K)

3

10

(f ←i)

Figure 4.5: (Color online) Branching ratio of the calculated recombination rates r3
as functions of the collision energy E for the positive-scattering length case. The line-styles solid, dashed,
dot, dash-dot and short-dashed indicate recombination rate to different recombination channels:
f = 1, 2, 3, 4, 5, respectively. The thick black lines, thinner red lines and thinnest blue lines indicate
recombination rates from three different incident channels: i = 1, 2, 3. The inset shows the branching ratio between the deep bound states only and it excludes the shallowest bound state (see thetext for further details).

56

0

-1

10

r

3

(f

i)

10

-2

10

1

10

2

10

E( K)

3

10

Figure 4.6: (Color online) Same as Fig. 4.5 but for the negative scattering length case.

57

4.4

Dominant recombination pathways
The extensive numerical three-body recombination rates presented in the preceding section

are now interpreted in order to extract the important recombination pathways. Once these are
identified, the surprising low-energy threshold behavior and the branching ratio regularities cited
above become clear.
Our model consists of carrying out a perturbation expansion of the S-matrix and then associating each term to a specific pathway. As a first step, Eq. (2.17) is recast in matrix formasTR + W − E1 F = 0,

(4.8)

1 d2
δµν
2µ3b dR2

(4.9)

where
(TR )µν = −
and


1
d
Wνµ = Uν δνµ −2Pνµ+ Qνµ .
2µ3b
dR

(4.10)The off-diagonal terms of W are treated perturbatively, suggesting that the hyperradial Green’s
function matrix should be defined as the solution ofe − E1 G R, R0 = δ R − R0 1
TR + U(4.11)

e is the diagonal submatrix of W , whose matrix elements coincide with Ueν in Eqs. (4.5)
where Uand (4.6). One can, therefore, write the hyperradial Green’s function asG R, R0 = −πif (R< ) h(+) (R> ) ,(4.12)

where f and h(+) are both diagonal matrices. The matrix elements of f are the solutions of Eq.
(4.11) regular at R = 0, and the outgoing Hankel solutions h are given byh(+) (R) = f (R) + ig (R) ,(4.13)

58
where g represent the corresponding irregular solutions. For the three-body break-up channels,
since the centrifugal barriers are dominant, the regular and irregular energy-normalized basis pairfi and gi are well approximated in terms of Bessel functions asfi (R) ≈pµ3b RJλi +2 (kR) ,(4.14)

gi (R) ≈pµ3b RYλi +2 (kR) .

(4.15)

The above hyperradial Green’s function can now be used to expand the S-matrix in a
distorted-wave Born series,
(0)

(1)

(2)

Sf i = Sf i + Sf i + Sf i + ...,

(4.16)

where(0)

Sf i = δf i = 0,

(4.17)(f 6= i). In Eq. (4.16), the first order expansion of the scattering matrix element is simply given
by,
Z ∞

(1)

Sf i = 2πidRff (R) Wf i (R)fi (R) ,(4.18)

0and the low energy behavior of the S-matrix elements can be easily determined by inspection. The√
integrand in Eq. (4.18) is only significant at small values of kR where fi (R) = RJλi +2 (kR) ∝k λi +2 . Therefore,(1)

Sf i ∝ k λi +2 .

(4.19)In terms of the pathways, the first order S-matrix element is the probability amplitude to transit
from the incident channel i and then tunnel through the centrifugal barrier and scatter into recombination channels at short distances (R ∝ r0 ). Therefore, if the recombination process were solely
described by this pathway, the low energy behavior of recombination would be given by(f ←i)

K3=

192π 2 (2J + 1) (1) 2
Sf i ∝ k 2λi = E λi ,µ3b k 4

recovering the usual threshold laws from Wigner’s analysis [113, 103].

(4.20)

59
The first-order result shown in Eq. (4.20) for the low energy behavior of recombination fails,
however, to explain our numerical coupled-channel results [see Figs. 3 and 5] implying that highorder perturbation terms in Eq. (4.16) are crucial in order to determine the actual threshold laws.
Hence we consider the second order partial-wave Born expansion, given by(2)

Sf i = −2π 2 (I1 + I2 ) ,(4.21)

where,

I1 =

X Z ∞

dRff (R) Wf m (R) fm (R)m6=i6=f 0
Z ∞




dR0 h(+)
R0 Wmi R0 fi R0 ,
m
R
X Z ∞
dRff (R) Wf m (R) h(+)
m (R)

×
I2 =

m6=i6=f

Z R
×(4.22)0
dR0 fm R0 Wmi R0 fi R0 .

(4.23)

0

The first integral I1 describes the quantum amplitude for a pathway in which the systemcoming inward in incident channel i to first scatter into an intermediate state m via a long-range
coupling and then scatters to the final channel f at short distances. I2 describes the amplitude fora different second-order pathway for which the system first scatters into an intermediate state m
at short distances and then scatters into the final channel f in a second step. Accordingly in our
analysis, the most important pathway for all incident channels is the one associated with the I1
term in Eq. (4.21), i.e., the pathways incorporated in I2 are much more strongly suppressed in the
low-energy limit.
Interestingly, the second-order correction for the S-matrix element associated with the lowest
three-body incidence channel (i = 1) can only be associated with an intermediate state (m) which isa deeply-bound molecular channel or else an excited three-body continuum channel. In both cases,
our analysis shows that these contributions are unimportant in the low-energy limit. Therefore,
the threshold law for the lowest three-body channel is still given by Eq. (4.20) [with i = 1, λ1 =600]. For recombination events starting from excited three-body channels (i > 1), however, thesituation is different. In this case the dominant pathway is the one that involves the lowest threebody continuum channel as an intermediate channel (m = 1), with a corresponding second ordercorrection:
(2)
Sf i

2

Z ∞

dRff (R) Wf 1 (R) f1 (R)
≈ −2π
Z ∞ 0



(+)
dR0 h1 R0 W1i R0 fi R0 .
×

(4.24)

R

The long-range coupling W1i between the lowest three-body break-up channel and a higher incident
channel is dominated by the P-matrix element between the two channels. For different i > 1, theP-matrix element P1i follows the same asymptotic behavior [54]P1i (R) ∝1,
R2(R → ∞) .

(4.25)

Using the above equation and definition of Wνµ in Eq. (4.10), the integral in the second line of Eq.
(4.24) has the property that
Z ∞
R

(+)

dR0 h1
R0 W1i R0 fi R0 ∝ k.

(4.26)The integral in the first line of Eq. (4.24) is the same as Eq. (4.18). Therefore, the second-orderscattering-matrix element for the i > 1 three-body break-up channels follows(2)

Sif ∝ k λ1 +2 k = k 3 ,(4.27)

which is larger than the first-order S-matrix for channels i > 1 in the small k limit [see Eq. (4.19)].
Therefore, based on the discussions above, the threshold behavior of the partial recombination rate
for any incident channel can be summarized as


 E0,
(i←f )
K3
∝

 E1,

i = 1,(4.28)
i > 1,

which is consistent with our numerical results shown in Figs. 3 and 5.

61
The present analysis, therefore, demonstrates that the important recombination pathway
for excited three-body incidence channels involves an intermediate transition to the lowest threebody incidence channel, controlled by a strong and long-range coupling between continuum channels[Eq. (4.25)]. This recombination pathway also explains why the relative recombination rate to reachthe same final recombination channel from different incident three-body channels is the same. For
ultracold collisions triggered from every excited three-body incidence channel (i > 1), our analysis
shows that the final state contribution for recombination is mainly controlled by the coupling
between the lowest three-body break-up channel at short distances. Therefore, the corresponding
relative final state contribution are independent of the initial excited three-body channel.

4.5

Summary
The methodology elaborated in this chapter is capable of calculating recombination rateand, similarly, any other three-body scattering observable for systems that possess many two-bodybound states. Our numerical study was performed for systems with up to 10 bound states, but it
can be extended to larger problems with a good level of accuracy. Although our calculations forlarger systems might be limited by memory usage and CPU time, our approach still allows for the
analysis of increasingly more complex systems.
A key outcome is an understanding of the modified threshold laws for partial channel contributions to three-body recombination of three identical bosons with angular momentum J = 0.
Our analysis for the important recombination pathways reveals that the threshold behavior of
the recombination rate for excited three-body incidence channels is significantly less suppressed at
low energy than a simple generalization of the Wigner’s threshold laws predicts. In addition, thebranching ratio of recombination into any given final state f is found to be the same for different
incident channels.

Chapter 5

Origin of the Three-body Parameter Universality in Efimov Physics

Three-body recombination has been recognized as one of the most important scattering observables in which features related to universal Efimov physics can be manifested. This chapter
focuses on studies of Efimov physics for a ≤ 0. The field-free recombination rate at the negativesite has an a4 scaling plus some resonance features. These resonances, as was pointed out initiallyby [35], are closely related with the Efimov states. In principle, there are an infinite number ofthese resonances corresponding to the infinite expected number of Efimov states. In addition, theratio between the positions of nearby resonances are universal, i.e., the ratio does not depend on
the details of the short-range two- and three-body interactions. This fact implies that all resonancepositions can be determined solely by a single parameter, known in the literature as the “three-bodyparameter”. While previous fundamental assumptions in the theory of the Efimov effect suggestthat the three-body parameter should not be universal, recent experiments support the opposite
conclusion. The theoretical investigation in this chapter resolves this apparent contradiction by
elucidating the unanticipated implications of two-body interactions. Our study shows that thethree-body parameter universality emerges because a universal effective barrier in the three-bodypotentials prevents the three particles from simultaneously getting close to each other. Our results
also show the limitations on this universality; it is more likely to occur for neutral atoms and lesslikely to extend to light nuclei.
Note that the material in this chapter has recently been published [115], and some of the text
is adapted from that publication almost verbatim.

63

5.1

Introduction
Efimov physics was discovered by V. Efimov in the 1970’s [26]. He predicted that when the

two-body s-wave scattering length a → ∞, there would be an infinite number of three-body boundstates for a system consisting of three identical bosons. The energies of these three-body statesobey the simple relation:
k+1
k
E3b
/E3b
= exp (−2π/s0 ),(5.1)

where the integer k denotes the energy levels in ascending order. The geometric scaling factor
exp (π/s0 ) ≈ 22.7 is universal, i.e., it does not depend on the details of the two-body interaction
nor additional three-body forces as long as both two- and three-body interactions are shortrange.
Equation (5.1) implies that we only need to know one additional “three-body parameter” todetermine the absolute values of all the weakly bound Efimov spectra. One of the most fundamental
assumptions underlying our theoretical understanding of this peculiar effect is that this three-bodyparameter should encapsulate all the details of the interactions at short distances [33]. For thisreason, the three-body parameter has been viewed as nonuniversal since its value for any specificsystem would depend on the precise details of the underlying two- and three-body interactions[116, 117, 118]. Indeed, in nuclear physics, this picture seems to be consistent, i.e., three-bodyweakly bound state properties seem to be sensitive to the nature of the two- and three-body shortrange interactions [117]. However, a recent experiment in ultracold atomic gases to study Efimov
effects gives the opposite conclusion: the three-body parameter is universal for several different
alkali atomic species.
Although the universality of the geometric scaling factor implies that the Efimov physics
can be observed in any system consisting of three identical bosons, ultracold atom gases provide
a perfect playground because of the extraordinary degree of control for such systems. One of the
key experimental tools is to use a Fano-Feshbach resonance to magnetically tune the scattering
length [119]. Tuning the scattering length a from infinity to finite and negative, but still muchlarger than the interaction range r0 , the Efimov energies cross the three-body threshold one by

64
th Efimov state) and cause resonances in
one at some scattering length a−k (corresponding to the k−
the three-body recombination. The resonance positions also obey a simple relation: a−k+1 /ak =−exp (π/s0 ). People usually use a−3b ≡ a0 corresponding to the Efimov ground state as the three0 or the corresponding wave
body parameter. Equivalently, the Efimov ground state energy E3b
0 /~2 1/2 , where m is the atomic mass, can also be used as the three-bodynumber κ∗ = ma E3b
aparameter.
Recently, Berninger et al. measured a−
3b near four different Fano-Feshbach resonances withultracold atomic gases of 133 Cs. Even though the short-range physics was expected to vary fromone resonance to another, Efimov resonances were found for values of the magnetic field at which
1a=a−3b =−9.1(2)rvdW , where rvdW is the van der Waals length [120].

Therefore, in each of thesecases, the three-body parameter was approximately the same, thus challenging a fundamental
assumption of the universal theory. Even more striking is the observation in Ref. [6] that the Efimov
resonance positions obtained for 39 K [7],2 7 Li [8, 9, 10], 6 Li [11, 12, 13, 14], and 85 Rb [15] are also
133 Cs. These observations provide strong evidence
consistent with values of a−3b /rvdW found forthat the three-body parameter has a universal character for spherically symmetric neutral atoms,
therefore suggesting that something else beyond the universal theory needs to be understood.

5.2

Van der Waals interaction and classical suppression
Our theoretical model for two-body interactions mimics the tunability of the interatomic

interactions via Fano-Feshbach resonances by directly altering the strength of the interparticle interactions and, consequently, leads to the desired changes in a. A Fano-Feshbach resonance isfundamentally a multichannel problem. However, for a broad resonance that is open-channel dominated, a single channel potential with van der Waals tail is a very good approximation [119]. This
work considers various model interactions to test the universality of our three-body calculations.
1

The van der Waals length is defined as rvdW ≡ (2µ2b C6 /~2 )1/4 /2 where C6 is the van der Waals coefficient andthe two-body reduced mass µ2b . Note also that in Ref [6] the results were quoted in terms of the mean scattering
length ā ≈ 0.9556rvdW as defined in Ref. [120].
2Here, we are speculating that the feature observed in this experiment at a = −11.02 a.u. might in fact bea three-body resonance, instead of a four-body resonance. The possibility of such a reassignment is by no means
proven, of course, and can only be answered through additional experimental studies.

65
The models we have used for the two-body interactions are:vsch (r) = −Dsech2 (r/r0 ) ,

(5.2)


vλa (r) = − Cr66 1 − λ6 /r6 ,

(5.3)


vλb (r) = − Cr66 exp −λ6 /r6 ,(5.4)

hs (r) = B Θ (r − r) − C6 Θ (r − r ) .
vvdW
hs
hs
hs
r6

(5.5)

The potential model in Eq. (5.2) is the modified Pöschl-Teller potential, where D determines thepotential depth. Equation (5.3) is the usual Lennard-Jones potential, and Eq. (5.4) is a dispersion
potential with a soft wall at short range. In Eq. (5.5), Θ(x) the step-function [Θ(x) = 0 for x < 0
and 1 elsewhere. In practice, however, we have used a smooth representation of the step-function
Θ in order to simplify our numerical calculations.] The potential in Eq. (5.5), therefore, consists
6 ) and a long-range dispersion −1/r 6 potentialof a hard-sphere potential for r < rhs (Bhs  C6 /rhs

for r > rhs . In the present study, the parameters D and λ in Eqs. (5.2)–(5.4) are adjusted to givethe desired a and number of bound states. For convenience, we denote the values of D and λ at
which there exist zero-energy bound states (|a| → ∞) as Dn∗ and λ∗n , where n corresponds to thenumber of s-wave bound states. For the potential model in Eq. (5.5), however, we adjusted rhs
to produce the changes in a, but we only performed three-body calculations near the first pole at
rhs ≈ 0.8828rvdW . While vsch (r) is usually used to mimic the interaction potential between nuclei,
potentials with van der Waals tails are very good models for interaction between neutral atoms.

5.2.1

Van der Waals interaction

One of the important properties of atomic gases is that the long-range two-body interaction
between two neutral atoms can be very well described by a van der Waals interaction: V2b (r) →−C6 /r6 , where r is the interatomic distance. The asymptotic parameter C6 is known quite well for1/4
most atomic pairs of interest, and determines the van der Waals length rvdW = 2µ2b C6 /~2/2,2
where µ2b is the two-body reduced mass. A corresponding energy scale EvdW = ~2 / 2µ2b rvdW,
called van der Waals energy, can also be defined. In general, when r  rvdW , the interaction is

66
described by the asymptotic van der waals interaction; when r  rvdW , V2b (r) is dominated bycomplicated short-range physics.
Flambaum et al. have shown that, for a short-range potential well with a −Cm /rm tail (wherem are integers), the low-energy two-body scattering properties have some “universal” features,providing that the short-range potential is “deep” enough to enable one to use semiclassical (or
WKB) approximation within the potential well [121, 120]. In this case, the s-wave scattering lengtha and effective range reff are related byreff = Fm −

Gm Hm
+ 2 ,
a
a

where Fm , Gm , and Hm only depend on the characteristic length scale rm = 2µCm /~2(5.6)
1/(m−2)/2(r6 ≡ rvdW for van der Waals potential) and ν = 1/ (m − 2):4 π
Γ (ν) Γ (4ν)
rm ,
ν 2ν
3 sin νπ[Γ (2ν)]2 Γ (3ν)

(5.7)

16 π
Γ (1 − 2ν) Γ (4ν) 2ν 4ν
r ,3 sin νπνΓ (ν) Γ (2ν) Γ (3ν) m(5.8)

Γ (1 − 3ν) Γ (1 − ν) Γ (4ν) 3
16 π
ν 6ν
rm .
3 sin νπ
ν 2 [Γ (ν)]2 [Γ (2ν)]2

(5.9)

Fm =
Gm =
Hm =

Equation (5.6) implies that the relation between effective range and scattering length only depends
on the long-range behavior of the potential and does not depend on the complicated short-range
physics. In another word, the relation is “universal”. For a van der Waals potential, the effective
range can be expressed as,
reff ≈ rvdW

r
 
rvdW
vdW 2
2.78946 − 5.33333+ 5.09856.
aa(5.10)

In particular, at the pole a → ∞, the effective range approximately equals to 2.78946 rvdW .

5.2.2

Classical (WKB) Suppression

The universal relation between effective range and scattering length for potentials with vander Waals tails relies on the fact that the two-body zero-energy wave function at short-range can bewell described by the WKB approximation. When the short-range potential well is deep, the WKB

67
approximation becomes very good, and the zero-energy wave function amplitude at short distances
becomes very small. This WKB suppression could derive from the usual classical suppression (orWKB suppression) of the probability for two particles to exist between r and r+dr in regions of highlocal velocity ~kL (r), which is proportional to [~kL (r)/mdr]−1 (m being the particle mass), thetime spent in that interval dr. It is possible that there could be an additional suppression as well,through quantum reflection from a potential cliff [122]. Systems supporting many bound states,
such as the neutral atoms used in ultracold experiments with their strong van der Waals attraction,clearly exhibit this suppression. In general, a finite-range two-body potential that supports manybound states decreases steeply with decreasing interparticle distance r, starting when r/rvdW . 1,
at which point the potential cliff plays a role analogous to a repulsive potential for low-energy
scattering. As will be seen latter, the origin of the universality of the three-body parameter is also
related to the suppression of the probability of finding two particles at distances r < rvdW . This
suppression at the two-body level is studied here with different numerical models.
The first two-body model utilized in this chapter is the Lennard-Jones potential vλa . Figure
5.1 shows three Lennard-Jones potentials with λ = λ∗2 , λ∗3 , λ∗4 so that the 2nd, 3rd and 4th s-wave
bound states are exactly at the threshold, and a → ∞ for all the three cases. It is clear that thesepotentials agree with −C6 /r6 for r  rvdW and are very different at short distances r  rvdW . The
inset shows the corresponding zero-energy wave functions. The wave functions again agree with eachother at large distances and differ at short distances. However, the differences at short distances
are very small, mainly because the wave function amplitude itself is very small. The amplitudeof the zero-energy wave function at short-range becomes small because of classical suppression.
Next, a detailed quantitative study of the suppression is carried out with all the model potentials
in Eqs. (5.2)–(5.5).
The potential models in Eqs. (5.3)–(5.5) have in common the same large r behavior givenby the van der Waals interaction with the characteristic length scale rvdW . Therefore, in order tocompare the results from these models to those for vsch , we define an equivalent rvdW for vsch throughthe relationship between rvdW and the effective range reff for vsch , namely, rvdW ≈ reff /2.78947

68

0
-200

1.0

0.8

-600

*
3
*

=

4

2

C r

|

-800

2

=

0.6

|

V (r)

-400

*

=

-

0.4

6

/

6

-1000
0.2

-1200
0.0
0.2

-1400

0.4

0.6

0.8

1.0

1.2

1.4

1.6

1.8

2.0

r
1

2

3

4

5

r
Figure 5.1: Lennard-Jones potential with the 2nd, 3rd and 4th s-wave bound state exactly at
the threshold (denoted by λ = λ∗2 , λ∗3 , λ∗4 ). The insert shows the corresponding zero-energy wave
function.

69[121, 120], which is valid for values of |a|  r0 . In fact, for vsch we have found that reff (|a| = ∞)strongly depends on the potential depth Dn∗ (shown in Fig. 5.2), while for van der Waals type ofinteractions, such as those in Eqs. (5.3)-(5.5), reff is fixed by reff ≈ 2.78947(2µ2b C6 /~2 )1/4 /2.
Next we explore the origin of this suppression at the two-body level with these model potentials. To gain some insight into the likelihood of finding particles inside the potential well, we startby defining the following quantities,ξpin (k)

1
r0

Z r0

|ψk (r)|2 dr,
Z r
1
out
ξp (k) = lim|ψk (r)|2 dr,r→∞ r − r0 r
0=

(5.11)

0(5.12)

where ψk (r) is the two-body scattering wave function at energy E2b = k 2 /2µ2b , defined such thatr→∞ sin(kr + δ)ψk (r) =sin δ,(5.13)

with δ ≡ δ(k) being the s-wave scattering phase shift. This definition for ψk (r), therefore, leads
to a zero-energy (k → 0) wave function of the form: ψ0 (r) = 1 − r/a. [Note that in the above
equations r0 is the characteristic range of the two-body interaction. For the potential model in
Eq. (5.2), r0 is just the quantity in the argument of the sech function, while for the potential
models in Eqs. (5.3)–(5.5) it is defined to be r0 = rvdW ].
The parameters ξpin and ξpout can be associated with the “average” amplitude of the wave
function inside and outside the potential well, respectively. The amplitude inside the well relative
to the amplitude outside the well is therefore defined as:
ξprel (k) =ξpin (k)= 2 ξpin (k) sin2 δ.
ξpout (k)(5.14)

This relative amplitude vanishes in the limit k → 0 (sin δ ≈ −ka) as a result of our choice forthe asymptotic solution in Eq. (5.13), except at |a| = ∞, when δ is an odd multiple of π/2. The
quantity ξpin (k), however, remains finite in the k → 0 limit. Accordingly, this analysis suggeststhat ξpin (k → 0) is the most relevant parameter for our analysis of the suppression inside the well.
Rigorously speaking, ξpin is not a probability, but it does measure the likelihood of finding two

70

10
9
8
7

r

eff

6
5
4
3
2
1

0

500

1000

1500

2000

2500

3000

3500

4000

Dn *
Figure 5.2: Effective ranges at the pole as a function of the potential depths Dn∗ for vsch .

71particles within r < r0 , where the short-range interactions are experienced. Figure 5.3 presents a
typical result for |a|/r0 and ξpin (k → 0) for the two-body potential in Eq. (5.2). It also shows that in
the universal regime near the poles in a, the wave function is suppressed (small ξpin ) and documents
the fact that this suppression becomes more efficient as the potential becomes deeper and more
two-body states are bound. The black filled circles, open circles and open squares in Fig. 5.3 whichshow the values of ξpin at |a| → ∞, a = 5r0 , and a = −5r0 , respectively, illustrate this trend. Note,however, that for values |a| . r0 , the parameter ξpin quickly increases, indicating a higher likelihoodof finding particles inside the potential well. Similar results are also obtained for the potentials vλa
and vλb [Eqs. (5.3) and (5.4), respectively].
The suppression in ξpin (for a fixed value of |a|  r0 ) can be understood to be a result of theusual semiclassical suppression of the wave function. The wave function inside the well is given by
ψWKB (r < r0 ) =C
1/2

kL (r)

Z r
sin


π
kL (r )dr +
,
4
0

0

(5.15)

where C is a normalization constant, and kL2 (r) = 2µ2b [E2b − v(r)] defines the local wave numberkL (r). Therefore, for deep potentials, the suppression of the wave function inside the potential wellis simply related to the factor kL (r)−1/2 that leads to amplitude suppression of the WKB wave
function [Eq. (5.15)] between r and r + dr. Physically, this can be interpreted as the increase of
the local velocity ~kL (r)/m (m being the particle mass) and the corresponding decrease in thetime spent in that interval dr, [mdr/~kL (r)]. Therefore, in the WKB approximation, one expects
ξpin to be simply proportional to 1/kL . This expectation is tested in the Fig. 5.3 plot of the valueof 1/kL (r0 ) (black dash-dotted line), showing that the suppression in ξpin (k → 0) is consistent withthe semiclassical suppression described above. The relation between the universality of the threebody parameter and this classical suppression can be studied in the three-body level through the
adiabatic hyperspherical representation.

72

20

10

10

1

0
0.1

-10
0.01

-20
2

4

6

8

10

12

14

16

Figure 5.3: The red solid curve represents the scattering length, a/r0 , while the green dashed curve
represents the parameter ξpin (k → 0). Both quantities are plotted as functions of the depth D of
the two-body interaction model vsch [Eq. (5.2)], whose values for which |a| = ∞ are indicated in
the figure as Dn∗ , where n is the number of s-wave states. The black circles, open circles, andopen squares are the values of ξpin at |a| → ∞, a = 5r0 , and a = −5r0 , respectively. Their trends
documents the suppression of the ξpin as the number of bound states increases. The results for ξpin
also show higher efficiency of the suppression inside the well for |a|/r0  1. The black dash-dottedcurve shows the suppression factor 1/kL (r0 ), confirming the classical origin of the suppression
mechanism.

73

5.3

Three-body adiabatic hyperspherical representation
Using the three-body potential as a summation of the two-body model potential from the lastsection, Eq. (2.12) is solved for the adiabatic potential and channel functions. Treating problems
with deep two-body interactions — necessary to see strong inside-the-well suppression — requires us
to solve Eq. (2.12) for two-body model interactions that support many bound states, a challenge formost theoretical approaches. Using our recently developed methodology in chapter 4 [95], however,we were able to treat systems with up to 100 two-body rovibrational bound states.
Figure 5.4(a) shows the adiabatic potentials Uν at |a| = ∞ obtained using the two-body
Lennard-Jones potential supporting 25 dimer bound states. At first glance, it is difficult to identify
any universal properties of these potentials. Efimov physics, however, occurs at a very smallenergy scale near the three-body breakup threshold. Indeed, a closer analysis of the energy range
2[Fig. 5.4 (b)] reveals the universal properties of the key potential curve controlling|E| < ~2 /mrvdW

Efimov states and universal scattering properties.
Figure 5.4(b) shows one of our most important pieces of theoretical evidence for universalityof the three-body parameter: the effective adiabatic potentials Wν obtained using vλa for more andmore two-body bound states converge to a single universal curve. [In some cases in Fig. 5.4 (b) wehave manually diabatized Wν near sharp avoided crossings in order to improve the visualization.
The details are discussed in Appendix E.] As one would expect, the usual Efimov behavior for theeffective potentials, Wν =−~2 (s20 + 1/4)/2µ3b R2 with s0 ≈ 1.00624, is recovered for R > 10rvdW . It
is remarkable, however, that Wν also converge to a universal potential for R < 10rvdW and, more
importantly, these effective potentials display a repulsive wall or barrier at R ≈ 2rvdW . This barrier
prevents the close collisions that would probe the short-range three-body physics, including threebody forces known to be important in chemistry, thus making the three-body parameter universal,
as we will confirm below. This is in fact our most striking result: a sharp cliff or attractionin the two-body interactions produces a strongly repulsive universal barrier in theeffective three-body interaction potential.

74

0

-1

0.4
0.2

-2

0
-0.2
-0.4
-0.6

0.5

1

1.5

2

2.5

3

...

-0.8

-3

-1

3.5
0

4
2

4

6

8

10

Figure 5.4: (a) Full energy landscape for the three-body potentials at a = ∞ for our vλa model
potential. (b) Effective diabatic potentials Wν relevant for Efimov physics for vλa with an increasingly large number of bound states (λ∗n is the value of λ that produces a = ∞ and n s-wave boundstates). The Wν converge to a universal potential displaying the repulsive barrier at R ≈ 2rvdWthat prevents particles access to short distances. (c)–(e) demonstrate the suppression of the wave
function inside the potential well through the channel functions Φν (R; θ, ϕ) for R fixed near theminima of the Efimov potentials in (b). (c) shows the mapping of the geometrical configurations
onto the hyperangles θ and ϕ. (d) and (e) show the channel functions, where the “distance” fromthe origin determines |Φν |1/2 , for two distinct cases: in (d) when there is a substantial probabilityof finding two particles inside the potential well (defined by the region containing the gray disks)
and in (e) with a reduced probability — see also our discussion in Fig. 5.5. In (d) and (e), we used
the potentials vsch and vλa , respectively, both with n = 3.

75

Figure 5.5: Density plot of the three-body probability density |Φν (R; θ, ϕ)|2 sin 2θ that determines
the three particle configuration [see Fig. 5.4 (c)] in the θ-ϕ hyperangular plane for a fixed R (sin 2θ
is the volume element). (a)–(d) show the results for an R near the minima of the Efimov potentials
in Fig. 5.4 (b) for the first four scattering length poles of the vλa model, as indicated. (a) showsthat there is a negligible probability of finding the particles at distances smaller than rvdW (outerdashed circle) and, of course, inside the 1/r12 repulsive barrier (inner solid circle). For higher poles,
i.e., as the strength of the hard-core part of vλa potential decreases, the potential becomes deeper
and penetration into the region r < rvdW is now classically allowed. Nevertheless, (b)–(e) showthat inside-the-well suppression still efficiently suppresses the probability of finding particle pairs
at distances r < rvdW , found to be in the range 2–4%.

76
Qualitatively, this universality derives from the reduced probability of finding particles inside
the attractive two-body potential well. This effect can be seen in terms of the channel functions
Φν , in Figs. 5.4 (c)–(d) and the hyperangular probability densities in Fig. 5.5. In the adiabatic
hyperspherical representation, the space forbidden to the particles fills an increasingly larger portionof the hyperangular volume as R decreases. This evolution can be visualized as the dashed lines
in Fig. 5.5 (a)–(d) expanding outward. In the process, the channel function Φν is squeezed into anincreasingly smaller volume, driving its kinetic energy higher and producing the repulsive barrier inthe universal Efimov potential. Moreover, this suppression implies that the details of the interaction
should be largely unimportant. Consequently, different two-body model potentials should give
similar three-body potentials. Figure 5.6 demonstrates this universality by comparing Wν obtained
from different potential models supporting many bound states. Perhaps more importantly, it
compares them with the results obtained from a two-body model that replaces the deep well by ahard wall, essentially eliminating the probability of observing any pair of atoms at short distances.
hs
in Supplementary Eq. (5.5)]. Quantitatively, however, the fact that the barrier occurs[See vvdW

only at R ≈ 2rvdW indicates that universality might not be as robust as in the cases discussed inRefs. [123, 124, 125, 126]. It is thus important to quantify the value of the three-body parameterto assess the size of nonuniversal effects. The next section focuses on the numerical study of thethree-body parameter with different two-body model potentials and quantifies the universality.
In principle, the three-body parameter could be defined in terms of any observable related
to Efimov physics [104]. Two of its possible definitions are [104]: the value of a = a−3b < 0 at which
the first Efimov resonance appears in three-body recombination (see for instance Ref. [35]) and0 |/~2 )1/2 , where E 0 is the energy of the lowest Efimov state at |a| → ∞.
κ∗ = (m|E3b3b5.4

Three-body parameters
Strictly speaking, when there are two or more s-wave bound states in the two-body model,

the Efimov states are not real bound states but quasi-bound states with a finite width. This widthis related to the decay processes to deeper dimer channels (than the Efimov channel) because of

77

0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0

2

4

6

8

10

Figure 5.6: The Efimov potential obtained from the different two-body potential models used
here. The reasonably good agreement between the results obtained using models supporting manyhsbound states (vsch , vλa and vλb ) and vvdW
[obtained by replacing the deep potential well with ahard wall but having only one (zero-energy) bound state] supports our conclusion that the insidethe-well suppression of the wave function is the main physical mechanism behind the universalityof the three-body effective potentials. The differences between these potentials are seen to cause
differences of a few percent in the three-body parameter.

78
the couplings at short distances. It is well known that when there is a quasi-bound state buriedin the continua, the scattering amplitude in the continuous state and the quasi-bound state caninterfere with each other and form a Fano resonance. In our case, the scattering matrix elements
Sii for deeper atom-dimer thresholds (denoted by i’s) was calculated as a function of the energy E.
Near the resonance due to an Efimov state, the quantity |1 − Sii |2 (which is closely related withthe elastic scattering cross section) is fitted with the formula|1 − Sii |2 = sa

(q + ε)2+ sb ,
1 + ε2(5.16)
0 / Γ0 /2 gives the resonance position
where sa , sb , and q are fitting parameters, and ε = E − E3b3b0 and the width Γ0 .
(Efimov energy) E3b
3b

Figure 5.7 shows an example of fitting the resonance energy and the width using Fano lineshapes. The system studied here consists of three identical bosons with a total angular momentumof J = 0 interacting with each other via a Lennard-Jones potential vλa . The λ is chosen to be equalto λ∗2 corresponding to the case where the second s-wave bound state is exactly at the threshold.
There are three deeper dimer channels corresponding to a g-wave dimer, a d-wave dimer, and a
deep s-wave dimer. The calculated quantities |1 − Sii |2 for each atom-dimer channel are shownin the black squares, red triangles, and blue circles for the g-wave, d-wave and s-wave channel,
correspondingly. The solid curves are the results of the fitting formula in Eq. (5.16), and show avery good agreement between the fitting and the numerical results. In this case, the Efimov energy
0 ≈ −0.0498E
0E3b
vdW and Γ3b ≈ 0.0078EvdW are deduced from the fitting.

This method can also be applied to finding higher Efimov states and to finite scattering length
cases. Figure 5.8 shows the Efimov scenario using vλb near the second pole λ = λ∗2 . The red curvewith solid circles corresponds to the Efimov ground state, and the black curve with solid squaresymbols represents the Efimov first-excited state. The blue dashed curve shows the shallowest
s-wave dimer energy. At the negative site, the Efimov energies cross the three-body break-up
threshold at some values of a−k , which can be deduced from extrapolations. For the case studied−− −here, a−0 ≈ −9.35 and a1 ≈ −158.37 are deduced from extrapolations. The ratio a1 /a0 ≈ 16.93

79

4.0
3.5

d

ii

|1-S |

2

3.0
2.5
2.0
1.5

g

1.0
0.5

s

0.0
-0.08

-0.07

-0.06

-0.05

-0.04

-0.03

-0.02

E /E

vdW

Figure 5.7: Fitting the Efimov resonance using a fano lineshape [Eq. (5.16)] for a system using thetwo-body model potential vλa with λ = λ∗2 . The blue circles, red triangles, and black squares are
the numerically calculated |1 − Sii |2 for the three deeper atom-dimer channels (a g-wave dimer, ad-wave dimer, and a deeper s-wave dimer, correspondingly.) The curves are fitting results fromusing Eq. (5.16).

80
deviates from the universal value 22.7 because of the finite range correction. When an Efimov statecrosses the three-body break-up threshold, a resonance shows up in the three-body recombination
rate.
Figure 5.9 shows the three-body recombination rate as a function of scattering length andexhibits a clear resonance structure. A zero-range model gives an analytical formula for the threebody recombination rate near an Efimov resonance at zero temperature (E → 0) [104]:
K3 = AK2sinsinh (2η)
.
s0 ln a/a−
+ sinh2 η
0(5.17)

This formula is used to fit the numerically calculated three-body recombination rate and to deduce
the resonance position a−0 ≈ −9.18. This value is very close to the position of an Efimov statecrossing the threshold given by the extrapolation method mentioned in the last paragraph (which
is about -9.35). The small difference might be due to the finite range correction and the finite
temperature corrections.
Directly calculating K3 and extrapolating the Efimov energy can give approximately thesame answer for a−k . This fact has also been checked numerically for several other two-body model
potentials with 2–4 s-wave bound states. However, calculating the Efimov energy is usually fasterand more stable. Therefore, in the next subsection, the extrapolation method is used to study−a−3b ≡ a0 for all the different two-body potential models.

5.4.1

Universality of the three-body parameter

To study the universality of the three-body parameter, different two-body potential models
are utilized. Our numerical results for κ∗ and a−
3b are summarized in Figs. 5.10(a) and (b), respectively, demonstrating their universality in the limit of many bound states. In fact, the values forhsκ∗ and a−3b in this limit differ by no more than 15% from the vvdW results — κ∗ = 0.226(2)/rvdW

and a−
3b = −9.73(3)rvdW [solid black line in Fig. 5.10(a) and (b)] — indicating, once again, thatthe universality of the three-body parameter is dependent upon the suppression of the probability
density within the two-body potential wells. Given this picture, we attribute the variations of κ∗

81

-5

10

-4

10

|

E|/E

vdW

-3

10

-2

10

-1

10

0

10

-0.6

-0.4

-0.2

0.0

0.2

sign(a)(a/l

vdW

)

0.4

0.6

1/4

Figure 5.8: Efimov scenario calculated for vλb with λ = λ∗2 . Both the ground and first excitedEfimov spectra as a function of sign(a)(a/rvdW )1/4 are shown.

82

1E-20

1E-22

3

6

K (cm /s)

1E-21

1E-23

1E-24
4

8

12

16

20

a r

| |/

vdW

Figure 5.9: The Efimov resonance corresponding to the Efimov ground state. The black squares
are numerically calculated results, and the solid curve is the fitting result from using Eq. (5.17).

83
and a−
3b in Fig. 5.10 to the small, but finite, probability of reaching short distances, thus introducing nonuniversal effects related to the details of the two- and three-body forces. Nevertheless, our
133 Cs [34, 6], 39 K [7],results for a−
3b are consistent with the experimentally measured values for
7 Li [8, 9, 10], 6 Li [11, 12, 13, 14],3

hs
and 85 Rb [15], all of which lie within about 15% of the vvdW

result. Curiously, if one simply averages the experimental values, then the result differs from thehs
vvdW
result by less than 3%.

Previous treatments have failed to predict the universality of the three-body parameter forvarious reasons. In treatments using zero-range interactions, for instance, the three-body parameter enters as a free parameter to cure the Thomas collapse [33], preventing any statement about
its universality. Finite-range models, such as those used in some of our own treatments [118] [corresponding to the results for vsch with n = 2 and 3 in Figs. 5.10 (a) and (b)], have failed for lack of
substantial suppression of the probability density in the two-body wells. This scenario, however,should reflect better the situation for light nuclei having few bound states and shallow attraction.
In contrast to Ref. [118], other models [127, 128, 129, 125, 130, 131] have found better agreementwith experiment. Our analysis of these treatments, however, indicates that the two-body modhs , therefore satisfying the prerequisite for aels used have many of the characteristics of our vvdW

universal three-body parameter.
A recent attempt [16] to explain the universality of the three-body parameter avoided explicit two-body models altogether, using instead an ad hoc hyperradial potential that bore little
resemblance to our numerical potentials in Fig. 5.4. This ad hoc three-body potential displayedstrong attraction at short distances in contrast to our key finding that a cliff of attraction for
two bodies produces a universal repulsive barrier in the three-body system. Consequently, eventhough a universal three-body parameter was found in Ref. [16], the fundamental understanding
provided by the approach is uncertain. We also construct a simple single channel model according
to our numerical findings. The differences between the two models and our single channel model
For the experiments with 6 Li [11, 12, 13, 14], we have determined a−
3b by using the definition of the mean
scattering length from: Wenz, A. N. et al., Phys. Rev. A 80, 040702(R) (2009).
3

84

0
0.2

0.7

-2
0.6

0.1

-4
0.5

-6

0
0

5

10

15

20

25

0.4

-8

0.3

-10

0.2

-12

0.1

-14
1

2

3

4

5

6

7

8

9

10

1

2

3

4

5

6

7

8

9

10

Figure 5.10: Values for the three-body parameter (a) κ∗ and (b) a−
3b as functions of the number nof two-body s-wave bound states for each of the potential model studied here. (c) Experimental
133 Cs [6] (red: ×, +, , and ∗), 39 K [7] (magenta: M), 7 Li [8] (blue: •) and [9, 10]values for a−
3b for(green:  and ◦), 6 Li [11, 12] (cyan: N and O) and [13, 14] (brown: H and ♦), and 85 Rb [15] (black:hs
results. The
). The gray region specifies a band where there is a ±15% deviation from the vvdW
ininset of (a) shows the suppression parameter ξp [Eq. (5.11)] which can be roughly understood as
the degree of sensitivity to nonuniversal corrections. Since ξpin is always finite — even in the large
n limit — nonuniversal effects associated with the details of the short-range interactions can stillplay an important role. One example is the large deviation in κ∗ found for the vsch (n = 6) model,
caused by a weakly bound g-wave state. For n > 10 we expect κ∗ and a−
3b to lie within the rangeof 15% established for n ≤ 10.

85
are presented in detail in the next subsection.

5.4.2

Single channel approximation

The strong multichannel nature of the problem can be illustrated by comparing the results
obtained from a single channel approximation, i.e., Wνν 0 (R) = 0 (ν 6= ν 0 ), with our solutions of the
fully coupled system of equations. Figure 5.11 (a) shows the three-body parameter κ∗ [related tothe energy of the lowest Efimov state through the relation κ∗ = (mE/~2 )1/2 ] obtained for the vλa
model obtained in the single channel approximation (open triangles) as well as our full numericalresults (open circles). The disagreement between these quantities increases with the number ofs-wave bound states n, meaning that the physics controlling the results becomes more and more
multichannel in nature. Nevertheless, we find that by imposing a simple change in the adiabatic
potentials near the barrier — to make the barrier more repulsive [see Fig. 5.11 (b)] — the single
channel approximation for κ∗ [filled circles in Fig. 5.11 (a)] reproduces the full numerical calculations
much better. This agreement indicates that most of the nonadiabaticity of the problem is related tothe exact shape of the barrier and that, to some extent, the effect of the nonadiabatic couplings is
to make the effective potential Wν more repulsive. For these reasons and, of course, the universalityof our full calculations (see for instance Fig. 5.10), we believe that the short-range barrier in thethree-body effective potentials indeed offers a physically valid explanation of the universality ofthe three-body parameter. We emphasize, though, that Fig. 4 only includes the results from ouressentially exact solutions of the full calculations. The single channel results discussed here areintended only as support of our physical interpretation.
It is within this context that we analyzed the model proposed in Ref. [16]. In Ref. [16],the three-body effective potential important for Efimov physics was estimated by considering thedifferent aspects controlling the physics at small and large distances. At distances comparable to
R = rvdW , it was assumed that the effective three-body potential is dominated by the contributions
from equilateral triangle geometries and included only on two-body interactions. Under theseassumptions, r12 = r23 = r31 = r giving R = 31/4 r (note that our definition for R differs from that

86

0.32
0.3
0.28
0.26
0.24
0.22
0.2
0.18
1

2

3

4

5

6

7

8

9

10

1

0.5

0

-0.5

-1
1

2

3

4

5

Figure 5.11: (a) This figure compares the energies (as characterized by the three-body parameter κ∗ )obtained from a single channel approximation with our full calculations. The three-body parameterκ∗ is shown for the vλa model in the single channel approximation (open triangles) as well as for
our full numerical results (open circles). The single channel approximation can be improved by
imposing a simple change in the adiabatic potentials near the barrier, as is shown in (b). Therehs
we smoothly connect the potential for vλa (red solid line) to the barrier obtained for vvdW
(blackasolid line), resulting in the potential labeled by v̄λ (green solid line). This new potential is actuallymore repulsive and has energies [filled circles in (a)] that are much closer to our full numericalcalculations.

87
used in Ref. [16]), the effective potential can be written as666Vm (R) = −C6 /r12
− C6 /r23− C6 /r31= −3C6 /r6 = −3 × 33/2 C6 /R6 .

(5.18)

This potential is expected to be valid for distances R < Ā = 4πΓ(1/4)−2 31/4 33/8 rvdW ≈ 1.9rvdW[120]. With our method, however, we have the means to check the validity of Eq. (5.18) by comparing it with our numerically calculated potentials. Figure 5.12 shows the three-body potentials
obtained using the vλa (λ = λ∗10 ) model supporting a total of 100 two-body bound states. The
potential from Eq. (5.18) is the black solid line passing near the series of avoided crossings andmight be loosely viewed as diabatically connecting the fully numerical potentials. This relation isreasonable given that this sequence of avoided crossings has been shown in Ref. [132] to be related
to the transition of the system from an equilateral triangle geometry to other geometries. This
figure therefore suggests that approximating the short range physics by Eq. (5.18) is not wholly
unphysical, but its validity depends on a strong assumption of diabaticity through a large numberof avoided crossings and is thus probably not an approximation satisfactory for a more quantitative
analysis.
The potential in Eq. (5.18), however, was not actually used in the calculations in Ref. [16].
Instead, it was used to motivate a claim that three-body quantum reflection plays an importantrole, allowing the short range behavior to be replaced by a repulsive potential for distances R <Ā ≈ 1.9rvdW . It is interesting to note that the value R ≈ 1.9rvdW obtained from Ref. [16] for the
position of the hard wall is quite close to the hyperradius where our potentials exhibit the universal
barrier (see Fig. E.1, for instance), indicating that Ā might have some physical meaning. It is worth
mentioning, however, that the barrier we observed in our calculations is model independent, i.e., it
doesn’t rely on the particular model used for the two-body interaction. For distances R > Ā, themodel in Ref. [16] assumed the three-body effective potentials to be given by the universal Efimov
formula,s2 + 1/4
VE (R) = −~2 0.
2µR2(5.19)

88

1
0
-1
-2
-3
-4
-5
-6
0.2

0.4

0.6

0.8

1

Figure 5.12: This figure shows the three-body potentials obtained using the vλa (λ = λ∗10 ) model
supporting a total of 100 bound states. Roughly speaking, the potential of Eq. (5.18) [16] (blacksolid line) can be seen as a diabatic potential since it passes near one of the series of avoided
crossings.

89
It is well known [?], however, that this potential is valid for rvdW  R  |a|, and R = 1.9rvdW iscertainly out of this range. In fact, one can see in Fig. 1 (b) that the use of the Efimov potential[Eq. (5.19)] for R < 10rvdW is a crude approximation. Nevertheless, using this model, Ref. [16]obtains a−
3b ≈ −9.48rvdW , a value consistent with experiments [34, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15]
as well as with our present calculations for a−
3b .
However, extending the model of Ref. [16] to the limit |a| = ∞, we obtained κ∗ ≈ 0.037/rvdW .
hs
Our result, by way of contrast, using vvdW
[Eq. (5.5)] is κ∗ = 0.226(2)/rvdW . We also have tested theeffects of finite a corrections on the model of Ref. [16] by replacing Eq. (5.19) with the three-bodypotential obtained with a zero-range model of the two-body interactions [133]. These corrections
are particularly important near R = |a|. With this modification, the model of Ref. [16] leads toa−3b ≈ −39.96rvdW . For these reasons, we believe that this model’s agreement with our results andexperimental data is fortuitous.
While multichannel couplings are needed to quantitatively describe this system, it is possible
to construct an effective hyperradial potential curve that correctly describes the behavior of threeatoms in the universal regime. Such a potential curve could be useful for simplified future studies.
The approximate form obtained from the present study is:
22µrvdWs20 + 1/4
2.334
1.348
44.52
4.0 × 104
u
W
(R)
≈
−
−
−
−
+
. (5.20)
ν
~2(R/rvdW )2 (R/rvdW )3 (R/rvdW )4 (R/rvdW )5 (R/rvdW )16

Here µ is the three-body reduced mass and rvdW is the two-body van der Waals length. For
comparison the speculative potential curve proposed by Chin [16] is shown, which does not resemble
the present result even qualitatively at small distances, as it exhibits far too little attraction in theregion R = 2 − 5rvdW .

5.5

Summary
In summary, our theoretical examination shows that the three-body parameter controllingmuch of universal Efimov physics can also be a universal parameter under certain circumstances
which should be realized in most ultracold neutral atom experiments. Provided the underlying

90

1

0

-1

-2

0

1

2

3

4

5

6

7

8

Figure 5.13: Comparison between the effective potential obtained from Ref. [16] and the potentialfrom Eq. (5.20).

91
two-body short-range interaction supports a large number of bound states, or it has some otherproperty leads to the suppression of the wave function at short distances, three-body properties
associated with Efimov physics can be expected to be universal. This surprising new scenario could
not have been, and was not, anticipated from the simple model calculations to date. Ironically,
increasing the complexity of the model simplified the outcome by effectively eliminating the impactof the deeply bound two- and three-body states on the low-energy bound and scattering three-bodyobservables. That is, the three-body parameter becomes largely universal.
While these arguments suggest universality also for the three-body parameter in heteronuclear
systems that exhibit Efimov physics with only resonant interspecies two-body interactions, verifying
this prediction is a high priority for future theory and experiment.
Equally important is the exploration of the relationship between a < 0 and a > 0 Efimov
features — currently a subject of a number of controversies [6] — under the new perspective ourpresent work offers.

Chapter 6

Efimov physics on the positive side

In the previous chapter, Efimov resonances of three-body recombination rates K3 for a < 0have been studied. The universal three-body parameter a−3b ≈ −9.73rvdW is found to be consistent
in recent experiments of ultracold atomic gases. Although through a very different mechanism,
Efimov physics also affects the three-body recombination rates on the positive side (a > 0), whichis the focus of this chapter. On the positive side, the three-body recombination rate shows aseries of minima on top of an a4 overall scaling. In principle, there are an infinite number of++these minima a+k that obey a universal relation ak+1 /ak ≈ 22.7 for three identical bosons, aclear Efimov feature predicted by [102, 35]. The universality in the three-body parameter thatwas discussed in the preceding chapter implies that the position of the first minimum a+
0 (andall other a+k ) should also be universal, and this is confirmed by the numerical studies reported
in the present chapter. The mechanism responsible for these minima is also studied. Reference[35] discovered that the mechanism causing these minima is a quantum interference phenomenon— the so called “Stückelberg minima”. This mechanism is examined through an analysis of thepartial recombination rates in all the final channels. In addition, two universal enhancements inthree-body recombination are found around a ≈ 0.956 rvdW . These two enhancements are found tobe associated with a universal d-wave dimer, crossing the three-body threshold, for systems with
two-body van der Waals potentials .
In this chapter, the Lennard-Jones potential vλa in Eq. (5.3) is adopted to study three-bodyrecombination rates on the positive side. The three-body recombination rates are calculated for93systems having 2 or 3 s-wave bound states for a large range of scattering lengths. Figure 6.1 shows
the recombination rates for the two cases (2 or 3 s-wave). The recombination rates for these twocases agree with each other very well. A Stückelberg minimum shows up at about a ≈ 28.09 rvdW
for both cases, suggesting that the three-body parameter is universal.

6.1

Stückelberg minima
The Stückelberg minima are the results of a destructive quantum interference phenomenon

within the adiabatic hyperspherical picture discussed in Ref. [35]. Since the scattering length a islarge and positive, the binding energy of the shallowest s-wave dimer Edimer ≈ −~2 / 2µ2b a2 is verysmall, forming a very shallow atom-dimer threshold in the three-body adiabatic potential curves. An
avoided crossing between the three-body entrance channel and this shallow s-wave recombination
channel occurs generally near R ≈ 3a. Such a Landau-Zener-type avoided crossing usually implies
a large coupling between the two channels and allows interference between two competing pathway,
leading to Stückelberg minima. Since in view of the Hellman-Feynman theorem, the couplingstrength is inversely proportional to the difference Ui (R) − Uf (R) (here i and f denote the entrance
and exit channels respectively) at some fixed R, the couplings are negligible for deeper channels.
Therefore the Stückelberg minima can only be observed in the partial recombination rates into
the shallowest s-wave channel. Comparing with the Efimov resonances on the negative side, theresonance behavior should show up in the partial recombination rates of every channel. Figure
6.2 shows the total and partial recombination rates near the Stückelberg minimum. This figure
demonstrates that only the shallowest s-wave channel has a Stückelberg minimum.

6.2

Three-body recombination resonances associated with d-wave interactions
This section focuses on two enhancement features of the three-body recombination rates atabout a ≈ 0.956 rvdW . They are found to be associated with a universal d-wave dimer that crosses
the threshold at about the same scattering length. In 2000, Bo Gao predicted that for interactions

94

-20

10

s-wave3 s-wave
2

-22

10

-24

3

6

K (cm /s)

10

-26

10

28.09

-28

10

a

4

-30

10

1

10

100

a/r

vdW

Figure 6.1: The three-body recombination rate K3 as a function of scattering length a. The black
curve with square symbols shows the results from a Lennard-Jones potential with two s-wave boundstates; the red curve with circles illustrates the results for three s-wave bound states. A Stückelberg
minimum appears at about a = 28.09rvdW for both cases; the minimum is indicated by a verticaldashed line. Two enhancement features are also shown for both cases near the small scattering
length of a ≈ 0.956 rvdW . The green dashed line is proportional to a4 , the overall scaling of K3 .

95

1E-23

3

6

K (cm /s)

1E-21

1E-25

1E-27
20

40

60

80

100

a/r

vdW

Figure 6.2: The total and partial three-body recombination rate K3 as a function of scattering
length a. The black curve with square symbols shows the total recombination rate, and the red
curve with circle symbols shows the partial recombination rate for the shallowest s-wave dimer
channel. The other curves shows the partial recombination rate for deeper atom-dimer channels.
The Stückelberg minimum only shows up in the shallowest s-wave dimmer channel, but not the
deeper atom-dimer channels.

96with a van der Waals tail, −C6 /r6 , there is always a d-wave dimer (and dimers with higher angularmomentum l = 4j + 2, where j = 1, 2, 3...) that becomes bound at a universal value of the s-wave
scattering length a = a∗ = 4π/[Γ(1/4)]2 rvdW ≈ 0.956rvdW [134]. The single-channel van der Waals
interaction is a good model for studying a broad Fano-Feshbach resonance in an ultracold atomic
gas that is open channel dominated [119]. We note that when a = a∗ , the effective range reff has
the value 2 (3π)−1 [Γ (1/4)]2 rvdW ≈ 2.789rvdW , which equals the effective range when the s-wave
scattering length goes to infinity.
The universal two-body s-wave scattering length a∗l (when a d-wave or i-wave dimer becomesbound) as a function of number of two-body s-wave bound states is studied numerically, using
Lennard-Jones potentials for the two-body interaction. Figure 6.3 shows the scattering length
when a d-wave (i-wave) dimer becomes bound. The more s-wave bound states, the deeper is thepotential; The change in depth implies different short-range physics. The horizontal dashed line
shows the universal prediction from Bo Gao’s work [134]. The numerical results obtained for a∗l,where l denotes different partial wave, agree well with Bo Gao’s prediction of two-body binding
energies [within 1% (6%) in the case of 10 s-wave bound states for a d-wave (i-wave)]. The results
show better and better agreement for a system with more and more s-wave bound states. In
addition, just after the dimer becomes bound (a ≤ a∗l ), the binding energy can be expressed as alinear function of the scattering length, i.e.,3El ≈ ~2 dl (a∗l − a) /(2µ2b rvdW
),

(6.1)

where d2 ≈ 5.6 and d4 ≈ 42 are approximately universal for all d-waves and i-waves respectively.
To identify the mechanism of these three-body recombination enhancements, we focus on theexample of Lennard-Jones potential with two s-wave bound states. In this numerical example, thed-wave dimers become bound at about 0.995rvdW and the i-wave dimers become bound at about
1.206rvdW . Figure 6.4 shows the three-body adiabatic hyperspherical potential curves for the caseof a ≈ 0.977rvdW and the final recombination channels. A study of the partial recombination
rates helps to discover the enhancement mechanism. Figure 6.5 shows the total recombination

97

1.3

d-wavei-wave

1.2

1.1

*

al /rvdW

1.0

0.9

0.8

0.7

0.6
2

3

4

5

Number of

6

7

8

9

10

s-wave bound states

Figure 6.3: The two-body s-wave scattering length a∗l when a d-wave (l = 2) dimer (black curvewith square symbols) or an i-wave (l = 6) dimer becomes bound (red curve with circle symbols),
shown as functions of the number of two-body s-wave bound states.

98

20

0

-20

i1

U (R)/E

vdW

s2

d2

-40
g1

-60

-80
d1

-100
s1

-120
1

2

3

4

5

R/r

vdW

Figure 6.4: Adiabatic potential curves for a ≈ 0.977. The 2 + 1 channels, which corresponding to a
dimer plus an atom at very large distance, are labeled by a combination of a letter and a number.
The letter denotes the angular momentum quantum number of the dimer, and the number labels
the channels for the same dimer angular momentum from low-to-high dimer energies.

99
rate (K3 ) and the partial recombination rates (K3f ) for every recombination channel f as functionsof the scattering length a, near the appearance of a new d-wave dimer threshold. The three-bodyrecombination rates show two enhancements with peaks A (aA ≈ 1.09rvdW ) and B (aB ≈ 0.98rvdW ).
An enhancement around the peak A appears in all the partial three-body recombination rates
right before the d-wave dimer becomes bound. After the d-wave becomes bound and forms a new
threshold, the K3 rate enhancement B is dominated by the partial rate in this new threshold, whichis indicated in the inset of Fig. 6.5. These features suggest that the peak A might be a resonancedue to that a three-body state crosses the three-body threshold, while the peak B seems to relate
to some threshold behavior. The behavior of the effective potential curves near these enhancements
gives qualitative insight into their origin.
Figure. 6.6 (I) shows the adiabatic hyperspherical potential curves at the resonance peak A.
Only the channels relevant to the resonances are shown in this figure. The red dashed-line showsthe effective hyperspherical potential. We have manually diabatized the potential near an avoid
crossing (see Appendix E) that doesn’t play an important role, and show it in the black solid curve.
The black solid curve shows an outer barrier, and a three-body state can be supported for thispotential. To check whether there is a three-body state near the threshold, a calculation of the
WKB phase at zero scattering energy is carried out:
Z Rb
φWKB =Ra~pdR,−2µ3b Wνν (R)(6.2)

where Ra and Rb are the classical turning points shown in Fig. 6.6 (I). The calculated WKB phase
φWKB is about 0.51π, which is a strong evidence of the existence of a three-body state. It is thisthree-body state associated with the d-wave dimer state crossing the three-body break-up thresholdthat causes the three-body recombination resonances. The physical picture is that the three-bodystate at threshold can cause the system to be within the short-range potential valley for a long time,hence making it more likely that the system will decay into deeper channels. Figure. 6.6 (II) showsthe effective potential Wνν (R) as a function of R around the enhancement peak B. Three points(a), (b), and (c) in the inset of Fig. 6.6 are chosen for the study of the Wνν corresponding to the

100

-28

8.0x10

A

Total

6

-28

3

10

-27

10

-29

-28

4.0x10

0.8

B

3

6

K (cm /s)

6.0x10

s2
g1
d1
s1
i1
d2

K (cm /s)

-28

10

0.9

1.0

1.1

1.2

1.3

1.4

a/r

vdW

-28

2.0x10

0.0
0.8

0.9

1.0

1.1

1.2

1.3

1.4

a/r

vdW

Figure 6.5: Total and partial three-body recombination rate as a function of two-body scattering
length a, shown on a linear scale, near a d-wave dimer that is just becoming bound. The inset
shows the same graph with a logarithmic scale for the y axis. The solid vertical line shows thes-wave scattering length when the d-wave dimer becomes bound. The two dashed vertical lines
show the two peaks of recombination rate enhancement, denoted A and B respectively.

101
solid curves, dashed curves, and dash-dotted curves. For each point, the thin curve is the entrance
channels that do not change much between the three points. The thick curve is the most important
recombination channel: the new d-wave channel. It is clear that from point (a) to (c), the barrier
and the threshold become lower and lower. When the maximum of the barrier is about to be lower
than the three-body threshold, three-body recombination starts decreasing. This decrease may be
due to the fact that when the maximum of the barrier is well above the three-body threshold, i.e.,
at point (a), the wave function of the exit channel has a large amplitude outside the barrier, andhas a large overlap with the entrance channel wave function. When the maximum of the barrier
becomes lower than the three-body break-up channel, i.e., at point (c), the wave function will havelarge amplitude inside the three-body potential valley and have much less overlap with the entrance
channel wave function. The WKB phase calculated for these three-points is also not close to π/2
indicating that this enhancement is not due to a three-body state resonance. It is therefore notshown in the partial recombination rate for channels other than the newly formed d-wave.

6.3

Three-body state associated with the d-wave dimmer
While the Efimov states can be viewed as three-body states associated with an s-wave dimernear threshold, the new three-body state discovered in the previous section is a three-body state thatcan be associated with a d-wave dimer. Because of the existence of deeper atom-dimer thresholds,
this three-body state is actually a quasi-bound state. Figure 6.7 shows the three-body quasi-boundstate energy as a function of the scattering length a (red dots), with the error bars showing thewidths of the quasi-bound states. The red line represents a fitting formula:22
E3b / ~2 / 2µ2b rvdW
= d3b (a∗3b − a) /rvdW + e3b (a∗3b − a)2 /rvdW
,

(6.3)

where a∗3b ≈ 1.09rvdW , d3b ≈ 8.21 and e3b ≈ 10.01. For comparison, the figure also shows the d-wave
energy as the black curve with solid square symbols. The three-body energy crosses the thresholdat a ≈ 1.09rvdW , which corresponds to peak A in K3 . Comparing with the d-wave dimer energy(the black curve with square symbols in Fig. 6.7), the three-body energy is shown to be always

102

(I)

4

W (R)/E

vdW

2

Ra

Rb

0

-2

-4

-6

2

4

6

R/r

8

10

vdW

b

c

(II)

( )

( )

3

0.4

K

Peak Ba

W (R)/E

vdW

( )

0.2

0.93 0.94 0.95 0.96 0.97 0.98 0.99 1.00

a

a

( )

0.0

b

( )

-0.2

c

( )

2

4

6

8

10

R/r

vdW

Figure 6.6: (I) Wνν (R) as a function of hyperradius R at the K3 resonance peak A. (II) W (R)νν asa function of hyperradius R around the K3 enhancement peak B. The inset shows the three points
where the potential curves correspond to.

103deeper than the dimer energy. The difference between the dimer and trimer energy also increases
when the scattering length becomes smaller, and the dimer becomes deeper. However, the widthof the three-body quasi-bound state does not seem to change appreciably.
From Gao’s analysis [134], it is known that there is also an i-wave dimer that becomes boundat a nearby scattering length. However, the i-wave dimer does not affect these two peaks. The
i-wave channel has a very sharp avoided crossing with all the other channels. As a consequence,the partial K3 rate is negligible as compared with the other partial rates. One explanation may
be that forming an i-wave (l = 6) dimer results in the exchange of a large amount of angularmomentum between the dimer and the additional atom. Therefore, although higher partial waves
l = 10, 14, 18... might also be formed at nearby scattering lengths, they are not expected to show
any strong features in the three-body recombination rate.
Finally, Fig. 6.8 compares the three-body enhancements for the model potential of LennardJones type, having either 2 or 3 s-wave bound states. The enhancement peaks of the two different
models differ by only 0.01 rvdW . This agrees with a small observed shift of the scattering length
where the d-wave dimer crosses the threshold. In this figure, the red (black) solid vertical line showswhere the d-wave dimer crosses the threshold for LJ with 2 (3) s-waves. The red (black) dashedvertical line indicates where the three-body state associated with the d-wave cross the threshold for
LJ with 2 (3) s-waves. This result suggests that the three-body state associated with the d-wave
dimer is also universal.

6.4

Summary
In summary, we have calculated the three-body recombination rates for Lennard-Jones po-

tential with 2 and 3 s-wave bound states at the positive side. A universal Stückelberg minimum forthe two cases is found at around a = 28.09 rvdW , and two universal enhancement peaks are foundat about a = 0.995 rvdW . In particular, one of the enhancement peaks is related with a universal
three-body state that is associated with a d-wave dimer.

104

0

Energy (vdw Unit)-1-2
dwave_2body
3body

-3

-4
0.80

0.85

0.90

0.95

1.00

1.05

1.10

a (vdw Unit)
Figure 6.7: Energy of the three-body bound state associated with a d-wave dimer as a function of
scattering length a.

105

-28

8.0x10

1.09

1.08
0.983

-28

0.995

6.0x10

LJ with 2 s-wave-28

4.0x10K

3

6(cm /s)

LJ with 3 s-wave-28

2.0x10

0.0
0.8

0.9

1.0

1.1

a/r

1.2

1.3

1.4

vdW

Figure 6.8: The enhancements for the total three-body recombination rate at about a = 0.995rvdW
for Lennard-Jones potential with 2 and 3 s-wave bound states.

Bibliography[1] C. M. Lindsay and B. J. McCall, J. Mol. Spectrosc. 210, 60 (2001).
[2] I. Dabrowski and G. Herzberg, Can. J. Phys. 58, 1238 (1980).
[3] G. Herzberg and J. K. G. Watson, Can. J. Phys. 58, 1250 (1980).
[4] G. Herzberg, H. Lew, J. J. Sloan, and J. K. G. Watson, Can. J. Phys. 59, 428 (1981).
[5] G. Herzberg, J. T. Hougen, and J. K. G. Watson, Can. J. Phys. 60, 1261 (1982).
[6] C. Chin, R. Grimm, P. S. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
[7] M. Zaccanti, B. Deissler, C. D’Errico, M. Fattori, M. Jona-Lasinio, S. Müller, G. Roati,
M. Inguscio, and G. Modugno, Phys. Rev. C 66, 064001 (2002).
[8] S. E. Pollack, D. Dries, and R. G. Hulet, Science 326, 1683 (2009).
[9] N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, Phys. Rev. Lett. 103, 163202(2009).
[10] N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, Phys. Rev. Lett. 105, 103203(2010).
[11] T. B. Ottenstein, T. Lompe, M. Kohnen, A. N. Wenz, and S. Jochim, Phys. Rev. Lett. 101,
203202 (2008).
[12] T. Lompe, T. B. Ottenstein, F. Serwane, K. Viering, A. N. Wenz, G. Zürn, and S. Jochim,
Phys. Rev. Lett. 105, 103201 (2010).
[13] J. H. Huckans, J. R. Williams, E. L. Hazlett, R. W. Stites, and K. M. O’Hara, Phys. Rev.
Lett. 102, 165302 (2009).
[14] J. R. Williams, E. L. Hazlett, J. H. Huckans, R. W. Stites, Y. Zhang, and K. M. O’Hara,
Phys. Rev. Lett. 103, 130404 (2009).
[15] J. R. Wild, P. Makotyn, J. M. Pino, E. A. Cornell, and D. S. Jin, Phys. Rev. Lett. 108,145305 (2012).
[16] C. Chin, arxiv: , 1111.1484 (2011).
[17] I. Newton, Philosophiæ Naturalis Principia Mathematica, 1687.

107[18] E. W. Brown, An Introductory Treatise on the Lunar Theory, Cambridge Universit Press,
1896.
[19] F. Haake, Quantum Signatures of Chaos, Edition: 2, Springer, 2001.
[20] K. Huang, C. N. Yang, and J. M. Luttinger, Phys. Rev. 105, 776 (1957).
[21] M. H. Anderson, J. R. Ensher, M. R. Matthews, C. E. Wieman, and E. A. Cornell, Science
269, 198 (1995).
[22] C. C. Bradley, C. A. Sacket, J. J. Tollett, and R. G. Hulet, Phys. Rev. Lett. 75, 1687 (1995).
[23] K. B. Davis, M. R. Andrews, N. J. van Druten, D. S. Durfee, D. M. Kurn, and W. Ketterle,
Phys. Rev. Lett. 75, 3969 (1995).
[24] L. M. Delves, Nucl. Phys. 9, 391 (1985).
[25] L. M. Delves, Nucl. Phys. 20, 275 (1960).
[26] V. Efimov, Sov. J. Nucl. Phys. 12, 589 (1971).
[27] J. Macek, J. Phys. B 1, 831 (1968).
[28] U. Fano, Phys. Today 29, 32 (1976).
[29] U. Fano, Rep. Prog. Phys. 46, 97 (1983).
[30] C. W. Clark and C. H. Greene, Phys. Rev. A 21, 1786 (1980).
[31] V. Kokoouline and C. H. Greene, Phys. Rev. A 68, 012703 (2003).
[32] N. P. Mehta, S. T. Rittenhouse, J. P. D’Incao, J. von Stecher, and C. H. Greene, Phys. Rev.
Lett. 103, 153201 (2009).
[33] L. H. Thomas, Phys. Rev. 47, 903 (1935).
[34] T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C. Chin, B. Engeser, A. D. Lange,
K. Pilch, A. Jaakkola, H.-C. Nägerl, and R. Grimm, Nature 440, 315 (2006).
[35] B. D. Esry, C. H. Greene, and J. P. Burke, Phys. Rev. Lett. 83, 1751 (1999).
[36] B. D. Esry and C. H. Greene, Nature 440, 289 (2006).
[37] E. Herbst and W. Klemperer, Astrophys. J. 185, 505 (1973).
[38] W. D. Watson, Astrophys. J. 183, L17 (1973).
[39] L. Trafton, D. F. Lester, and K. L. Thompson, Astrophys. J. 343, L73 (1989).
[40] J. E. P. Connerney, R. Baron, T. Satoh, and T. Owen, Science 262, 1035 (1993).
[41] J. J. Thomson, Philos. Mag. 21, 225 (1911).
[42] T. Oka, Phys. Rev. Lett. 45, 531 (1980).

108[43] B. J. McCall, A. J. Huneycutt, R. J. Saykally, T. R. Geballe, N. Djuric, G. H. Dunn, J. Semaniak, O. Novotny, A. Al-Khalili, A. Ehlerding, F. Hellberg, S. Kalhori, A. Neau, R. Thomas,
F. Osterdahl, and M. Larsson, Nature 422, 500 (2003).
[44] A. J. Huneycutt, Studies of Atmospheric and Interstellar Molecules Using Cavity Ringdown
Spectroscopy, PhD thesis, University of California, Berkeley, 2003.
[45] R. J. Saykally, E. A. Michael, J. Wang, and C. H. Greene, J. Chem. Phys. 133, 234302(2010).
[46] J. G. Frey and B. J. Howard, Chem. Phys. 99, 415 (1985).
[47] H. Suno, B. D. Esry, C. H. Greene, and J. P. Burke, Phys. Rev. A 65, 042725 (2002).
[48] R. C. Whitten and F. T. Smith, J. Math. Phys. 9, 1103 (1968).
[49] B. R. Johnson, J. Chem. Phys. 73, 5051 (1980).
[50] B. Lepetit, Z. Peng, , and A. Kuppermann, Chem. Phys. Lett. 166, 572 (1990).
[51] B. K. Kendrick, R. T. Pack, R. B. Walker, and E. F. Hayes, J. Chem. Phys. 110, 6673 (1999).
[52] B. K. Kendrick, R. T. Pack, R. B. Walker, and E. F. Hayes, J. Chem. Phys. 110, 6673 (1999).
[53] M. Aymar, C. H. Greene, and E. Luc-Koenig, Rev. Mod. Phys. 68, 1015 (1996).
[54] J. P. Burke, Theoretical investigation of cold alkali atom collisions, PhD thesis, Universityof Colorado at Boulder, 1999.
[55] E. Nielsen, D. V. Fedorov, A. S. Jensen, and E. Garrido, Phys. Rep. 347, 373 (2001).
[56] R. B. Nelson, AIAA Journal 14, 1201 (1976).
[57] E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. D. Cros, A. Greenbaum, S. Hammerling, A. Mckenney, S. Ostouchov, and D. Sorensen, LAPACK Users Guide, SIAM, 1992.
[58] U. Naumanm and O. Schenk, Combniatorial Scientific Computing, Champman and Hall,
2012.
[59] J. V. Lill, G. A. Parker, and J. C. Light, Chem. Phys. Lett. 89, 483 (1982).
[60] T. N. Rescigno and C. W. McCurdy, Phys. Rev. A 62, 032706 (2000).
[61] D. E. Manolopoulos and R. E. Wyatt, Chem. Phys. Lett. 152, 23 (1998).
[62] K. L. Baluja and P. G. Burke, Comp. Phys. Comm. 27, 299 (1982).
[63] J. Wang and C. H. Greene, Phys. Rev. A 82, 022506 (2010).
[64] H. Helm, Phys. Rev. Lett. 56, 42 (1986).
[65] H. Helm, Phys. Rev. A 38, 3425 (1988).
[66] A. E. Orel, I. F. Schneider, and A. Suzor-Weiner, Philos. Trans. R. Soc. London, Ser. A 385,
2445 (2000).

109[67] V. Kokoouline, C. Greene, and B. Esry, Nature (London) 412, 891 (2001).
[68] C. Jungen and S. T. Pratt, Phys. Rev. Lett. 102, 023201 (2009).
[69] M. Vervloet and J. K. G. Watson, J. Mol. Spectrosc. 217, 255 (2003).
[70] M. J. Seaton, Rep. Prog. Phys. 46, 167 (1983).
[71] C. H. Greene and C. Jungen, Adv. At. Mol. Phys. 21, 51 (1985).
[72] H. Takagi, J. Phys. B 26, 4815 (1993).
[73] O. Motapon, F. O. W. Tamo, X. Urbain, and I. F. Schneider, Phys. Rev. A 77, 052711(2008).
[74] A. Giusti, J. Phys. B 13, 386 (1980).
[75] O. I. Tolstikhin, S. Watanabe, and M. Matsuzawa, J. Phys. B: At. Mol. Opt. Phys. 29, L389
(1996).
[76] S. F. Santos, V. Kokoouline, and C. H. Greene, J. Chem. Phys. 127, 124309 (2007).
[77] W. Cencek, J. Rychlewski, R. Jaquet, and W. Kutzelnigg, J. Chem. Phys. 108, 2831 (1998).
[78] R. Jaquet, W. Cencek, W. Kutzelnigg, and J. Rychlewski, J. Chem. Phys. 108, 2837 (1998).
[79] M. Pavanello, L. Adamowicz, A. Alijah, N. F. Zobov, I. I. Mizus, O. L. Polyansky, J. Tennyson, T. Szidarovszky, A. G. Császár, M. Berg, A. Petrignani, and A. Wolf, Phys. Rev.
Lett. 108, 023002 (2012).
[80] I. Mistrı́k, R. Reichle, U. Müller, H. Helm, M. Jungen, and J. A. Stephens, Phys. Rev. A 61,
033410 (2000).
[81] A. Staib and W. Domcke, Z. Phys. D: At., Mol. Clusters 16, 275 (1990).
[82] E. E. Eyler and F. M. Pipkin, Phys. Rev. A 27, 2462 (1983).
[83] G. Herzberg and C. Jungen, J. Chem. Phys. 77, 5876 (1982).
[84] C. Jungen, I. Dabrowski, G. Herzberg, and D. J. W. Kendall, J. Chem. Phys. 91, 3926
(1989).
[85] N. R. Badnell and M. J. Seaton, J. Phys. B. 32, 3955 (1999).
[86] J. D. Augspurger and C. E. Dykstra, J. Chem. Phys. 88, 3817 (1988).
[87] J. D. Augspurger and C. E. Dykstra, J. Chem. Phys. 91, 1384 (1989).
[88] A. Fridman, Plasma Chemistry, Cambridge University Press, Cambridge, 2008.
[89] J. Glosı́k, O. Novotný, A. Pysanenko, P. Zakouril, R. Plašil, P. Kudrna, and V. Poterya,
Plasma Sources Sci. Technol. 12, S117 (2003).
[90] V. Kokoouline and C. H. Greene, Phys. Rev. A 69, 032711 (2004).

110[91] J. A. Stephens and C. H. Greene, J. Chem. Phys. 102, 1579 (1995).
[92] G. Herzberg, Annu. Rev. Phys. Chem. 38, 27 (1987).
[93] J. Glosı́k, I. Korolov, R. Plasil, O. Novotny, T. Kotrik, P. Hlavenka, J. Varju, I. A. Mikhailov,
V. Kokoouline, and C. H. Greene, J. Phys. B – At. Mol. Opt. 41, 191001 (2008).
[94] D. R. Bates and S. P. Khare, Proc. Phys. Soc. , 231 (1965).
[95] J. Wang, J. P. D’Incao, and C. H. Greene, Phys. Rev. A 84, 052721 (2011).
[96] J. Cubizolles, S. J. J. M. F. K. T. Bourdel, G. V. Shlyapnikov, and C. Salomon, Phys. Rev.
Lett. 91, 240401 (2003).
[97] S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, C. Chin, J. H. Denschlag, and R. Grimm,
Phys. Rev. Lett. 91, 240402 (2003).
[98] M. W. Zwierlein, C. A. Stan, C. H. Schunck, S. M. F. Raupach, Z. H. S. Gupta, and W. Ketterle, Phys. Rev. Lett. 91, 250401 (2003).
[99] J. P. D’Incao and B. D. Esry, Phys. Rev. Lett. 100, 163201 (2008).
[100] B. Borca, J. W. Dunn, V. Kokoouline, and C. H. Greene, Phys. Rev. Lett. 91, 070404 (2003).
[101] T. Weber, J. Herbig, and H.-C. N. M. Mark, Phys. Rev. Lett. 91, 123201 (2003).
[102] E. Nielsen and J. H. Macek, Phys. Rev. Lett. 83, 1566 (1999).
[103] J. P. D’Incao and B. D. Esry, Phys. Rev. Lett. 94, 213201 (2005).
[104] E. Braaten and H.-W. Hammer, Phys. Rep. 428, 259 (2006).
[105] J. P. D’Incao and B. D. Esry, Phys. Rev. Lett. 103, 083202 (2009).
[106] F. Ferlaino and R. Grimm, Physics 3, 9 (2010).
[107] H. Suno and B. D. Esry, Phys. Rev. A 78, 062701 (2008).
[108] H. Suno and B. D. Esry, Phys. Rev. A 80, 062702 (2009).
[109] Y. Wang, J. P. D’Incao, and B. D. Esry, Phys. Rev. A 83, 032703 (2011).
[110] H. Suno, J. Chem. Phys. 134, 064318 (2011).
[111] G. A. Parker, R. B. Walker, B. K. Kendrick, and R. T. Pack, J. Chem. Phys. 117, 6083
(2002).
[112] F. D. Colavecchia, F. Mrugala, G. A. Parker, and R. T. Pack, J. Chem. Phys. 118, 10387(2003).
[113] B. D. Esry, C. H. Greene, and H. Suno, Phys. Rev. A 65, 010705 (2001).
[114] E. B. S. G. C. Maitland, M. Rigby and W. A. Wakeham, Intermolecular Forces, Clarendon
Press, Oxford, 1981.

111[115] J. Wang, J. P. D’Incao, B. D. Esry, and C. H. Greene, Phys. Rev. Lett. 108, 263001 (2012).
[116] P. Soldán, M. T. Cvitas, and J. M. Hutson, Phys. Rev. A 67, 054702 (2003).
[117] E. Epelbaum, A. Nogga, W. Glöckle, H. Kamada, U.-G. Meißner, and H. Witala, Phys. Rev.
C 66, 064001 (2002).
[118] J. P. D’Incao, C. H. Greene, , and B. D. Esry, J. Phys. B 42, 044016 (2009).
[119] C. Chin, R. Grimm, P. S. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
[120] V. V. Flambaum, G. F. Gribakin, and C. Harabati, Phys. Rev. A 59, 1998 (1993).
[121] G. F. Gribakin and V. V. Flambaum, Phys. Rev. A 48, 546 (1993).
[122] R. Côté, H. Friedrich, and J. Trost, Phys. Rev. A 56, 1781 (1997).
[123] Y. Wang, J. P. D’Incao, , and C. H. Greene, Phys. Rev. Lett. 106, 233201 (2011).
[124] D. S. Petrov, Phys. Rev. Lett. 93, 143201 (2004).
[125] P. Massignan and H. T. C. Stoof, Phys. Rev. Lett. 93, 143201 (2004).
[126] Y. Wang, J. P. D’Incao, , and B. D. Esry, Phys. Rev. A 83, 042710 (2011).
[127] P. Naidon and M. Ueda, Comptes Rendus Physique 12, 13 (2011).
[128] P. Naidon, E. Hiyama, and M. Ueda, arxiv , 1109.5807.
[129] S. Nakajima, M. Horikoshi, T. Mukaiyama, P. Naidon, and M. Ueda, Phys. Rev. Lett. 106,
143201 (2011).
[130] M. D. Lee, T. Köhler, and P. S. Julienne, Phys. Rev. A 76, 012720 (2007).
[131] M. Jona-Lasinio and L. Pricoupenko, Phys. Rev. Lett. 104, 023201 (2010).
[132] D. Blume, C. H. Greene, and B. D. Esry, J. Chem. Phys. 113, 2145 (2000).
[133] S. T. Rittenhouse, N. P. Mehta, and C. H. Greene, Phys. Rev. A 82, 022706 (2010).
[134] B. Gao, Phys. Rev. A 62, 050702(R) (2000).
[135] D. A. Varshalovich, A. N. Moskalev, and V. K. Khersonskii, Quantum Theory of AngularMomentum, World Scientific, singapore, 1988.
[136] H. C. Longuet-Higgins, Advances in Spectroscopy, Interscience, New York, 1961.
[137] N. Gemelke, C.-L. Hung, X. Zhang, and C. Chin, arxiv: , 0812.1317 (2008).

Appendix A

Symmetry of R-matrix

The R-matrix R is a real and symmetric matrix. This symmetry is determined by theproperties of the radial Schrödinger equation [Eq. (2.17)], which can be rewritten as


d2
d
2µ3b
F =
(U − E1) − Q F − 2P
F,
2
2
dR
~
dR

(A.1)

where 1 denotes an identity matrix. Inserting the expansion of the wave function [Eq. (2.16)] into
Eq. (2.20) givesdFe =F + P × F.
dR
The definition of the R-matrix can then expressed as,dR
F + P × F = F.
dR
Taking derivatives of both sides of this equation leads to




 2



d
d
d
d
d
d
R
F +P ×F +R
F
+
P
F
+
P
F
=
F.
dR
dR
dR2
dR
dR
dR
Replacing the second derivative of F by Eq. (A.1) gives us





d
d
2µ3b
d
2
R
F +P ×F +R(U − E1) − P F = (R × P + 1)
F,
2
dR
dR
~
dR

(A.2)

(A.3)

(A.4)

(A.5)

d
with the help of the relation dR
P = −P 2 + Q. Finally, after multiplying both sides of Eq. (A.5)
 d−1with dR
F +P ×F
and using the definition of the R-matrix Eq. (A.3), we arrive at




−1
d
2µ3b
d
d
2
R = −R(U − E1) − P R + (R × P + 1)
F
F +P ×F
dR
~2
dR
dR


2µ3b
= −R
(U − E1) − P 2 R + (R × P + 1) (1 − P × R) .
(A.6)
~2

113
Since P is antisymmetric and P 2 is symmetric, this equation implies that the symmetry of thematrix R at one value of R immediately guarantees symmetry at all R values. In most problems,
the matrix R is diagonal in some region of space, often at small distances (R → 0), which is enoughto prove that it will remain real and symmetric at all hyperradii.
The symmetry of R guarantees the symmetry of the reaction matrix K and the unity of
scattering matrix S. In the asymptomatic region at large hyperradii, the P matrix vanishes. The
relation between R and K is given by, 
R f 0 − g 0 K = f − gK

(A.7)




 0
f − g 0 K = R−1 f − gK .

(A.8)

or equivalently by

Transposing on both sides of Eq.(A.7) gives f − KT g = f 0 − KT g 0 R.

(A.9)

Together, Eqs. (A.8) and Eq. (A.9) give the followings



 

f − KT g f 0 − g 0 K = f 0 − KT g 0 f − gK ,(A.10)

which can be transformed into
K

T



dg
df
f −g
dR
dR





df
dg
= g
−
f K.
dR dR

(A.11)

Since f and g are linearly independent solutions of the Schrödinger equation, their Wronskian has
the good property of
dg
df
dg
df
2µ3b
f −g
=f
−
g=
1,
dR
dR
dR dR
~π

(A.12)

whoich guarantees that KT = K. Therefore, we prove that K is real and symmetric, and consequently, S = (1 + iK) (1 − iK)−1 is unitary.

Appendix B

R-matrix propagation method with the traditional adiabatic method.

The model described in the main text uses the traditional adiabatic approach combined with
an R-matrix propagation method for large values of R, where the P and Q matrices are smoothfunctions of R. One advantage of using this representation is that instead of calculating the valuesof the P and Q matrices at every mesh point in R, we can solve the hyperangular part of the
Hamiltonian at relatively fewer grid points. (The number of grid points is generally set by thecharacteristic wavelength associated with the collision energy.) This strategy, therefore, allows
the use of interpolation and/or extrapolation methods to generate the required much denser gridwithout memory storage problems. Appendix B describes the approach in more detail.
In the traditional method, ψν 0 is expanded asψν 0 (R) =X

cjν,ν 0 πj (R) Φν (R; Ω).

(B.1)

jν

A comparison of Eq. (B.1) with Eq. (2.52) shows that the main differences between the traditional
adiabatic method and the SVD method derive from using different three-body numerical basis sets.
(Notice that in Eq. (B.1), the Φν (R; Ω) are channel functions evaluated at R, while in Eq. (2.52),
the Φν (Rj ; Ω) are channel functions evaluated at Rj .) However, one can show that the expansion
coefficients cjν,ν 0 are the same for the two expansions if Φν (R; Ω) is smooth so that the DVR
approximation can be applied,Z
dRπi (R) Φν (R; Ω) πj (R) ≈ Φν (Ri ; Ω) δij .

(B.2)

115
Equation (B.2) implies that the traditional adiabatic method and the SVD method are equivalent within the DVR approximation. Therefore, when the P and Q matrices change rapidly andare hard to evaluate numerically, it is highly advantageous to choose the SVD method; when the
P and Q matrices are smooth, however, the traditional method is simpler and benefits from lower
memory storage requirements.
Next, the details of the traditional approach are elaborated and the R-matrix propagationfrom a point b1 to another point b2 is explained. Insertion of Eq. (B.1) into Eq. (2.19) yields
Fνν 0 (b1 ) =X

cjν,ν 0 πj (b1 ),

(B.3)

j

Fνν 0 (b2 ) =X

cjν,ν 0 πj (b2 ),

j

X

0
Feνν 0 (b1 ) =δνµ Fµν0 (b1 ) + Pνµ (b1 ) Fµν 0 (b1 ) ,
µ

X

0
Feνν 0 (b2 ) =δνµ Fµν0 (b2 ) + Pνµ (b2 ) Fµν 0 (b2 ) .
µAs the next step, rewrite Eq. (2.11) in the basis of Eq. (B.1) in the matrix form of Eq. (2.59),
with the matrix elements

Hiµ,jν

=

1
2µ3b

Z b2

1
−
2µ3b
Liµ,jν

="πi0 (R) πj0 (R) dRδµν +

b1

Z b2

#
2 (R )Pµν
jUν (Rj ) δµν −
δij
2µ3b

πi (R) Pµν (R) πj0 (R) − πi0 (R) Pνµ (R) πj (R) dR,

(B.4)

b1

 b2
1 
πi (R) δµν πj0 (R) + πi (R) Pµν (R) πj (R)
.
2µ3b
b1

(B.5)

Use of these matrix elements and replacing a1(a2) with b1(b2), the same procedure as
Eqs. (2.63-2.73) accomplishes the matrix propagation.

Appendix C

Permutation symmetry of the basis functions

For convenience, we first consider trial basis functions as+IΦtry = eim2 ϕ RN
K + m+ (α, β, γ) ΦgI uj (θ) ,(C.1)

and use them to construct basis functions with proper permutation symmetry Eqs.(3.5–3.6). The
continuity condition for Smith-Whitten hyperspherical coordinates [52, 47] requires that,Φtry (θ, ϕ, α, β, γ) = Φtry (θ, ϕ + 2π, α, β, γ + π) .

(C.2)

This boundary condition leads to the condition that K + /2 + m2 must be integral. Hence, if K + iseven, m2 must be integral; if K + is odd, m2 must be half integral. We remark here that the parity
Π+ is determined by K + only; Π+ = +1, if K + is even, and Π+ = −1, if K + is odd [52, 47].
The permutation symmetries for the basis functions chosen for each degree of freedom areshown in Table 4. Here, ΦIgI is the nuclear-spin basis function defined as in Ref. [31]. gI equals+

zero for the ortho state, and ±1 for the para state. The rotational part RN
K + m+ (α, β, γ) is givenby,+RN
K + m+ (α, β, γ) =ri∗
2N + + 1 h N +
D
(α,
β,
γ)
m+ K +
8π 2

(C.3)

+

N
where Dm
+ K + are the Wigner D functions of the Euler angles. The phase of the Wigner function

is chosen as in Varshalovich et al. [135]. N + is the total angular momentum of the ion, K + is the
projection of N + onto the body frame z-axis, and m+ is the projection onto the laboratory frame
Z-axis. We also use a set of fifth-order basis splines uj (θ) to expand the wave function in θ. Sinceuj (θ) is unaffected by permutations, it is not shown in Table 4.

117
Table C.1: Permutation symmetry for basis functions in the different degrees of freedom.
Permutation
Operation
P12+

eim2 ϕRN
K + m+

ei4π/3 e−im2 ϕ(−)N +K RN−K + m+

+

+

ΦIgI
+

ei4πgI /3 ΦI−gI

N+

N+

+

ei2πgI /3 ΦI−gI
ei2πgI ΦI−gI+

ei2πgI /3 ΦIgI
ei4πgI /3 ΦI−gI

P23
P31

ei2π/3 e−im2 ϕ
ei2π e−im2 ϕ(−) R−K + m++++(−)N +K RN−K + m+P12 P31
P12 P23

ei2π/3 eim2 ϕ
ei4π/3 eim2 ϕ(−)K RN
K + m++RN+
K m+Application of the antisymmetrization operator A in Eq.(3.7) to Eq.(C.1) leads to:AΦtry =

h
i
+
2π
4π
uj (θ) 1 + (−)K ei 3 (m2 +gI ) + ei 3 (m2 +gI )hi+
N + −imϕ N +Ii 2π(m2 +gI )I
3× eimϕ RN
Φ
−
e
(−)
e
R
Φ
K + m+ gI−K + m+ −gI(C.4)

Following the fact that m2 is integral (half integral) if K + is even (odd), and that gI equals 0 or
±1, we see that the factors in the second line of the right hand side of Eq.(C.4) vanish unless,m2 + gI = 3n if K+ iseven,(C.5)

m2 + gI = 3n + 3/2 if K+ isodd,(C.6)

where n is any integer. Under the conditions Eqs.(C.5–C.6), the factors in the first line of the righthand side of Eq.(C.4) vanish if m2 = 0, gI = 0, K + = 0 when N + is even. Therefore, the basis
functions that obey the permutation symmetry are:
+

+

+

m gI
I
ΦN
= uj (θ) eim2 ϕ RN
K + m + Φ gI ,jm2 K +

(C.7)

if m2 = 0, K + = 0, gI = 0, and N + is odd, otherwise,hi+ m+ g1
N + +K + −im2 ϕ N +im2 ϕ N +
I
I
I
√
u
(θ)
e
R
Φ
−
(−)
e
R
Φ
ΦN
=
+
+
+
+
j
+
gI
−gI ,
K m
−K m
jm2 K
2

(C.8)

where m2 , gI and K + obey the relations Eqs. (C.5–C.6).
+

+

mgI
m gI
The explicit form of Φ̃Nin Eq. (3.15) is closely related to ΦN. Insertion of
jm2 N + K + Λ
jm2 K +

Eq.(3.13) into Eq.(3.10) yields
+

+

mgI
N K
Φ̃N
= uj (θ) (−)l−Λ eim2 ϕ ΦIgI Cl,−Λ;N,K
RN
Km ,jm2 N + K + Λ

(C.9)

118
with K = K + + Λ, if m2 = 0, gI = 0, K + = 0 and N + is odd. Otherwise,hiuj (θ)mgIl−Λ
N + +K + −im2 ϕ I
N + −K +im2 ϕ I
N +K+
N
N
√
Φ̃N
=
(−)
e
Φ
C
R
−
(−)
e
Φ
C
R
,
+
+
g
Km
−g
l,−Λ;N,K
I
I
jm2 N K Λl,−Λ;N,K̃ K̃m
2(C.10)
where K̃ = −K + + Λ.

Appendix D

Body-frame quantum defect matrix elements

The form of reaction matrix K describing the Jahn-Teller coupling of the p-wave electron
can be written as





K01
K0−1 
 K00




KΛΛ0 =  K10
.
K
K
11
1−1




K−10 K−11 K−1−1(D.1)

In perturbation theory, the matrix elements of the reaction matrix near the equilibrium position
obey,
ZKΛΛ0 ≈

−π dr0 f¯εlΛ r0 f¯εlΛ0 r0
Z



∗
θe0 , ϕ0e ∆Ve Q; re 0 YlΛ0 θe0 , ϕ0e ,× dωe0 YlΛ

(D.2)

where ∆Ve = Ve (Q; re 0 ) − Ve (Q = 0; re 0 ), and Ve (Q, re0 ) is the interaction between the ion and theRydberg electron excluding the Coulomb potential, Q are the vibrational symmetry coordinates,
and re 0 = {r0 , θe0 , ϕ0e } are the spherical coordinates of the electron in the body-frame. f¯εlΛ (r0 ) is
the regular phase-renormalized Coulomb wave function defined by, f¯εlΛ = fεl cos (πµΛΛ (Q = 0)) −gεl sin (πµΛΛ (Q = 0)), where {fεl , gεl } are the usual Coulomb wave functions with energy ε andangular momentum l. µΛΛ (Q = 0) are the constant zero-order term of the diagonal quantum
defect matrix elements.
To explore the symmetry properties of KΛΛ0 , we write the Taylor expansion in terms of

120
Q = {Q1 , Q+ , Q− }, where Q+ = ρe+iφ and Q− = ρe−iφ ,KΛΛ0 =Xν
KΛΛ
0 Qν +

ν

1 X νµK 0 Qν Qµ + . . .
2 νµ ΛΛ(D.3)

where each summation is over the subscripts 1, + and −. However, the expansion Eq. (D.3) is
not valid in our body frame coordinates because ∆Ve (Q; re 0 ) is not single-valued at Q+ = Q− = 0.
This fact is demonstrated in the following discussion of our body frame coordinates.
The body frame coordinates used here are the same as the coordinates defined in Eq. (7)of Ref. [47], but we denote them as x0 y 0 z 0 instead of xyz to make our notations consistent. The
z 0 -axis is perpendicular to the plane defined by the three nuclei and the x0 -axis is associated withthe smallest moment of inertia. After manipulating Eqs.(2), (3), and Eq.(7) of Ref. [47] with somealgebra, we can write down the cartesian coordinates of the positions of the ith nuclei (x0i , yi0 , zi0 ) as


ϕ

θ π
−
cos
+ ϑi ,
2 4
2


ϕ

2
θ π
yi0 = − R sin−sin+ ϑi ,
3d
2 4
2
x0i =

2
R cos
3d



zi0 = 0.

(D.4)
(D.5)
(D.6)

where ϑ1 = 5π/6, ϑ2 = −π/2 and ϑ3 = π/6. When θ is very small, following Eqs. (2.5) and Eqs.
(3.19) the two sets of coordinates {Q1 , ρ, φ} and {R, θ, ϕ} have the following relationship,

Q1 = 31/4 f (R − R0 ) ,(D.7)

ρ = 31/4 f Rθ/2,(D.8)

φ = ϕ − 2π/3.

(D.9)

where R0 = 31/4 requi . Therefore, when Q+ = Q− = 0, and hence ρ = 0, the positions of the
ith nuclei (x0i , yi0 , zi0 ) can be written as,

x0i (Q1 , Q+ = 0, Q− = 0) =

√ 



2
Q1
φ
π
R0 + 1/4cos+ ϑi +
,
3d
2
3
3 f

(D.10)

121
√ 



2
Q1
φ
π
0
yi (Q1 , Q+ = 0, Q− = 0) =R0 + 1/4,
sin+ ϑi +
3d2
3
3 f(D.11)zi (Q1 , Q+ = 0, Q− = 0) = 0.

(D.12)

These equations show that the positions of the three nuclei can not be expressed by Q1 ,Q+ and Q− when ρ = 0, and hence the expansion Eq. (D.3) is not valid as ∆Ve (Q; re 0 ) is not
single-valued at Q+ = Q− = 0 and therefore is not infinitely differentiable. However, if we define
another set of xe0 , ye0 , ze0 axis by a passive rotation through φ/2 about z 0 -axes, the positions of nuclei
have defined values when Q+ = Q− = 0 in this new frame:
√ 



2
Q1
2π
R0 + 1/4,
cos ϑi +
3d
3
3 f
√ 



2π
2
Q1
0
e
yi (Q1 , Q+ = 0, Q− = 0) =sin ϑi +R0 + 1/4
,
3d
3
3 f

xe0i (Q1 , Q+ = 0, Q− = 0) =zei0 (Q1 , Q+ = 0, Q− = 0) = 0.

(D.13)(D.14)(D.15)

Hence, the following expansion is valid,e ΛΛ0 =
K

X
ν

X νµ
e ν 0 Qν + 1
e 0 Qν Qµ + . . . ,KKΛΛ2 νµ ΛΛ(D.16)and therefore using the analysis of Longuet-Higgins [136] in which KΛΛ is expanded to third orderin Q1 and second order in ρ, gives:
1
11 2
111 3
e ΛΛ = K
e ΛΛ
e ΛΛ
e ΛΛ
e +− ρ2 .
K
Q1 + K
Q1 + K
Q1 + KΛΛ

(D.17)e ΛΛ±1 = Ke Λ±1Λ = 0,K

(D.18)

1e ΛΛ±2 = Ke Λ±2Λ = Ke ΛΛ±2Kρe∓iφ ,

(D.19)We also haveande ΛΛ0 also obey
to first order in ρ. The matrix elements of Ke ΛΛ0 ≈KZ

−π dr0 f¯εlΛ r0 f¯εlΛ0 r0
Z





0 ∗ e0
0
0
0
0
e
0
× de
ωe YlΛ θe , ϕ
ee Ve Q, re YlΛ θe , ϕee ,(D.20)

122
where we use the tilde notation to stress that these are expressed in terms of coordinates in the
new frame. The rotation to the new frame from our original one has the following effects: θee = θeand ϕee = ϕe − φ/2. Hence a comparison with Eq. (D.2) givese ΛΛ0 = KΛΛ0 ei(Λ−Λ0 )φ/2 .
K(D.21)

Because the reaction matrix K and quantum defect matrix µ are related by µ = µ (Q = 0) +
arctan (K) /π ≈ µ (Q = 0) + K/π + O K 3 /π + . . ., at least to first order, the K matrix and µmatrix have the same symmetry properties, and hence we write the form of µ matrix as Eq.(3.18)
and Eqs.(3.23-3.25).

Appendix E

Effective adiabatic potentials

The effective potentials shown in Fig. 1 (a) of the main text are very complicated, making
identification of the important physics challenging. For that reason, this section presents somedetails of our work that not only give support to our physical interpretation of the nature of thethree-body parameter but also show how some important nonuniversal aspects appear in the hyperspherical adiabatic representation. In fact, Fig. E.1 shows some of the most drastic nonadiabatic
effects found in our calculations. The model proposed in Ref. [16] is also considered, and we showthat this model offers an interesting qualitative picture; but since incorporation of corrections tothat model diminishes its accuracy in the three-body parameter universality, its agreement withexperiment might be fortuitous.
We first consider the validity of the single-channel adiabatic hyperspherical approximationand point out the manner in which some important nonuniversal features manifest themselves.
Figure E.1 shows the results for Uν (R) and Wν (R) obtained from three different two-body potential models. Figures E.1 (a) and (b) show the results for the potential models vλa (λ = λ∗2 ) andvλa (λ = λ∗10 ) [Eq. (5.3)], respectively, while Fig. E.1 (c) shows the results obtained for vsch (D = D6∗ )[Eq. (5.2)]. The most striking aspect of these figures is that Uν (R) and Wν (R) are substantially different, meaning that the nonadiabatic couplings Pνν 0 (R) and Qνν 0 (R) are important near R = rvdW .
Therefore, it is clear that one needs to go beyond a single channel approximation in order to describe the three-body observables. It is worth noting that the nonadiabatic couplings originate fromthe hyperradial kinetic energy. Their large values are thus consistent with our physical picture in

124
4
2
0
-2
-4
0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

4
2
0
-2
-4

4
2
0
-2
-4

Figure E.1: Comparison of Uν (R) (green dashed curves) and Wν (R) (red solid curves) demonstrating the importance of nonadiabatic effects introduced by Qνν (R). (a) and (b) show the results forthe potential models vλa (λ = λ∗2 ) and vλa (λ = λ∗10 ) [Eq. (5.3)], respectively, and in (c) we show theresults obtained for vsch (D = D6∗ ) [Eq. (5.2)]. (a) and (b) also show the effect of the diabatization
scheme used in order to prepare some of our figures in the main text (dash-dotted curves). The
goal of the diabatization is to eliminate the sharp features resulting from Qνν (R) which should notcontribute substantially to the three-body observables. The case shown in (c), however, does notallow us to easily trace the diabatic version of the potentials Wν . In this case, however, the “less”
sharp features have a larger contribution due to the crossing with a three-body channel describing
a collision between a g-wave molecular state and a free atom, giving rise to the anomalous n = 6
point for the vsch model in Fig. 4 of the main text. Although such cases are relatively infrequentin our calculations (and occur mostly for the vsch model), they do illustrate nonuniversal effects
that can affect the three-body parameter. Nevertheless, the three-body observables obtained forcases like the one shown in (c) are still within the 15% range we claimed for the universality of thethree-body parameter.

125
which the local kinetic energy increases and generates the repulsive barrier in our effective potentials. It is also worth mentioning that since the three-body repulsive barrier prevents particles fromapproaching to small distances, the question of whether or not the short range physics actuallychanges as a function of the external magnetic field (as in experiments in ultracold quantum gases)[137] can not be directly answered by observing features related to Efimov physics.