††thanks: Present address: Joint Quantum Institute, University of Maryland and
NIST, College Park, Maryland, 20742

# Universal three-body recombination via resonant $d$-wave interactions

Jia Wang  Department of Physics and JILA, University of Colorado, Boulder, CO
80309, USA  Department of Physics, University of Connecticut, Storrs, CT
06269, USA  J. P. D’Incao  Department of Physics and JILA, University of
Colorado, Boulder, CO 80309, USA  Yujun Wang  Department of Physics and JILA,
University of Colorado, Boulder, CO 80309, USA  Chris H. Greene  Department of
Physics and JILA, University of Colorado, Boulder, CO 80309, USA  Department
of Physics, Purdue University, West Lafayette, IN 47907, USA

###### Abstract

For a system of three identical bosons interacting via short-range forces,
when two of the atoms are about to form a two-body $s$-wave dimer, the Efimov
effect takes place leading to the formation of an infinite number of three-
body (Efimov) states. The lowest Efimov state crosses the three-body break-up
threshold when the $s$-wave two-body scattering length is $a\approx-9.73r_{\rm
vdW}$, $r_{\rm vdW}$ being the van der Waals length. This article focuses on a
generalized version of this Efimov scenario, where two of the atoms are about
to form a two-body $d$-wave dimer, resulting in strong $d$-wave interactions.
Bo Gao has argued, in [Phys. Rev. A. 62, 050702(R) (2000)], that for broad
resonances the $d$-wave dimer is always formed when $a\approx 0.956r_{\rm
vdW}$. Our results demonstrate that a single universal three-body state
associated with the $d$-wave dimer is also formed near the three-body break-up
threshold at $a\approx 1.09r_{\rm vdW}$, or alternatively $a_{2}=0.902r_{\rm
vdW}$, where $a_{2}$ is the two-body $d$-wave scattering length. Such a
universal three-body state is signaled experimentally by an enhancement of the
three-body recombination rate. The three-body effective potential curves that
are crucial for understanding the recombination dynamics are also calculated
and analyzed. An improved method to calculate the couplings, effective
potential curves, and recombination rate coefficients is presented.

##  I Introduction

In recent years, the investigation of Efimov physics Efimov (1971) has
attracted much interest and has sparked a tremendous advance in our
understanding of its fundamental aspects from both experimental Kraemer et al.
(2006); Berninger et al. (2011); Pollack et al. (2009); Gross et al. (2009,
2010); Ottenstein et al. (2008); Lompe et al. (2010); Huckans et al. (2009);
Williams et al. (2009); Wild et al. (2012); Barontini et al. (2009) and
theoretical Petrov (2004); Wang et al. (2011a); Braaten and Hammer (2006);
Hammer and Platter (2007); von Stecher et al. (2009); Wang et al. (2011b)
viewpoints. Efimov physics describes a universal phenomenon in three-body
systems: when the atoms are about to form a two-body $s$-wave dimer, i.e., the
$s$-wave scattering length $a$ goes to infinity, and an infinite number of
three-body Efimov states exists. The ratio between the energies of nearby
Efimov states is given by a universal scaling factor
$E_{n+1}/E_{n}=e^{-2\pi/s_{0}}$ , where $s_{0}$ depends only on the mass
ratios, number of resonant interactions and quantum statistics of the atoms,
e.g., $e^{\pi/s_{0}}\approx 22.7$ for three identical bosons.

Therefore, a single parameter, called the three-body parameter, is needed to
determine the whole Efimov spectrum and other scattering observables. One
possible definition of the three-body parameter is $a_{-}^{*}$, the value of
$a$ at which the first resonance in the three-body recombination rate $K_{3}$
appears, i.e., at which the lowest Efimov state energy crosses the three-body
break-up threshold ($E=0$). One fundamental assumption underlying all Efimov
physics is that the three-body parameter encapsulates the short-range details
of two- and three-body interactions, and thus should not be universal.
Surprisingly, however, ultracold experiments with alkali atoms near Fano-
Feshbach resonances Kraemer et al. (2006); Berninger et al. (2011); Pollack et
al. (2009); Gross et al. (2009, 2010); Ottenstein et al. (2008); Lompe et al.
(2010); Huckans et al. (2009); Williams et al. (2009); Wild et al. (2012);
Barontini et al. (2009) have observed a universal value, namely
$a_{-}^{*}\approx-9.1r_{\rm vdW}$ in homonuclear atomic systems, where $r_{\rm
vdW}$ is the van der Waals length. The van der Waals length $r_{\rm
vdW}=\left({2\mu_{2b}C_{6}}\right)^{1/4}/2$ is a characteristic length scale
for the van der Waals interaction $-C_{6}/r^{6}$ between two neutral atoms
with two-body reduced mass $\mu_{2b}$. A relevant energy scale, called the van
der Waals energy, can be defined as $E_{\rm
vdW}=\hbar^{2}/\left({2\mu_{2b}r_{\rm vdW}^{2}}\right)$.

This newly found universality in the three-body parameter was subsequently
studied in several different theoretical models Schmidt et al. (2012); Chin
(2011); Sørensen et al. (2012); Wang et al. (2012a); Naidon et al. (2012). In
particular, Refs. Wang et al. (2012a); Naidon et al. (2012) conclude that the
universality of the three-body parameter comes from a universal barrier in the
three-body potential curves arising from interatomic distances around 2
$r_{\rm vdW}$. This barrier originates from the suppression of the probability
to find two atoms at distances shorter than $r_{\rm vdW}$ due to the increase
of the classical local velocity and it is a universal property of van der
Waals interactions between neutral atoms. Our previous studies Wang et al.
(2012a, b) have indicated that the Lennard-Jones potential is an excellent
model potential to study the universality of three-body physics in ultracold
atomic gases, and for this reason it is also adopted in this analysis.

The present study, however, employs the Lennard-Jones potential model to study
three-body systems with atom pairs close to a $d$-wave resonance.
Interestingly, Gao Gao (2000) predicted that $d$-wave dimers form at a
universal value of $a$, $a\approx 0.956r_{\rm vdW}$, for potentials with a van
der Waals tail. Evidently, since that prediction was obtained for a single
channel potential model its validity is expected to hold only for broad Fano-
Feshbach resonances, and possibly is limited only to single-channel two-body
interactions. Nevertheless, one natural question is: will there also be a
universality of three-body physics near that value of scattering length, i.e.,
when $d$-wave interactions are expected to be strong? Here, the three-body
recombination rate is studied to answer this question, and it is seen to
exhibit two universal enhancements. Analysis of the three-body effective
potential curves yields an intuitive understanding of these enhancements. One
of these enhancements is caused by the simple fact that when the two-body
$d$-wave state becomes barely bound an additional decay channel is formed,
leading to an enhancement in recombination to $d$-wave dimers. The second
enhancement, however, has a different nature. It corresponds to the formation
of a universal three-body state that crosses the collision threshold for three
free particles. Our analysis indicates that other three-body states associated
with excited 2-body angular momenta are also possible, however, their
signatures in recombination will be suppressed in the low energy limit due to
the Wigner threshold laws for recombination Esry et al. (2001). It is
important to mention that since in experiments using ultracold atoms $s$\- and
$d$-wave states are coupled by the magnetic anisotropic dipole interaction the
$s$-wave scattering actually diverges when the $d$-wave dimer becomes bound.
Although the connection between these three-body states associated with
$d$-wave resonant interactions and the multichannel physics in ultracold atoms
still remains to be understood more deeply, the present study can offer an
alternative parameterization in terms of the properties of the $d$-wave
scattering length. Moreover, since the three-body effective potential curves
are an important tool for understanding universality in few-body systems at
ultracold temperatures Wang et al. (2012a); Naidon et al. (2012); Wang et al.
(2012b), an improved method to calculate the coupling matrix elements is also
developed in this paper.

The remainder of the paper is organized as follows. Section II tests the Gao
prediction numerically for a Lennard-Jones potential, and demonstrates very
good agreement. The three-body recombination rate is then calculated using the
Lennard-Jones two-body potential with 2 $s$-wave bound states in the vicinity
of the $d$-wave resonance. The improved numerical method is discussed in the
subsection II.A. Analysis of the recombination rate using effective potential
curves is made in subsection II.B. Section III checks the universality using
the Lennard-Jones potential with 3 $s$-wave dimer bound states.

##  II Universal two-body $d$-wave dimer

One of the key experimental tools used to study ultracold atomic gases is the
Fano-Feshbach resonance, which can be utilized to magnetically tune the
scattering length. Although the multichannel character of the hyperfine
interactions leads to a great deal of complication, the single-channel van der
Waals interactions has been shown to offer a good model for studying a broad
Fano-Feshbach resonance in ultracold atomic gases Chin et al. (2011). In 2000,
Bo Gao predicted that for single-channel interactions with a van der Waals
tail, $-C_{6}/r^{6}$, there is always a $d$-wave dimer (and dimers with higher
angular momentum $l=4j+2$, where $j=1,2,3...$) that becomes bound at a
universal value of the two-body $s$-wave scattering length
$a=a^{*}=4\pi/[\Gamma(1/4)]^{2}r_{\rm{vdW}}\approx 0.956r_{\rm{vdW}}$ Gao
(2000). [Note that $a^{*}$ has the same value as the so-called mean scattering
length $\bar{a}$ as defined by Gribakin and Flambaum Gribakin and Flambaum
(1993); Chin et al. (2011), which in turn should not be confused with the
usual scattering length $a$.] Another special property of the two-body system
at $a=a^{*}$ is that the effective range $r_{\rm{eff}}$ has the value
$2\left[{\Gamma\left({1/4}\right)}\right]^{2}/\left({3\pi}\right)r_{\rm{vdW}}\approx
2.789r_{\rm{vdW}}$, which equals the effective range when $s$-wave scattering
length is infinity. In Gao’s work Gao (2000), the universality of $a^{*}$ is
based on an $l$-independent quantum-defect theory, where the short-range
physics can be described by a single short-range parameter, $K^{c}$, which is
itself approximately $l$-independent under proper conditions that are
satisfied by most diatomic systems. In combination with solutions for the
long-range interactions, $K^{c}$ uniquely determines both the scattering and
the bound state properties of diatomic systems.

Figure 1: (Color online) The values of the two-body $s$-wave scattering length
$a^{*}_{l}$ at the point where a $d$-wave ($l=2$) dimer (black curve with
square symbols) and where an $i$-wave ($l=6$) dimer just becomes bound (red
curve with circular symbols) are shown as functions of the number of two-body
$s$-wave bound states.

In this section, the universal value of the two-body $s$-wave scattering
length $a_{l}^{*}$ [where a $d$-wave ($l=2$) or $i$-wave ($l=6$) dimer becomes
bound] is studied numerically as a function of the number of two-body $s$-wave
bound states and using the Lennard-Jones potential model for the two-body
interaction:

$v(r)=-\frac{C_{6}}{r^{6}}\left(1-\frac{\lambda^{6}}{r^{6}}\right).$ (1)

Note that in the present work the parameter $\lambda$ is used to adjust the
values of $a$ as well as the number of bound states. (Some additional details
about the Lennard-Jones potential are discussed in Appendix A.) In Fig. 1, the
black square symbols (red circular symbols) indicates the numerical values of
the two-body $s$-wave scattering length $a^{*}_{2}$ ($a^{*}_{6}$) at the point
where a $d$-wave ($i$-wave) dimer becomes bound. The more $s$-wave bound
states supported by the system, the deeper is the potential, implying
different short-range physics. The universal prediction from Gao’s work,
$a^{*}\approx 0.956r_{\rm vdW}$ Gao (2000) (horizontal dashed line in Fig. 1),
however, suggests that the actual form of the potential at short-distances is
not important. In fact, our numerical results obtained for $a^{*}_{l}$ agree
well with the Gao prediction of two-body binding energies [within $1\%$
($6\%$) in the case of 10 $s$-wave bound states for a $d$-wave ($i$-wave)] and
this agreement tends to improve as one increases the number of $s$-wave
states. In the vicinity of the point where the dimer is just barely bound,
i.e., when $a\leq a^{*}_{l}$, the binding energy can be expressed as a linear
function of the scattering length , i.e.,

$E_{l}/\left[{\hbar^{2}/(2\mu_{2b}{r_{\rm{vdW}}^{2}})}\right]\approx
d_{l}\left({a^{*}_{l}-a}\right)/r_{\rm{vdW}},$ (2)

where $d_{2}\approx 5.8$ and $d_{4}\approx 43$, for $d$-wave and $i$-wave
states, respectively, are approximately universal. It is also interesting to
note that near $a=a^{*}_{2}$ one expects strong $d$-wave interactions. In
order to parametrize $d$-wave interactions in the ultracold regime, scattering
properties must also be explored for two atoms colliding in the $l=2$ channel.
In this case, due to the $1/r^{6}$ dependence of the Lennard-Jones potential,
the $d$-wave elastic scattering phase-shift $\delta_{2}\left(k\right)$ is
known to have the following low energy expansion:

$\tan\delta_{2}\left(k\right)=-\lambda_{2}k^{4}-a_{2}^{5}k^{5},$ (3)

instead of the usual $k^{2l+1}$ dependence found for $l<2$. Here,
$\lambda_{2}=\left[{\pi\Gamma\left(5\right)\Gamma\left({1/2}\right)}\right]/\left[{4\Gamma^{2}\left(3\right)\Gamma\left({11/2}\right)}\right]r_{\rm
vdW}^{4}\approx 0.160r_{\rm vdW}^{4}$ is a constant that only depends on the
van der Walls length, and $a_{2}$ diverges when a $d$-wave state is just about
to be bound Willner and Gianturco (2006). Therefore, $a_{2}$ is denoted the
“$d$-wave scattering length” in this study (which should not be confused with
$a_{2}^{*}$), and is also used in order to characterize our findings in terms
of the $d$-wave interactions. This parameterization, therefore, allows a
better comparison to experiments in ultracold quantum gases. It is important,
however, to emphasize that the complicated multichannel structure of alkali
atoms used in ultracold quantum gases can strongly affect the relations
explored above between the three-body states found for $d$-wave interactions
and the single channel parameters that characterize the present study. This is
particularly true near narrow Feshbach resonances Wang et al. (2011a).

##  III Three-body recombination

This section focuses on the three-body recombination rate near the region
where the $d$-wave dimer is about to be bound. In particular, the example of a
Lennard-Jones potential with two $s$-wave bound states is elaborated. In this
numerical example, the $d$-wave dimer becomes bound at about
$a=0.995r_{\rm{vdW}}$ and the $i$-wave dimer becomes bound at about
$1.206r_{\rm{vdW}}$. Later, the three-body recombination rate with more two-
body $s$-wave bound states will be computed to test the universality of those
results.

###  III.1 Hyperspherical approach

The hyperspherical approach is utilized here to study the three-body systems.
The version of that approach adopted here has been discussed in detail in
other recent studies Wang et al. (2011c, 2012a). To briefly summarize, the
three-body system is described in hyperspherical coordinates: a hyperradius
$R$ that describes the overall size of the three-body system, two hyperangles
$\theta$ and $\phi$ that describe the shape of the three-body system, and the
three usual Euler angles $\alpha,\beta,\gamma$ that describe rigid-body
rotations. In these coordinates, the three-body Schrödinger equation is
written as,

$\left[{-\frac{{\hbar^{2}}}{{2\mu_{3b}}}\left({\frac{{\partial^{2}}}{{\partial
R^{2}}}-\frac{{\Lambda^{2}+15/4}}{{R^{2}}}}\right)+V\left({R,\theta,\varphi}\right)-E}\right]\psi_{E}=0,$
(4)

where $\Lambda^{2}\equiv\Lambda^{2}(\Omega)$
[$\Omega\equiv\\{\theta,\varphi,\alpha,\beta,\gamma\\}$] is the “grand
angular-momentum operator” and $\mu_{3b}=m/\sqrt{3}$ is the three-body reduced
mass of three identical atoms with mass $m$ Kendrick et al. (1999); Suno et
al. (2002).

Figure 2: (Color online) Adiabatic three-body potential curves for $a\approx
0.977r_{\rm vdW}$. The $2+1$ channels, which corresponding to a dimer plus a
free atom at very large distance, are labeled by a combination of a letter and
a number. The letter denotes the angular momentum quantum number $l$ of the
dimer, and the number labels the channels for the same dimer angular momentum
from low-to-high dimer binding energies.

To solve this Schrödinger equation, the first step in our method is to
diagonalize the adiabatic Hamiltonian $H_{\rm{ad}}\left({R,\Omega}\right)$ and
find the adiabatic potentials $U_{\nu}(R)$ and channel functions
$\Phi_{\nu}\left({R;\Omega}\right)$:

$H_{\rm{ad}}\left({R,\Omega}\right)\Phi_{\nu}\left({R;\Omega}\right)=U_{\nu}\left(R\right)\Phi_{\nu}\left({R;\Omega}\right),$
(5)

whose solutions depend parametrically on $R$. The adiabatic Hamiltonian,
containing all angular dependence and interactions, is defined as

$H_{\rm{ad}}\left({R,\Omega}\right)=\left[{\frac{{\hbar^{2}\Lambda^{2}}}{{2\mu_{3b}R^{2}}}+\frac{{15\hbar^{2}}}{{8\mu_{3b}R^{2}}}+V\left({R,\theta,\varphi}\right)}\right].$
(6)

For illustration, Fig. 2 shows the three-body adiabatic hyperspherical
potential curves for the case of $a\approx 0.977r_{\rm{vdW}}$. [Note that in
this figure the labels for the final recombination channels are also indicated
and are useful for further analysis.]

The adiabatic potentials $U_{\nu}(R)$ and channel functions
$\Phi_{\nu}(R;\Omega)$ obtained by solving Eq. (5) for fixed values of $R$
contain all the correlations relevant to this problem. For each $R$, the set
of $\Phi_{\nu}\left({\Omega;R}\right)$ is orthogonal,

$\int{d\Omega\Phi_{\mu}\left({R;\Omega}\right)^{*}}\Phi_{\nu}\left({R;\Omega}\right)=\delta_{\mu\nu},$
(7)

and complete

$\sum\limits_{\tau}{\Phi_{\tau}\left({R;\Omega}\right)\Phi_{\tau}\left({R;\Omega^{\prime}}\right)^{*}=\delta\left({\Omega-\Omega^{\prime}}\right)}.$
(8)

In practice, calculation of all the channel functions is time consuming and
impractical. However, numerical studies have shown that only a small number of
channels are needed as a truncated basis-set to expand the total wave
function, e.g.,

$\psi_{E}\left({R,\Omega}\right)=\sum\limits_{\nu=1}^{N_{c}}{F_{\nu}^{E}\left(R\right)\Phi_{\nu}\left({\Omega;R}\right)},$
(9)

where $N_{c}$ is the number of channels adopted. Insertion of Eq. (9) into Eq.
(4) leads to a set of coupled one-dimensional equations:

$\left[{-\frac{\hbar^{2}}{{2\mu_{3b}}}\frac{d^{2}}{{dR^{2}}}+U_{\nu}\left(R\right)}-E\right]F_{\nu\nu^{\prime}}^{E}\left(R\right)-\frac{\hbar^{2}}{{2\mu_{3b}}}\sum\limits_{\mu}{\left[{2P_{\nu\mu}\left(R\right)\frac{d}{{dR}}+Q_{\nu\mu}\left(R\right)}\right]}F_{\mu\nu^{\prime}}^{E}\left(R\right)=0,$
(10)

where $\underline{P}$ and $\underline{Q}$ are the coupling matrices defined
below, and $\nu^{\prime}$ denotes the $\nu^{\prime}-$th independent solution.
(Hereafter, unless otherwise specified, we use an underline to denote the
matrix form, e.g., $\underline{P}$ denotes a matrix with matrix element
$P_{\nu\mu}$.) In the equation above, the non-adiabatic coupling matrices are
defined as

$P_{\nu\mu}\left(R\right)=\int{d\Omega}\Phi_{\nu}\left({R;\Omega}\right)^{*}\frac{\partial}{{\partial
R}}\Phi_{\mu}\left({R;\Omega}\right),$ (11)
$Q_{\nu\mu}\left(R\right)=\int{d\Omega}\Phi_{\nu}\left({R;\Omega}\right)^{*}\frac{{\partial^{2}}}{{\partial
R^{2}}}\Phi_{\mu}\left({R;\Omega}\right),$ (12)

and are responsible for inelastic scattering processes as well as the finite
width of three-body resonant states. In practice, only the

$P_{\nu\mu}^{2}\left(R\right)=-\int{d\Omega}\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R;\Omega}\right)^{*}\frac{\partial}{{\partial
R}}\Phi_{\mu}\left({R;\Omega}\right),$ (13)

component of $Q_{\nu\mu}$ is needed to solve the coupled equations Burke
(1999). The relation between $\underline{Q}$ and $\underline{P}^{2}$ is given
by $\frac{d}{{dR}}\underline{P}=-\underline{P^{2}}+\underline{Q}$. From the
definition of the $\underline{P}$ and $\underline{Q}$ matrices, it is easy to
see that the coupling matrices have the following properties:
$P_{\nu\mu}=-P_{\mu\nu}$ and $P_{\nu\mu}^{2}=P_{\mu\nu}^{2}$, which leads to
$P_{\nu\nu}=0$, and $Q_{\nu\nu}=P_{\nu\nu}^{2}$. In the hyperspherical
adiabatic representation, the effective potentials are usually defined as

$\widetilde{U}_{\nu}\left(R\right)\equiv
W_{\nu\nu}\left(R\right)=U_{\nu}\left(R\right)-\frac{{\hbar^{2}}}{{2\mu_{3b}}}P_{\nu\nu}^{2}\left(R\right),$
(14)

and differs for the hyperspherical potential $U_{\nu}(R)$ by the addition of
the diagonal correction $P_{\nu\nu}^{2}(R)$. The effective potentials are
usually more physical than the raw Born-Oppenheimer potentials because they
include hyperadial kinetic energy contributions through the $P_{\nu\nu}^{2}$
term. For example, the effective potential gives physical asymptotic behaviors
of the system at large $R$ with the correct coefficient of $R^{-2}$. In some
cases, one can even neglect the nonadiabatic couplings between different
adiabatic channels and approximately solve the Schrödinger equation as a
single channel problem:

$\left[{-\frac{{\hbar^{2}}}{{2\mu_{3b}}}\frac{{d^{2}}}{{dR^{2}}}+\tilde{U}_{\nu}\left(R\right)-E}\right]F_{\nu\nu}^{E}\left(R\right)=0,$
(15)

called the adiabatic approximation Kokoouline and Greene (2003). Although a
semi-quantitative picture can emerge within this approximation, the fully
coupled system of equation has to be solved to obtain a more quantitative
description of the properties of the three-body system including inelasticity
in particular.

Quantitative solutions of Eq. (10) hinge on the ability to calculate non-
adiabatic couplings accurately. A traditional method for calculating the
coupling matrices is to apply a simple finite difference scheme for the
derivative of $\Phi_{\mu}\left({R;\Omega}\right)$, i.e.,

$\frac{\partial}{{\partial
R}}\Phi_{\mu}\left({R;\Omega}\right)\approx\frac{{\Phi_{\mu}\left({R+\Delta
R;\Omega}\right)-\Phi_{\mu}\left({R-\Delta R;\Omega}\right)}}{{2\Delta R}}.$
(16)

This scheme, however, is only accurate up to the first order in $\Delta R$.
Here, we use instead an improved method to calculate the
$\partial\Phi_{\mu}\left({R;\Omega}\right)/\partial R$. This method only
relies on $H_{\rm{ad}}\left({R,\Omega}\right)$ and $\partial
H_{\rm{ad}}\left({R,\Omega}\right)/\partial R$ at the desired value of $R$ and
is shown to be in principle numerically exact (see Appendix B for a detailed
discussion).

###  III.2 Three-body recombination rate

The numerical method described in the previous subsection yields more accurate
coupling matrix elements and hence improved effective potential curves. These
effective potentials are not only crucial for intuitively understanding the
system, but also for quantitative calculations of three-body observables like
the recombination rate. Together with the $R$-matrix propagation method
discussed in Ref. Wang et al. (2011c), the three-body recombination rate for
$J^{\pi}=0^{+}$ symmetry, where $J$ is the three-body total angular momentum
and $\pi$ the parity, are calculated and shown in Fig. 3. The solid black
curve with square symbols shows the total recombination rate as a function of
the scattering length $a$, near the $d$-wave resonance, i.e., for values of
$a$ near $a^{*}_{2}$. The total rate exhibits two clear enhancements labeled
by “A” ($a_{\rm A}\approx 1.09r_{\rm{vdW}}$) and “B” ($a_{\rm B}\approx
0.98r_{\rm{vdW}}$). The enhancements “A” and “B”, indicated in Fig. 3 by the
vertical dashed lines, occur right after and right before the $d$-wave dimer
becomes bound (see vertical solid line in Fig. 3), respectively. Analysis of
the partial recombination rates helps to elucidate the mechanism responsible
for these enhancements. The enhancement around the peak “A” appears in all the
partial three-body recombination rates. After the $d$-wave becomes bound and
forms a new decay channel, the partial $K_{3}$ rates near the enhancement “B”
have a different characteristic. The total rate is now dominated only by the
partial rate into this new $d$-wave channel [see the inset of Fig. 3]. These
features suggest that peak “A” might be a resonance due to a three-body state
that crosses the three-body threshold, while peak “B” seems to relate to some
other threshold behavior related to the formation of the new $d$-wave decay
channel. The effective potentials near these enhancements, however, can yield
insights into their origin.

Figure 3: (Color online) Total and partial $J^{\pi}=0^{+}$ three-body
recombination rates as functions of the two-body scattering length $a$, shown
on a linear scale, near a $d$-wave resonance associated with the $d$-wave
dimer formation. The inset shows the same graph with a logarithmic scale for
the y-axis. The units of $K_{3}$ are converted to $\rm{cm^{6}/s}$ by using the
van der Waals length $r_{\rm vdW}=101.0$ bohr and mass $132.905429$ amu of
133Cs. The solid vertical line indicates the value of the $s$-wave scattering
length at the point where the $d$-wave dimer becomes bound. The two dashed
vertical lines indicate the two values of $a$ at which recombination rate is
enhanced, denoted A and B respectively.

First, consider the physical origin of the enhancement A. Figure. 4 (a) shows
the adiabatic hyperspherical potential curves at peak A. Only the channels
relevant to the resonances are shown in this figure. The red dashed-line shows
the effective hyperspherical potential. We have diabatized the potential near
an avoided crossing since it plays no role, and display it as the black solid
curve. The black solid curve shows an outer barrier, and a three-body short-
ranged state could be supported for this potential. 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:

$\phi_{\rm
WKB}=\int_{R_{a}}^{R_{b}}{\frac{\hbar}{{\sqrt{-2\mu_{3b}W_{\nu\nu}\left(R\right)}}}}dR,$
(17)

where $R_{a}$ and $R_{b}$ are the classical turning points shown in Fig. 4
(a). The calculated WKB phase $\phi_{\rm WKB}$ is about $0.51\pi$, which is
strong evidence for the existence of a three-body state. [Below, we calculate
the actual energy of this state and describe some of its properties.] It is
this three-body state associated with the $d$-wave dimer state crossing the
three-body break-up threshold that causes the enhancement A in the three-body
recombination rate. The physical picture is that the three-body state at
threshold can trap the system within the potential well for a long time, hence
making it more likely that the system will decay into deeper channels.
Therefore, this decay process only depends on the initial channel, which is
consistent with the fact noted above that near the enhancement A all partial
rates display the resonance effect.

Next consider the enhancement mechanism responsible for peak B. Figure 4 (b)
shows the effective potentials $W_{\nu\nu}(R)$ as functions of $R$ for three
different scattering lengths around the enhancement peak B as solid curves,
dashed curves, and dash-dotted curves. The corresponding $K_{3}$ rates are
indicated by the three points $(i)$, $(ii)$, and $(iii)$ in the inset of Fig.
4. Here, the dominant pathway at low energy is from the lowest entrance
channel to the newly formed $d$-wave channel. For each scattering length, the
thin curves in Fig. 4 are the entrance channels, and as one can readily see,
they do not change substantially in all three cases. This means that the
incident channel potential curve depends on the atom-atom scattering length
only weakly. The thick curves in Fig. 4 are the most important recombination
channel associated with the new $d$-wave channel. It is clear that from point
$(i)$ to $(iii)$, the potential barrier decreases appreciably and can explain
the increase of recombination as it gradually allows the three-body wave
function to approach shorter distances, where non-adiabatic couplings are
expected to be stronger.

From Ref. Wang et al. (2011c), the scattering matrix element describing the
recombination from the lowest incident channel to the newly formed $d$-wave
channel can be approximated by,

$S_{fi}\approx 2\pi
i\int_{0}^{\infty}{dRF_{f}\left(R\right)W_{fi}\left(R\right)F_{i}\left(R\right)}$
(18)

where $F_{i}$ ($F_{f}$) are the regular solutions of Eq. (15) — single channel
approximation — for the initial (final) channel $i$ ($f$). As mentioned above,
the incident channel is dominated by the centrifugal barrier which, in turn,
is independent of the scattering length. Hence the energy-normalized $F_{i}$
are well approximated by Bessel functions:

$F_{i}\left(R\right)\approx\sqrt{\mu_{3b}R}J_{2}\left({k_{i}R}\right),$ (19)

which is, of course, also independent of scattering length. The coupling
matrix element $W_{fi}\left(R\right)$ is also found, numerically, to depend on
scattering length only weakly for the three cases we studied here. Therefore,
we assume that the product of $W_{fi}\left(R\right)F_{i}\left(R\right)$ does
not depend appreciably on the scattering length. In addition, this product is
shown numerically to be significant only at small hyperradius (around 1.3
$r_{\rm vdW}$ to 2.8 $r_{\rm vdW}$). Therefore, the scattering length
dependence of $S_{fi}$ is fully determined by $F_{f}\left(R\right)$ at short
distances. There are two competing mechanisms controlling the amplitude of
$F_{f}\left(R\right)$: one is the tunneling through the outer barrier, which
makes the amplitude become increasingly larger, as the barrier becomes
increasingly lower from $(i)$ to $(iii)$; the other is the overall amplitude
for an energy normalized wave function, which makes the amplitude becomes
increasingly smaller, since the threshold energy associated to the energy of
the $d$-wave dimer become lower and lower, and the out-going wave vector
become larger and larger. These two competing mechanisms help to form peak
structure B. The WKB phase calculated for these three cases is smaller than
$\pi/2$, which is another evidence that this enhancement is not due to a
three-body state resonance. This is also consistent with the fact that such an
enhancement is not shown in the partial recombination rates for channels other
than the newly formed $d$-wave.

Figure 4: (a) (Color online) $W_{\nu\nu}(R)$ as a function of hyperradius $R$
at the $K_{3}$ resonance peak A. (b) (Color online) $W(R)_{\nu\nu}$ as a
function of hyperradius $R$ around the $K_{3}$ enhancement peak B. The inset
shows the three points where the potential curves correspond to.

##  IV Three-body state associated with the $d$-wave dimer

While the Efimov states can be viewed as three-body states associated with a
$s$-wave dimer near the three-body breakup threshold, the new three-body state
discussed in the previous section is a three-body state that can 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 5
shows the three-body state energy as a function of the scattering length $a$
(red dots), with the error bars indicating the width of the quasi-bound state.
The red line represents a fitting formula to this energy:

$\frac{E_{3b}}{\hbar^{2}/\left({2\mu_{2b}r_{\rm{vdW}}^{2}}\right)}=d_{3b}\frac{\left({a^{*}_{3b}-a}\right)}{r_{\rm
vdW}}+e_{3b}\frac{\left({a^{*}_{3b}-a}\right)^{2}}{r_{\rm vdW}^{2}},$ (20)

where $a^{*}_{3b}\approx 1.09r_{\rm{vdW}}$, $d_{3b}\approx 8.21$ and
$e_{3b}\approx 10.01$ are fitting parameters. For comparison, the figure also
shows the $d$-wave energy as the black curve with solid square symbols. The
three-body energy crosses the three-body breakup threshold at $a\approx
1.09r_{\rm{vdW}}$, which corresponds to peak A in $K_{3}$. As one can see, the
three-body state is formed before the $d$-wave state becomes bound and its
energy is always below the $d$-wave dimer energy (the black curve with square
symbols in Fig. 5). The difference between the dimer and trimer energy also
increases when the scattering length becomes smaller, and the dimer becomes
deeper. However, the width of the three-body quasi-bound state does not seem
to change appreciably.

Figure 5: (Color online) Energy of the three-body bound state associated with
a $d$-wave dimer as a function of scattering length $a$, both in van der Waals
units.

The Gao analysis in Ref. Gao (2000), suggests that there is also an $i$-wave
dimer that becomes bound at 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
$K_{3}$ rate into the $i$-wave dimer channel is negligible compared with the
other partial rates. One explanation may be that forming a $i$-wave ($l=6$)
dimer results in the exchange of a large amount of angular momentum between
the dimer and the free 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 rates
at ultracold energies.

Figure 6: (Color online) The enhancements for the total three-body
recombination rates at about $a=0.995r_{\rm vdW}$ for Lennard-Jones potential
with 2 and 3 $s$-wave bound states. $K_{3}$ is convert to $\rm{cm^{6}/s}$ by
using van der Waals length $r_{\rm vdW}=101.0$ bohr and mass $132.905429$ amu
of 133Cs Figure 7: (Color online) The total three-body recombination rates as
a function of $r_{\rm vdW}^{5}/a_{2}^{5}$. The units of $K_{3}$ are the same
as Fig. 6.

Finally, to test the universality of our results, the three-body recombination
enhancements for the model potentials of Lennard-Jones type, having either 2
or 3 $s$-wave bound states, are shown in Fig. 6. In this figure, the red
(black) solid vertical line shows where the $d$-wave dimer crosses the
threshold for LJ with 2 (3) $s$-waves. The red (black) dashed vertical line
indicates where the three-body state associated with the $d$-wave level
crosses the threshold for LJ with 2 (3) $s$-wave states. The values of the
$s$-wave scattering length for such enhancements from these two different
models differs by only 0.01 $r_{\rm vdW}$. This small difference is consistent
with the small shift of the scattering length value where the $d$-wave dimer
crosses the three-body threshold (see Fig. 1). As discussed before, the
$d$-wave scattering length is a better parameter to describe $d$-wave
interactions. We have therefore analyzed the three-body recombination rate as
a function of $1/a_{2}^{5}$, which is shown in Fig. 7. The agreement between
calculations using different potential models in fact looks better in this
parameterization, where the three-body state crosses the threshold at
$a_{2}\approx 0.902r_{\rm vdW}$ for both LJ with 2 and 3 $s$-waves bound
states.

This result suggests (although it does not prove rigorously) that the three-
body state associated with the $d$-wave dimer is also universal. We also note
that at the threshold for formation of the three-body state, the effective
potential relevant for that state does also display an inner repulsive barrier
[see Fig. 4 (a)] preventing atoms from approaching to hyperradii smaller than
$1.5r_{\rm vdW}$. This indicates that the properties of such a three-body
state should also be insensitive to the details of nonadditive short-range
three-body interactions.

##  V Summary

In summary, we have calculated the three-body recombination rates for Lennard-
Jones potentials with 2 and 3 $s$-wave bound states for small positive values
of the scattering length. Two universal enhancement peaks are found at about
$a=0.995$ $r_{\rm vdW}$. In particular, one of the enhancement peaks
corresponds to a new universal three-body state that is associated with a
$d$-wave dimer and is expected to exist in general for potentials having a van
der Waals tail.

## Acknowledgements

The authors want to thank Eric A. Cornell, John Bohn, Cheng Chin, Robin Côté
and Doerte Blume for helpful discussions. This work was supported in part by
NSF and the AFOSR-MURI.

##  Appendix A Lennard-Jones potential

The two-body model potentials applied here are the Lennard-Jones potential,
which can be expressed as a function of interatomic distances $r$:

$V\left(r\right)=-\frac{{C_{6}}}{{r^{6}}}\left({1-\frac{{\lambda^{6}}}{{r^{6}}}}\right).$
(21)

This form of potential has a van der Waals tail $-C_{6}/r^{6}$ at long range
and a repulsive core near the origin controlled by the parameter $\lambda$. In
the present study, the values of $\lambda$ are adjusted to give the desired
$s$-wave scattering length $a$ and number of bound states. Figure 8 shows two
Lennard-Jones potentials with a zero-energy $d$-wave state and their
corresponding spectrum of bound state with even $l$. (The zero energy $d$-wave
state is not shown in Fig. 8). The black thicker curve corresponds to
$\lambda\approx 0.503r_{\rm vdW}$ and has two $s$-wave bound states with
$a\approx 0.995r_{\rm vdW}$. The red thinner curve corresponds to
$\lambda\approx 0.414r_{\rm vdW}$ and has three $s$-wave bound states with
$a\approx 0.983r_{\rm vdW}$. It is shown that although the potentials at
short-range are very different for the two different parameters, the
scattering length and the spectrum near threshold are similar.

Figure 8: (Color online) Lennard-Jones potential with $\lambda\approx
0.503r_{\rm vdW}$ (black thicker curve) and $\lambda\approx 0.414r_{\rm vdW}$
(red thinner curve). The horizontal lines show the bound state spectrum with
even angular momentum quantum number $l$. The line style solid, dashed, dash-
dotted, dash-dot-dotted, and short-dashed indicates $l=0,2,4,6,8$
correspondingly.

##  Appendix B Improved method for calculating nonadiabatic coupling matrices

Equation (16) is only accurate up to the first order in $\Delta R$. In
addition, the value chosen for $\Delta R$ in a realistic numerical calculation
can sometimes be tricky. When $\Phi_{\mu}\left({R;\Omega}\right)$ changes very
rapidly, e.g., near a sharp avoided crossing, we need to choose a very small
step size $\Delta R$. In contrast, when $\Phi_{\mu}\left({R;\Omega}\right)$
changes very slowly, e.g., at very large distances $R$, we need to choose a
relatively larger step size $\Delta R$, or else $\Phi_{\mu}\left({R+\Delta
R;\Omega}\right)-\Phi_{\mu}\left({R-\Delta R;\Omega}\right)$ would be too
small, and the accuracy would be limited by the machine precision.

One way to improve the accuracy is to apply the Hellmann-Feynman theorem. This
theorem can give us analytical formulas for the coupling matrices if we know
the derivative of the adiabatic Hamiltonian $\frac{\partial}{{\partial
R}}H_{\rm{ad}}$ from the following derivation. First, taking the derivative of
both sides of Eq. (5) leads to

$\displaystyle\left[{H_{\rm{ad}}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}\right]\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R;\Omega}\right)=$
$\displaystyle-\left[{\frac{\partial}{{\partial
R}}H_{\rm{ad}}\left({R,\Omega}\right)-\frac{\partial}{{\partial
R}}U_{\nu}\left(R\right)}\right]\Phi_{\nu}\left({R;\Omega}\right).$ (22)

Next, multiplying $\Phi_{\mu}\left({R;\Omega}\right)^{*}$ on both sides of Eq.
(B) and integrating over $\Omega$ gives

$\displaystyle P_{\mu\nu}$ $\displaystyle=$
$\displaystyle\int{d\Omega\Phi_{\mu}\left({R;\Omega}\right)^{*}}\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R;\Omega}\right)$ (23) $\displaystyle=$
$\displaystyle-\frac{\int{d\Omega\Phi_{\mu}\left({R;\Omega}\right)^{*}\left[{\frac{\partial}{{\partial
R}}H_{\rm{ad}}\left({R,\Omega}\right)}\right]\Phi_{\nu}\left({R;\Omega}\right)}}{\left[{U_{\mu}\left(R\right)-U_{\nu}\left(R\right)}\right]},$

where $\mu\neq\nu$, and

$\frac{\partial}{{\partial
R}}U_{\nu}\left(R\right)=\int{d\Omega\Phi_{\nu}\left({R;\Omega}\right)^{*}}\left[{\frac{\partial}{{\partial
R}}H_{\rm{ad}}\left({R,\Omega}\right)}\right]\Phi_{\nu}\left({R;\Omega}\right)$
(24)

after some manipulation of algebra. The Hellmann–Feynman theorem is in
principle exact. The matrix elements for $P^{2}$ can be obtained by

$P_{\mu\nu}^{2}=\sum\limits_{\tau=1}^{N_{c}}{P_{\mu\tau}P_{\tau\nu}},$ (25)

where $N_{c}$ is the number of channels. However, numerical studies show that
the convergence of $P_{\mu\nu}^{2}$ with respect to number of channels is very
slow, making this method impractical.

We now introduce an improved method to calculate $\frac{\partial}{{\partial
R}}\Phi_{\mu}\left({R;\Omega}\right)$. It is in principle numerically exact.
The first hint of the derivation of this method is that Eq. (B) seems
plausible to be directly solved for $\frac{\partial}{{\partial
R}}\Phi_{\mu}\left({R;\Omega}\right)$ by

$\displaystyle\frac{\partial}{{\partial R}}\Phi_{\nu}\left({R;\Omega}\right)=$
$\displaystyle-\left[{H_{\rm{ad}}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}\right]^{-1}\times$
$\displaystyle\left[{\frac{\partial}{{\partial
R}}H_{\rm{ad}}\left({R,\Omega}\right)-\frac{\partial}{{\partial
R}}U_{\nu}\left(R\right)}\right]\Phi_{\nu}\left({R;\Omega}\right).$

However, this solution is forbidden since
${H_{\rm{ad}}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}$ is singular:
$\left|{H_{\rm{ad}}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}\right|=0$,
meaning that ${H_{\rm{ad}}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}$ is
not invertible. The singularity can also be understood from the fact that the
equation

$\displaystyle\left[{H_{ad}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}\right]\chi_{\nu}\left({R;\Omega}\right)=$
(27) $\displaystyle-\left[{\frac{\partial}{{\partial
R}}H_{ad}\left({R,\Omega}\right)-\frac{\partial}{{\partial
R}}U_{\nu}\left(R\right)}\right]\Phi_{\nu}\left({R;\Omega}\right)$

does not have a unique solution, $\chi_{\nu}\left({R;\Omega}\right)$. In fact,
any functions with the form of

$\chi_{\nu}\left({R;\Omega}\right)=\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R;\Omega}\right)+c\Phi_{\nu}\left({R;\Omega}\right),$ (28)

(where $c$ is an arbitrary number) can be a solution of Eq. (27). The
singularity of matrix
${H_{\rm{ad}}\left({R,\Omega}\right)-U_{\nu}\left(R\right)}$ can be removed by
considering the additional condition that

$\int{d\Omega\Phi_{\nu}\left({R,\Omega}\right)^{*}}\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R,\Omega}\right)=0,$ (29)

which can be derived from the normalization condition of
$\Phi_{\nu}\left({R,\Omega}\right)$, as shown in Ref. Nelson (1976).
Nevertheless, our numerical studies show that even without removing the
singularity, the use of numerical solution packages such as “Linear Algebra
PACKage” (LAPACK) Anderson et al. (1992) or PARDISO Naumann and Schenk (2012)
to solve Eq. (27) directly still gives an accurate solution
$\chi_{\nu}\left({R;\Omega}\right)$ in the form of Eq. (28) with an unknown
$c$. And once we have the numerical solution
$\chi_{\nu}\left({R;\Omega}\right)$, $c$ can be calculated by

$c=\int{d\Omega\Phi_{\nu}\left({R;\Omega}\right)^{*}}\chi_{\nu}\left({R;\Omega}\right).$
(30)

Finally, the derivative of $\Phi_{\nu}\left({R;\Omega}\right)$ can be written
as

$\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R;\Omega}\right)=\chi_{\nu}\left({R;\Omega}\right)-c\Phi_{\nu}\left({R;\Omega}\right),$
(31)

which can be inserted into Eq. (11) and Eq. (13) 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 $\frac{\partial}{{\partial
R}}\Phi_{\nu}\left({R;\Omega}\right)$ are numerically accurate. Therefore,
this method can give numerically very accurate coupling matrices. The major
limitation of this approach is that it becomes computationally more expensive
to obtain nonadiabatic couplings as the number of requested channels
increases. Nevertheless, for our typical calculations the number of relevant
three-body channels is not so large and this method has been shown to be quite
efficient.

## References

  * Efimov (1971) V. Efimov, Sov. J. Nucl. Phys. 12, 589 (1971).
  * Kraemer et al. (2006) T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C. Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, H.-C. Nägerl, et al., Nature 440, 315 (2006).
  * Berninger et al. (2011) M. Berninger, A. Zenesini, B. Huang, W. Harm, H. C. Nägerl, F. Ferlaino, R. Grimm, P. S. Julienne, and J. M. Hutson, Phys. Rev. Lett. 107, 120401 (2011).
  * Pollack et al. (2009) S. E. Pollack, D. Dries, and R. G. Hulet, Science 326, 1683 (2009).
  * Gross et al. (2009) N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, Phys. Rev. Lett. 103, 163202 (2009).
  * Gross et al. (2010) N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, Phys. Rev. Lett. 105, 103203 (2010).
  * Ottenstein et al. (2008) T. B. Ottenstein, T. Lompe, M. Kohnen, A. N. Wenz, and S. Jochim, Phys. Rev. Lett. 101, 203202 (2008).
  * Lompe et al. (2010) T. Lompe, T. B. Ottenstein, F. Serwane, K. Viering, A. N. Wenz, G. Zürn, and S. Jochim, Phys. Rev. Lett. 105, 103201 (2010).
  * Huckans et al. (2009) J. H. Huckans, J. R. Williams, E. L. Hazlett, R. W. Stites, and K. M. O’Hara, Phys. Rev. Lett. 102, 165302 (2009).
  * Williams et al. (2009) 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).
  * Wild et al. (2012) R. J. Wild, P. Makotyn, J. M. Pino, E. A. Cornell, and D. S. Jin, Phys. Rev. Lett. 108, 145305 (2012).
  * Barontini et al. (2009) G. Barontini, C. Weber, F. Rabatti, J. Catani, G. Thalhammer, M. Inguscio, and F. Minardi, Phys. Rev. Lett. 103, 043201 (2009).
  * Petrov (2004) D. S. Petrov, Phys. Rev. Lett. 93, 143201 (2004).
  * Wang et al. (2011a) Y. Wang, J. P. D’Incao, and B. D. Esry, Phys. Rev. A 83, 042710 (2011a).
  * Braaten and Hammer (2006) E. Braaten and H.-W. Hammer, Phys. Rep. 428, 259 (2006).
  * Hammer and Platter (2007) H. W. Hammer and L. Platter, Eur. Phys. J. A 32, 113 (2007).
  * von Stecher et al. (2009) J. von Stecher, J. P. D’Incao, and C. H. Greene, Nature Physics 5, 417 (2009).
  * Wang et al. (2011b) Y. Wang, J. P. D’Incao, , and C. H. Greene, Phys. Rev. Lett. 106, 233201 (2011b).
  * Schmidt et al. (2012) R. Schmidt, S. P. Rath, and W. Zwerger, arxiv: p. 1201.4310 (2012).
  * Chin (2011) C. Chin, arxiv: p. 1111.1484 (2011).
  * Sørensen et al. (2012) P. K. Sørensen, D. V. Fedorov, A. S. Jensen, and N. T. Ainner, arxiv: p. 1206.2274 (2012).
  * Wang et al. (2012a) J. Wang, J. P. D’Incao, B. D. Esry, and C. H. Greene, Phys. Rev. Lett. 108, 263001 (2012a).
  * Naidon et al. (2012) P. Naidon, S. Endo, and M. Ueda, arxiv: p. 1208.3912 (2012).
  * Wang et al. (2012b) Y. Wang, J. Wang, J. P. D’Incao, and C. H. Greene, arxiv: p. 1207.6439 (2012b).
  * Gao (2000) B. Gao, Phys. Rev. A 62, 050702(R) (2000).
  * Esry et al. (2001) B. D. Esry, C. H. Greene, and H. Suno, Phys. Rev. A 65, 010705 (2001).
  * Chin et al. (2011) C. Chin, R. Grimm, P. S. Julienne, and E. Tiesinga, Rev. Mod. Phys. 107, 120401 (2011).
  * Gribakin and Flambaum (1993) G. F. Gribakin and V. V. Flambaum, Phys. Rev. A 48, 546 (1993).
  * Willner and Gianturco (2006) K. Willner and F. A. Gianturco, Phys. Rev. A 74, 052715 (2006).
  * Wang et al. (2011c) J. Wang, J. P. D’Incao, and C. H. Greene, Phys. Rev. A 84, 052721 (2011c).
  * Kendrick et al. (1999) B. K. Kendrick, R. T. Pack, R. B. Walker, and E. F. Hayes, J. Chem. Phys. 110, 6673 (1999).
  * Suno et al. (2002) H. Suno, B. D. Esry, C. H. Greene, and J. P. Burke, Phys. Rev. A 65, 042725 (2002).
  * Burke (1999) J. P. Burke, Ph.D. thesis, University of Colorado at Boulder (1999).
  * Kokoouline and Greene (2003) V. Kokoouline and C. H. Greene, Phys. Rev. A 68, 012703 (2003).
  * Nelson (1976) R. B. Nelson, AIAA Journal 14, 1201 (1976).
  * Anderson et al. (1992) E. Anderson, Z. Bai, C. Bischof, J. Demmel, J. Dongarra, J. D. Cros, A. Greenbaum, S. Hammerling, A. Mckenney, S. Ostouchov, et al., _LAPACK Users Guide_ (SIAM, 1992).
  * Naumann and Schenk (2012) U. Naumann and O. Schenk, _Combinatorial Scientific Computing_ (Champman and Hall, 2012).