# Accurate simulation of Efimov physics in ultracold atomic gases with
realistic three-body multichannel interactions

J. van de Kraats<EMAIL_ADDRESS>Eindhoven University of Technology, P.
O. Box 513, 5600 MB Eindhoven, The Netherlands  S.J.J.M.F. Kokkelmans
Eindhoven University of Technology, P. O. Box 513, 5600 MB Eindhoven, The
Netherlands

###### Abstract

We give a detailed and self-contained description of a recently developed
theoretical and numerical method for the simulation of three identical bosonic
alkali-metal atoms near a Feshbach resonance, where the Efimov effect is
induced. The method is based on a direct construction of the off-shell two-
body transition matrix from exact eigenfunctions of the embedded two-body
Hamiltonians, obtained using realistic parameterizations of the interaction
potentials which accurately reproduce the molecular energy levels. The
transition matrix is then inserted into the appropriate three-body integral
equations, which may be efficiently solved on a computer. We focus especially
on the power of our method in including rigorously the effects of multichannel
physics on the three-body problem, which are usually accounted for only by
various approximations. We demonstrate the method for 7Li, where we recently
showed that a correct inclusion of this multichannel physics resolves the
long-standing disagreement between theory and experiment regarding the
Efimovian three-body parameter. We analyze the Efimovian enhancement of the
three-body recombination rate on both sides of the Feshbach resonance,
revealing strong sensitivity to the spin structure of the model thus
indicating the prevalence of three-body spin-exchange physics. Finally, we
discuss an extension of our methodology to the calculation of three-body
bound-state energies.

##  I Introduction

The Efimov effect constitutes a central result of few-body quantum physics,
describing the sudden appearance of an infinite tower of three-body bound
states in any system with short-range but resonant pairwise interactions [1,
2, 3, 4, 5, 6]. In the simplest case of zero-range contact interactions,
parameterized by an $s$-wave scattering length $a$, the associated spectrum of
Efimov states has discrete scale invariance relative to a characteristic
scaling factor $\lambda=e^{\pi/s_{0}}$, where $s_{0}=1.00624$ is a universal
dimensionless number. In configuration space this invariance originates from
an induced attractive interaction scaling precisely as $1/R^{2}$, where $R$ is
the size of the three-body system. For more realistic interactions with a
finite range, the scale invariance of the Efimov spectrum is broken through
the appearance of a characteristic length scale, now typically referred to as
the three-body parameter. Initially, the three-body parameter was expected to
be highly sensitive to short-range details of the system, thus making any
reasonably universal prediction for its value from two-body physics unlikely
unless the interaction could be characterized by a large universal effective
range [7, 8].

Interest in the Efimov effect found a resurgence once it was realized that it
could be experimentally accessed in trapped ultracold atomic gases, where
precision control of the scattering length is enabled through the use of
magnetic Feshbach resonances [9, 10, 11]. In this scheme, an Efimov state can
be experimentally characterized by scanning the scattering length on both the
repulsive and attractive side of the resonance, whilst monitoring the loss of
atoms from the trap. Then, Efimovian enhancement of three-body recombination
processes leads to a log-periodic modulation of the atomic loss rate as a
function of scattering length, at a characteristic frequency set by the Efimov
constant $s_{0}$ [12, 3, 4]. The three-body parameter, typically expressed in
terms of the characteristic scattering lengths $a_{\pm}$ and trimer widths
$\eta_{\pm}$, can then be directly extracted from the absolute phase of the
oscillations.

By exploiting Feshbach resonances, the three-body parameter has been measured
in a wide array of different atomic species amicable to laser cooling.
Surprisingly, these experiments have revealed that the three-body parameter
for broad Feshbach resonances is universally determined as $a_{-}\approx-9.7\
r_{\mathrm{vdW}}$, where $r_{\mathrm{vdW}}$ is the range of the two-body
interaction [13, 14, 15]. This van der Waals universality contradicts the
earlier expectation that the three-body parameter would be sensitive also to
the nonuniversal atomic structure, thus reigniting theoretical interest in the
Efimov effect. Detailed studies of the three-body problem subsequently
revealed that van der Waals universality originates from a characteristic
repulsive barrier in the three-body potential, which shields the particles
from probing the non-universal short range [16, 17]. Later, this framework was
extended further to show how the three-body parameter obeys distinct
universality classes dependent on the shape of the two-body interaction [18].

The studies of van der Waals universality provide a satisfactory picture of
the Efimov effect for broad Feshbach resonances, where the effective range of
the two-body interaction is small [11]. In the opposite limit of a narrow
Feshbach resonance, with large effective range, the three-body parameter is
also expected to be universally determined [8, 19], although experimental
verification in this regime is lacking due to difficulties in working with
very small magnetic field ranges. Experimental data is available however, in
the much more challenging regime of intermediate resonance widths, where no
single dominant two-body length scale can be defined. Here, theoretical models
based on multichannel two-body interactions generally predict an increase of
the value of $\absolutevalue*{a_{-}}$ as the resonance narrows, which then
smoothly connects to the analytic regime of very narrow resonances where
$a_{-}$ starts to scale directly with the effective range [20, 21, 22, 23].
Experimental measurements of $a_{-}$ obtained near intermediate Feshbach
resonances in 39K are generally in line with this prediction [24, 25, 26],
giving credence to the fundamental assumptions about the relevant microscopic
physics on which these models are based. Puzzling however, has been the case
of 7Li, where the value of $\absolutevalue{a_{-}}$ has been measured to be
smaller than the universal van der Waals value, even though the associated
Feshbach resonances are actually narrower than those used in potassium [27,
28, 29, 30]. Furthermore, most known theoretical models, even those that
attempted to include non-universal details of the two-body interaction, could
not reproduce these findings to any satisfactory degree [31, 32, 33]. This
lithium few-body puzzle, became additionally interesting when new measurements
revealed the existence of an Efimov state above the atom-dimer scattering
threshold, which was subsequently found to originate from a universal reshape
of the three-body interaction for intermediate to narrow Feshbach resonances
[34]. At roughly the same time, the anomalous value of $a_{-}$ was finally
reproduced in a numerical calculation, which included the effects of strong
three-body spin-exchange interactions, thus resolving the long-standing
discrepancy between theory and experiment [35].

Following up on the developments above, in this article we discuss in detail a
theoretical and numerical method for the three-body problem which allows for
the inclusion of all spin degrees of freedom of the interacting atoms, as was
required in resolving the lithium three-body puzzle regarding the three-body
parameter [35]. We subsequently discuss the extension of this method to also
calculate the energy of Efimov resonances, which are finding increased
experimental relevance [34], and apply the method for a study of three-body
recombination in 7Li on both sides of the Feshbach resonance, including also
the regime of positive scattering length where a shallow dimer state exists.
The paper will be structured as follows. In Sec. II we discuss in detail our
theoretical and numerical methods for obtaining both the three-body
recombination rates and the energies of Efimov states in systems of three
interacting ultracold alkali-metal atoms. Then in Sec. III we give and discuss
the results we obtain upon applying this model to the 7Li system specifically.
Finally we conclude this paper in Sec. IV, and give some outlook for further
research.

##  II Method

In this section we describe in detail our multichannel method for calculating
both the three-body recombination rate and three-body bound state energy near
a Feshbach resonance, for three bosonic alkali-metal atoms with mass $m$
labelled with Latin indices $j=1,2,3$. We work in Jacobi coordinates [36],
meaning that the three-body system is subdivided into a pair + a third
particle. We specify the choice of pair by a Greek index $\alpha$, equal to
the index $j$ of the third particle. For example if $\alpha=3$, and particles
have individual positions $\mathbf{x}_{j}$ relative to some arbitrary origin,
then the set of Jacobi coordinates is defined as,

$\displaystyle\begin{split}\mathbf{r}_{3}&=\mathbf{x}_{1}-\mathbf{x}_{2},\\\
{\bf\it\rho}_{3}&=\mathbf{x}_{3}-\frac{1}{2}(\mathbf{x}_{1}+\mathbf{x}_{2}),\\\
\mathbf{R}&=\frac{1}{3}(\mathbf{x}_{1}+\mathbf{x}_{2}+\mathbf{x}_{3}).\end{split}$
(1)

Similarly, if particles have individual momenta $\tilde{\mathbf{k}_{j}}$,

$\displaystyle\begin{split}\mathbf{k}_{3}&=\frac{1}{2}\left(\tilde{\mathbf{k}}_{1}-\tilde{\mathbf{k}}_{2}\right),\\\
\mathbf{p}_{3}&=\frac{2}{3}\tilde{\mathbf{k}}_{3}-\frac{1}{3}(\tilde{\mathbf{k}}_{1}+\tilde{\mathbf{k}}_{2}),\\\
\mathbf{K}&=\tilde{\mathbf{k}}_{1}+\tilde{\mathbf{k}}_{2}+\tilde{\mathbf{k}}_{3}.\end{split}$
(2)

We refer to $\mathbf{r}_{3},\mathbf{k}_{3}$ as the dimer separation/momentum
and ${\bf\it\rho}_{3},\mathbf{p}_{3}$ as the atom-dimer separation/momentum.
We assume translational invariance of the total three-body system, meaning
that the trapping potential is essentially flat in the regime of interest,
such that it is valid to set $\mathbf{R}=0$ or $\mathbf{K}=0$. We also note
that the particles are indistinguishable, so any choice of Jacobi set is
arbitrary, and we will see later that all equations are by nature symmetrized
with respect to different choices of $\alpha$.

###  II.1 Hamiltonian

In absence of interactions, a system of three alkali-metal atoms moving in an
external magnetic field $\mathbf{B}$ is described by the sum of single-
particle Hamiltonians [22, 37],

$\displaystyle\begin{split}\hat{H}_{0}=\sum_{j=1}^{3}\left[\hat{T}_{j}+\hat{H}_{j}^{\mathrm{hf}}+\hat{H}_{j}^{\mathrm{Z}}(\mathbf{B})\right].\end{split}$
(3)

Here $\hat{T}_{j}$ is the kinetic energy operator for particle $j$, and
$\hat{H}_{j}^{\mathrm{hf}},\hat{H}_{j}^{\mathrm{Z}}(\mathbf{B})$ are the
single-particle hyperfine and Zeeman Hamiltonians which account for the spin
state of the particle. They are defined explicitly as,

$\displaystyle\begin{split}\hat{H}_{j}^{\mathrm{hf}}=A_{j}^{\mathrm{hf}}\hat{\mathbf{s}}_{j}\cdot\hat{\mathbf{i}}_{j},\qquad\hat{H}_{j}^{\mathrm{Z}}(\mathbf{B})=\left(\gamma_{j}^{n}\hat{\mathbf{i}}_{j}+\gamma_{j}^{e}\hat{\mathbf{s}}_{j}\right)\cdot\mathbf{B},\end{split}$
(4)

where $\hat{\mathbf{s}}_{j}$ and $\hat{\mathbf{i}}_{j}$ are the electronic and
nuclear spin operators, $A_{j}^{\mathrm{hf}}$ is the hyperfine energy, and
$\gamma_{j}^{\mathrm{n/e}}$ are the nuclear and electronic gyromagnetic
ratios. Upon diagonalizing
$\hat{H}_{j}^{\mathrm{hf}}+\hat{H}_{j}^{\mathrm{Z}}$ one obtains a spectrum of
single particle eigenstates with magnetic field dependent energy
$\varepsilon_{c_{j}}$, as shown in Fig. 1(a) for the case of 7Li. By
convention, these states are labelled in alphabetical order with respect to
increasing energy, or alternatively by the hyperfine state
$\ket*{f_{j},m_{f_{j}}}$ to which they adiabatically connect when
$\absolutevalue{\mathbf{B}}\rightarrow 0$. Here
$\hat{\mathbf{f}}_{j}=\hat{\mathbf{s}}_{j}+\hat{\mathbf{i}}_{j}$, and
$m_{f_{j}}$ is the associated projection on the magnetic field direction.

Figure 1:  Characteristics of the coupled-channels few-body problem [38]. In
(a) we show the spectrum $\varepsilon_{c_{j}}$ of channel states for a single
7Li atom, using the labelling convention as explained in the main text. States
highlighted in red correspond to the initial states studied in this work,
which possess Feshbach resonances at the magnetic field values indicated by
the vertical dashed lines. In (b) we plot the two-body singlet and triplet
potential energy surfaces $V_{\alpha}^{S}(r_{\alpha})$, as formulated for 7Li
in Ref. [39].

The total Hamiltonian including interactions is written as,

$\displaystyle\begin{split}\hat{H}=\hat{H}_{0}+\sum_{\alpha}\hat{V}_{\alpha},\end{split}$
(5)

where $\hat{V}_{\alpha}$ represents the two-body interaction for the pair
$\alpha$. For neutral atoms at large interparticle separations
($r_{\alpha}\rightarrow\infty$), the dominant contribution to the interaction
is given by the attractive van der Waals interaction, i.e.
$V_{\alpha}(r_{\alpha})\rightarrow-C_{6}/r_{\alpha}^{6}$. Here the dispersive
coefficient $C_{6}$ defines the characteristic van der Waals length scale
$r_{\mathrm{vdW}}=\frac{1}{2}\left(\frac{mC_{6}}{\hbar^{2}}\right)^{\frac{1}{4}}\approx
32.4863a_{0}$. At shorter separations, where $r_{\alpha}\ll r_{\mathrm{vdW}}$,
the overlapping wave functions of the valence electrons are projected into
either singlet ($S=0$) or triplet ($S=1$) spin configurations. In this regime,
the interaction accordingly decomposes into two components,

$\displaystyle\begin{split}\hat{V}_{\alpha}=\hat{V}_{\alpha}^{0}\hat{\mathcal{P}}_{\alpha}^{0}+\hat{V}_{\alpha}^{1}\hat{\mathcal{P}}_{\alpha}^{1},\end{split}$
(6)

where $\hat{\mathcal{P}}_{\alpha}^{S}$ projects the pair into the electronic
spin state $S$. The detailed shape of the singlet and triplet molecular
potentials $\hat{V_{\alpha}}^{S}$ is generally very complicated, depending on
non-universal detail of the atomic structure. We use a realistic
parameterization of the molecular potential energy surface obtained by fitting
two-body coupled-channels calculations to low-energy scattering and
spectroscopic data [39], as plotted in Fig. 1(b).

###  II.2 Three-body basis

In a uniform system, the state of the paired particles can be represented by
the plane wave state $\ket*{\mathbf{k}_{\alpha}}$ [40], where,

$\innerproduct*{\mathbf{r}_{\alpha}}{\mathbf{k}_{\alpha}}=\frac{1}{(2\pi)^{\frac{3}{2}}}e^{i\mathbf{k}_{\alpha}\cdot\mathbf{r}_{\alpha}},\qquad\int
d^{3}\mathbf{k}_{\alpha}\outerproduct*{\mathbf{k}_{\alpha}}{\mathbf{k}_{\alpha}}=\mathbb{I}$
(7)

These basis states are not very efficient for the present computations since
they do not exploit the spherical symmetry of the pairwise interaction.
Instead, we work in a partial wave basis
$\ket*{k_{\alpha},l_{\alpha},m_{l_{\alpha}}}$, where $l_{\alpha}$ is the
orbital angular momentum quantum number of the dimer subsystem, and
$m_{l_{\alpha}}$ the associated projection on the quantization axis. The
normalization reads [41],

$\displaystyle\begin{split}&\innerproduct*{k_{\alpha}^{\prime},l_{\alpha}^{\prime},m_{l_{\alpha}}^{\prime}}{k_{\alpha},l_{\alpha},m_{l_{\alpha}}}=\frac{1}{4\pi}\frac{\delta(k_{\alpha}^{\prime}-k_{\alpha})}{k_{\alpha}^{\prime}k_{\alpha}}\delta_{l_{\alpha}l_{\alpha}^{\prime}}\delta_{m_{l_{\alpha}}m_{l_{\alpha}}^{\prime}},\\\
&\sum_{l_{\alpha},m_{l_{\alpha}}}\int_{0}^{\infty}dk\ 4\pi
k_{\alpha}^{2}\ket{k_{\alpha},l_{\alpha},m_{l_{\alpha}}}\bra{k_{\alpha},l_{\alpha},m_{l_{\alpha}}}=\mathbb{I}.\end{split}$
(8)

As follows from the Wigner-Eckart theorem, any spherically symmetric operator
in this basis is independent of $m_{l_{\alpha}}$ and diagonal in $l_{\alpha}$.
The position and momentum projections are,

$\displaystyle\begin{split}\innerproduct*{\mathbf{r}_{\alpha}}{k_{\alpha},l_{\alpha},m_{l_{\alpha}}}&=\frac{i^{l_{\alpha}}}{\sqrt{2}\pi}\frac{\hat{j}_{l_{\alpha}}(k_{\alpha}r_{\alpha})}{k_{\alpha}r_{\alpha}}Y_{l_{\alpha}}^{m_{l_{\alpha}}}(\hat{\mathbf{r}}_{\alpha}),\\\
\innerproduct*{\mathbf{k}_{\alpha}^{\prime}}{k_{\alpha},l_{\alpha},m_{l_{\alpha}}}&=\frac{1}{2\sqrt{\pi}}\frac{\delta(k_{\alpha}-k_{\alpha}^{\prime})}{k_{\alpha}k_{\alpha}^{\prime}}Y_{l_{\alpha}}^{m_{l_{\alpha}}}(\hat{\mathbf{k}}_{\alpha}^{\prime}),\end{split}$
(9)

where $\hat{j}_{l_{\alpha}}(k_{\alpha}r_{\alpha})$ is the Riccati-Bessel
function [40], and
$Y_{l_{\alpha}}^{m_{l_{\alpha}}}(\hat{\mathbf{k}}_{\alpha})$ a spherical
harmonic.

The extension to the three-body sector follows naturally by defining an
equivalent basis state
$\ket*{p_{\alpha},\lambda_{\alpha},m_{\lambda_{\alpha}}}$ for the atom-dimer
degree of freedom. As before, $p_{\alpha}$ is the atom-dimer momentum,
$\lambda_{\alpha}$ the atom-dimer angular momentum, and $m_{\lambda_{\alpha}}$
its projection. Total rotational symmetry of the system then motivates
coupling the dimer and atom-dimer angular momenta to obtain the total angular
momentum $\mathbf{L}=\mathbf{l}_{\alpha}+{\bf\it\lambda}_{\alpha}$ [41]. The
associated basis states are formed as,

$\displaystyle\begin{split}\ket*{l_{\alpha},\lambda_{\alpha},L,M_{L}}=\sum_{m_{l_{\alpha}}m_{\lambda_{\alpha}}}\innerproduct*{l_{\alpha},m_{l_{\alpha}},\lambda_{\alpha},m_{\lambda_{\alpha}}}{L,M_{L}}\ket*{l_{\alpha},m_{l_{\alpha}},\lambda_{\alpha},m_{\lambda_{\alpha}}}.\end{split}$
(10)

Here
$\innerproduct*{l_{\alpha},m_{l_{\alpha}},\lambda_{\alpha},m_{\lambda_{\alpha}}}{L,M_{L}}$
is a Clebsch-Gordan coefficient. Spherical symmetry dictates that $L,M_{L}$
are both good quantum numbers. Including also the absolute momenta, we now
write the total three-body partial wave basis as
$\ket*{k_{\alpha}p_{\alpha}(l_{\alpha}\lambda_{\alpha})LM_{L}}$ [41]. It is an
eigenstate of the kinetic energy operator,

$\displaystyle\begin{split}\sum_{j}T_{j}\ket*{k_{\alpha}p_{\alpha}(l_{\alpha}\lambda_{\alpha})LM_{L}}=\left(\frac{\hbar^{2}k_{\alpha}^{2}}{m}+\frac{3}{4}\frac{\hbar^{2}p_{\alpha}^{2}}{m}\right)\ket*{k_{\alpha}p_{\alpha}(l_{\alpha}\lambda_{\alpha})LM_{L}}.\end{split}$
(11)

Going forward we assume $L=M_{L}=0$, as this term gives the dominant
contribution to the three-body recombination rate at low energies [42]. This
immediately enforces that $\lambda_{\alpha}=l_{\alpha}$, so we only need a
single orbital quantum number $l_{\alpha}$ to specify the three-body state.

To describe the spin state of particles, we do not work directly with the
hyperfine states $\ket*{f_{j},m_{f_{j}}}$ as defined in the previous section,
but instead use the basis of atomic or Zeeman single-particle states
$\ket*{\xi_{j}}=\ket*{m_{s_{j}},m_{i_{j}}}$, where $m_{s_{j}},m_{i_{j}}$ are
the magnetic projections of the electronic and nuclear spin respectively. In
the three-body sector, the only rigorously conserved quantum number is the
total magnetic spin projection $M_{F}$,

$\displaystyle
M_{F}=\sum_{j=1}^{3}m_{s_{j}}+\sum_{j=1}^{3}m_{i_{j}}\equiv\sum_{j=1}^{3}m_{f_{j}}.$
(12)

This means that, for any particular choice of Jacobi set $\alpha$, the spin
projection for the dimer particles equals $M_{F}-m_{f_{\alpha}}$, which is
known to be preserved by $V_{\alpha}$. Hence the spin state of the third
particle, $\ket*{\xi_{\alpha}}$, will determine the set of coupled spin states
for the dimer subsystem, henceforth notated as $\ket*{\Xi_{\alpha}}$. For
identical particles, the dimension of this set can be reduced by
symmetrization,

$\displaystyle\begin{split}\ket*{\bar{\Xi}_{\alpha}}=\frac{1+(-1)^{l}\hat{P}_{\alpha}}{\sqrt{2\left(1+\delta_{\alpha}\right)}}\ket*{\Xi_{\alpha}},\end{split}$
(13)

Here $\hat{P}_{\alpha}$ denotes a permutation operator that exchanges the two
paired particles, and $\delta_{\alpha}$ equals unity when the two particles
are in identical states and zero otherwise. The factor $(-1)^{l}$ corrects for
the spatial parity due to the angular momentum. Now, in the three-body sector,
we form the matrix
$\sum_{j=1}^{3}\matrixelement*{\bar{\Xi}_{\alpha}^{\prime},\xi_{\alpha}^{\prime}}{\hat{H}_{j}^{\mathrm{hf}}+\hat{H}_{j}^{\mathrm{Z}}}{\bar{\Xi}_{\alpha},\xi_{\alpha}}$,
which is diagonalized to obtain eigenstates
$\ket*{\bar{C}_{\alpha},c_{\alpha}}$, referred to as channel states, with
energies
$E_{\bar{C}_{\alpha},c_{\alpha}}=\varepsilon_{\bar{C}_{\alpha}}+\varepsilon_{c_{\alpha}}$.
In Fig. 2 we illustrate this three-body channel spectrum for 7Li, given two
distinct values of the total projection $M_{F}$. Here one observes clearly how
$M_{F}$ sets the total number of coupled spin states, and accordingly
influences the complexity of the three-body calculation. For the remainder of
this section, we will set $E=0$ at the three-body scattering threshold, given
by the lowest three-body channel energy as highlighted in red in Fig. 2.

Figure 2:  Three-body channel energies $E_{\bar{C},c}$ for 7Li as a function
of magnetic field, for two values of the total magnetic quantum number $M_{F}$
relevant to this work. The state with lowest energy, setting the three-body
scattering threshold, is drawn in red. Vertical dashed lines show the location
of the Feshbach resonance studied in this work.

This concludes the definition of the three-body basis,
$\ket*{k_{\alpha}p_{\alpha}l_{\alpha};\bar{C}_{\alpha}c_{\alpha}}\equiv\ket*{k_{\alpha}p_{\alpha}(l_{\alpha}l_{\alpha})00}\otimes\ket*{\bar{C}_{\alpha}c_{\alpha}}$.
For notational brevity, we will suppress the Jacobi index $\alpha$ on momenta
and quantum numbers going forward, as the appropriate Jacobi index can always
be inferred from the surrounding equation.

###  II.3 Off-shell two-body transition matrix

As input to the three-body calculation, we require knowledge of the two-body
transition operator for a given partial wave $l$ and third particle state $c$,
denoted as $\hat{t}_{\alpha}^{lc}(z)$. Here $z$ is the two-body energy, which
will have to be evaluated off-shell as the two-body subsystem may lose or gain
energy by interacting with the third particle [36]. In general, the transition
operator can be defined through the Lippmann-Schwinger equation, which has the
formal solution [40],

$\displaystyle\hat{t}_{\alpha}^{lc}(z)=\hat{V}_{\alpha}+\hat{V}_{\alpha}\hat{g}_{\alpha}^{lc}(z)\hat{V}_{\alpha},$
(14)

with $\hat{g}_{\alpha}^{lc}(z)=(z-\hat{h}_{\alpha}^{lc})^{-1}$ the two-body
Green’s operator, and two-body Hamiltonian $\hat{h}_{\alpha}^{lc}$. Defining
two-body eigenstates as
$\hat{h}_{\alpha}^{lc}\ket*{\psi_{lc}^{\nu}}=\varepsilon_{lc}^{\nu}\ket*{\psi_{lc}^{\nu}}$,
with state index $\nu$, the Green’s operator can be expanded as,

$\displaystyle\begin{split}\hat{g}_{\alpha}^{lc}(z)=\sum_{\nu}\frac{\ket*{\psi_{lc}^{\nu}}\bra*{\psi_{lc}^{\nu}}}{z-\varepsilon_{lc}^{\nu}},\end{split}$
(15)

Inserting this expansion into Eq. (14) and evaluating in the two-body partial
wave and spin basis gives,

$\displaystyle\matrixelement*{k^{\prime}l;\bar{C}^{\prime}}{\hat{t}_{\alpha}^{lc}(z)}{kl;\bar{C}}$
$\displaystyle=\matrixelement*{k^{\prime}l;\bar{C}^{\prime}}{\hat{V}_{\alpha}}{kl;\bar{C}}+\sum_{\nu}\matrixelement*{k^{\prime}l;\bar{C}^{\prime}}{\hat{V}_{\alpha}\frac{\ket*{\psi_{lc}^{{\nu}}}\bra*{\psi_{lc}^{{\nu}}}}{z-\varepsilon_{lc}^{\nu}}\hat{V}_{\alpha}}{kl;\bar{C}}.$
(16)

Here we have integrated away the angular degrees of freedom, such that the
inconsequential quantum number $m_{l}$ drops out. We will now use Eq. (16) to
construct the transition matrix from the two-body eigenfunctions
$\ket*{\psi_{lc}^{{\nu}}}$, which we obtain by solving the associated quantum
mechanical scattering problem in position space. First, we exploit spherical
symmetry of $\hat{V}_{\alpha}$ by expanding the wave function as,

$\displaystyle\innerproduct*{r}{\psi_{lc}^{{\nu}}}=\sum_{\bar{C}}\frac{u_{lc\bar{C}}^{\nu}(r)}{r}\ket*{\bar{C}},$
(17)

where $u_{lc\bar{C}}^{\nu}(r)$ is the radial wave function for each individual
two-body scattering channel. It obeys the radial Schrödinger equation,

$\displaystyle\begin{split}\left[\frac{\hbar^{2}}{m}\left(-\derivative[2]{r}+\frac{l(l+1)}{r^{2}}\right)+V_{\alpha}(r)\right]u_{lc\bar{C}}^{\nu}(r)=\left(\varepsilon_{lc}^{\nu}-\varepsilon_{\bar{C}}\right)u_{lc\bar{C}}^{\nu}(r).\end{split}$
(18)

We solve Eq. (18) via a discrete variable representation (DVR) [43], with a
transformed grid in $r$ for increased efficiency [44, 45]. Although such
methods are typically used to calculate molecular binding energies,
constructing the transition matrix requires not just the bound states but also
the scattering states, for which
$\varepsilon_{lc}^{\nu}-\varepsilon_{\bar{C}}>0$. In the DVR these states are
not exact, due to the discretization of the coordinate space via a hard-box
boundary condition at some (large) two-body separation $r_{b}$, which should
far exceed all other length scales. For a detailed study of the effects of
this box on the resulting transition matrix and three-body observables, we
refer to ref. [46]. Based on the findings of that reference, we will use
$r_{b}=2000r_{\mathrm{vdW}}$ in all our numerical calculations for 7Li.

Upon obtaining $u_{lc\bar{C}}^{\nu}(r)$ and $\varepsilon_{lc}^{\nu}$ we
calculate the matrix elements
$\matrixelement*{k^{\prime}l;\bar{C}^{\prime}}{\hat{t}_{\alpha}^{lc}(z)}{kl;\bar{C}}$
directly from Eq. (16). The momentum grid $k$ is composed using $N_{k}$ zeros
of the spherical Bessel function, such that $j_{l}(kr_{b})=0$, up to a maximum
two-body momentum $k_{\mathrm{max}}$. Once the transition matrix is obtained,
we diagonalize to obtain the spectral expansion,

$\displaystyle\begin{split}\matrixelement*{k^{\prime}l;\bar{C}^{\prime}}{\hat{t}_{\alpha}^{lc}(z)}{kl;\bar{C}}=\sum_{n}\innerproduct*{k^{\prime};\bar{C}}{\chi_{lc}^{n}(z)}\tau_{lc}^{n}(z)\innerproduct*{\chi_{lc}^{n}(z)}{k;\bar{C}},\end{split}$
(19)

with eigenvalues $\tau_{lc}^{n}(z)$ and form factors
$\innerproduct*{k;\bar{C}}{\chi_{lc}^{n}(z)}$. Once more we suppress the
Jacobi index $\alpha$ on the right-hand side. Note that we represent
eigenfunctions of the two-body Hamiltonian with index $\nu$ and eigenfunctions
of the transition matrix with index $n$, and will maintain this convention
going forward.

The number of terms $n$ for each $l$, $c$ and a given value of the energy $z$
equals the number of grid points in $k$ times the number of two-body channels.
As this number is typically too large to store the full transition matrix in
memory, we sort the spectral expansion (19) with respect to decreasing
magnitude of $\tau_{lc}^{n}(z)$, and subsequently truncate the expansion at a
term $n_{\mathrm{cut}}$, which sets the precision to which we include the
transition matrix. In that case, Eq. (19) represents a separable approximation
to the formal transition matrix, which will allow us to write the three-body
equations in the next sections into a one-dimensional form that can be
efficiently solved on a computer [36]. For further discussion of different
separable approximations to the transition matrix and their consequences for
studying Efimov physics, we refer the interested reader to Ref. [47].

###  II.4 Three-body recombination rate

We are interested in the transition rate from a free zero-energy three-body
state $\ket*{\Psi_{\mathrm{in}}}$ to an outgoing atom-dimer state
$\ket*{(\mathbf{p},c),\varphi_{\nu}}$. The latter describes a state in which
the pair is in the molecular state $\ket*{\varphi_{\nu}}$, while the third
particle is in a plane wave state $\ket*{\mathbf{p};c}$ with atom-dimer
momentum $\mathbf{p}$ and spin state $c$. Such a process induces a particle
density loss $\mathrm{d}n/\mathrm{d}t=-K_{3}\ n^{3}$, where the recombination
rate constant $K_{3}$ is determined as [48, 49, 50],

$\displaystyle K_{3}=\frac{12\pi m}{\hbar}(2\pi\hbar)^{6}\lim_{E\rightarrow
0}\sum_{{\nu},c}\int d\hat{\mathbf{p}}\
p\absolutevalue*{\matrixelement*{(\mathbf{p},c),\varphi_{\nu}}{\hat{U}_{\alpha
0}(E)}{\Psi_{\mathrm{in}}}}^{2}.$ (20)

Here $\hat{U}_{\alpha 0}(E)$ is the appropriate three-body transition
operator. By inserting a complete set of three-body partial-wave basis states,
one finds,

$\displaystyle\begin{split}\int
d\hat{\mathbf{p}}\absolutevalue*{\matrixelement*{(\mathbf{p},c),\varphi_{\nu}}{\hat{U}_{\alpha
0}(E)}{\Psi_{\mathrm{in}}}}^{2}&=4\pi\sum_{l}\absolutevalue*{\matrixelement*{\varphi_{l}^{\nu}pl;c}{\hat{U}_{\alpha
0}(E)}{kp0;\bar{C}_{\mathrm{in}}c_{\mathrm{in}}}}^{2}.\end{split}$ (21)

Here we have defined explicitly the $s$-wave incoming state
$\ket*{\Psi_{\mathrm{in}}}=\lim_{E\rightarrow
0}\ket*{kp0;,\bar{C}_{\mathrm{in}}c_{\mathrm{in}}}$, and have decomposed the
bound state $\ket*{\varphi_{\nu}}$ into partial-wave components
$\ket*{\varphi_{l}^{\nu}}$ which capture the state of the pair. Accordingly,
the full three-body state can then be written as
$\ket*{\varphi_{l}^{\nu}pl;c}$.

To evaluate the matrix elements of $\hat{U}_{\alpha 0}(E)$, we solve the
appropriate Alt-Grassberger-Sandhas (AGS) equation for three identical
particles [51, 41],

$\displaystyle\hat{U}_{\alpha
0}(E)=\hat{G}_{0}^{-1}(E)+\hat{P}\hat{\mathcal{T}}_{\alpha}(E)\hat{G}_{0}(E)\hat{U}_{\alpha
0}(E).$ (22)

Here $\hat{G}_{0}(E)=(E-\hat{H}_{0})^{-1}$ is the uncoupled three-body Green’s
function, and $\hat{P}=\hat{P}_{+}+\hat{P}_{-}$ is the sum of cyclic and
anticyclic permutation operators of the three particle indices. The two-body
problem of the previous section enters through the generalized two-body
transition matrix, defined as,

$\displaystyle\hat{\mathcal{T}}_{\alpha}(E)=\hat{V}_{\alpha}+\hat{V}_{\alpha}\hat{G}_{\alpha}(E)\hat{V}_{\alpha},$
(23)

where $\hat{G}_{\alpha}(E)=(E-\hat{H}_{0}-\hat{V}_{\alpha})^{-1}$. Note that
by this definition $\hat{\mathcal{T}}_{\alpha}(E)$ is diagonal in the
coordinates of the third particle, and its action on the pair is equivalent to
the two-body transition matrix. The AGS equation can be rewritten to,

$\displaystyle\hat{U}_{\alpha
0}(E)=\frac{1}{3}\hat{G}_{0}^{-1}(E)(1+\hat{P})+\hat{P}\hat{\mathcal{T}}_{\alpha}(E)\hat{G}_{0}(E)\hat{U}_{\alpha
0}(E).$ (24)

Since the first term vanishes on the energy shell, we may as well define a new
operator $\hat{\bar{U}}_{\alpha
0}(E)=\hat{P}\hat{\mathcal{T}}_{\alpha}(E)\hat{G}_{0}(E)\hat{U}_{\alpha
0}(E)$, which then obeys,

$\displaystyle\hat{\bar{U}}_{\alpha
0}(E)=\frac{1}{3}\hat{P}\hat{\mathcal{T}}_{\alpha}(E)(1+\hat{P})+\hat{P}\hat{\mathcal{T}}_{\alpha}(E)\hat{G}_{0}(E)\hat{\bar{U}}_{\alpha
0}(E).$ (25)

Finally, for identical particles it is valid to write $\hat{P}=2\hat{P}_{+}$,
such that,

$\displaystyle\begin{split}\hat{\bar{U}}_{\alpha
0}(E)=\frac{2}{3}\hat{P}_{+}\hat{T}_{\alpha}(E)\left(1+2\hat{P}_{+}\right)+2\hat{P}_{+}\hat{\mathcal{T}}_{\alpha}(E)\hat{G}_{0}(E)\hat{\bar{U}}_{\alpha
0}(E).\end{split}$ (26)

This equation may now be used to evaluate the matrix element
$\matrixelement*{\varphi_{l}^{\nu}pl;c}{\hat{U}_{\alpha
0}(E)}{kp0;\bar{C}_{\mathrm{in}}c_{\mathrm{in}}}$, by inserting complete sets
of three-body basis states $\ket*{kpl;\bar{C}c}$ between all operators. The
required matrix elements of the transition operator read,

$\displaystyle\begin{split}&\matrixelement*{kpl;\bar{C}c}{\hat{\mathcal{T}}_{\alpha}(E)}{k^{\prime}p^{\prime}l^{\prime};\bar{C}^{\prime}c^{\prime}}=\frac{1}{4\pi}\frac{\delta(p-p^{\prime})}{pp^{\prime}}\delta_{ll^{\prime}}\delta_{cc^{\prime}}\sum_{n}\innerproduct*{k;\bar{C}}{\chi_{l}^{n}\left(Z_{c}(p)\right)}\tau_{l}^{n}(Z_{c}(p))\innerproduct*{\chi_{l}^{n}\left(Z_{c}(p)\right)}{k^{\prime};\bar{C}^{\prime}},\end{split}$
(27)

where we have inserted the spectral expansion in Eq. (19). We have also
defined,

$\displaystyle\begin{split}Z_{c}(p)=E-\frac{3}{4}\frac{\hbar^{2}p^{2}}{m}-\varepsilon_{c},\end{split}$
(28)

which represents the non-interacting energy of the pair. For brevity we
suppress the third particle index $c$ for $\tau$ and $\chi$, as it can be
inferred from the energy argument $Z_{c}$. To write the AGS equation we also
need the matrix elements of $\hat{P}_{+}$, referred to as the recoupling
coefficients [41]. We follow Ref. [22] and decompose the operator as
$\hat{P}_{+}=\hat{P}_{+}^{s}\hat{P}_{+}^{c}$, where $\hat{P}_{+}^{s}$ acts
only on the spin indices and $\hat{P}_{+}^{c}$ acts only on coordinates. The
matrix elements of $\hat{P}_{+}^{c}$ are well known in literature and read
[41],

$\displaystyle\begin{split}&\matrixelement*{kpl}{\hat{P}_{+}^{c}}{k^{\prime}p^{\prime}l^{\prime}}=\frac{1}{32\pi^{2}}(-1)^{l^{\prime}}\sqrt{(2l+1)(2l^{\prime}+1)}\int_{-1}^{1}dx\frac{\delta\left(\sqrt{p^{\prime
2}+\frac{1}{4}p^{2}+pp^{\prime}x}-k\right)}{k\sqrt{p^{\prime
2}+\frac{1}{4}p^{2}+pp^{\prime}x}}\frac{\delta\left(\sqrt{p^{2}+\frac{1}{4}p^{\prime
2}+pp^{\prime}x}-k^{\prime}\right)}{k^{\prime}\sqrt{p^{2}+\frac{1}{4}p^{\prime
2}+pp^{\prime}x}}\\\ &\hskip 227.62204pt\times
P_{l}\left[\frac{\frac{1}{2}p^{2}+pp^{\prime}x}{p\sqrt{p^{\prime
2}+\frac{1}{4}p^{2}+pp^{\prime}x}}\right]P_{l^{\prime}}\left[\frac{\frac{1}{2}p^{\prime
2}+pp^{\prime}x}{p^{\prime}\sqrt{p^{2}+\frac{1}{4}p^{\prime
2}+pp^{\prime}x}}\right],\end{split}$ (29)

where $P_{l}(x)$ is the Legendre polynomial. Upon inserting into the AGS
equation and integrating away the $\delta$-functions, one finds that the
recoupling of angular momenta can be conveniently expressed by the function,

$\displaystyle\begin{split}\mathfrak{P}_{ll^{\prime}}(p,q,x)=(-1)^{l^{\prime}}\sqrt{(2l+1)(2l^{\prime}+1)}P_{l}\left[\frac{\frac{1}{2}p^{2}+pqx}{p\pi(q,p,x)}\right]P_{l^{\prime}}\left[\frac{\frac{1}{2}q^{2}+pqx}{q\pi(p,q,x)}\right],\end{split}$
(30)

where we have adopted the following notation from Ref. [41],

$\displaystyle\begin{split}\pi(p,q,x)=\sqrt{p^{2}+\frac{1}{4}q^{2}+pqx}.\end{split}$
(31)

The spin permutation operator may be calculated following the method of Ref.
[22]. For example, if we choose $\alpha=3$ it is given as,

$\displaystyle\begin{split}\matrixelement*{\bar{C}_{3}c_{3};l}{\hat{P}_{+}^{\mathrm{s}}}{\bar{C}_{3}^{\prime}c_{3}^{\prime};l^{\prime}}&=\frac{1}{\sqrt{2(1+\delta_{c_{1}c_{2}})}}\frac{1}{\sqrt{2(1+\delta_{c_{1}^{\prime}c_{2}^{\prime}})}}\bigg{[}\delta_{c_{1}c_{2}^{\prime}}\delta_{c_{2}c_{3}^{\prime}}\delta_{c_{3}c_{1}^{\prime}}+(-1)^{l}\delta_{c_{2}c_{2}^{\prime}}\delta_{c_{1}c_{3}^{\prime}}\delta_{c_{3}c_{1}^{\prime}}\\\
&\hskip
142.26378pt+(-1)^{l^{\prime}}\delta_{c_{1}c_{1}^{\prime}}\delta_{c_{2}c_{3}^{\prime}}\delta_{c_{3}c_{2}^{\prime}}+(-1)^{l+l^{\prime}}\delta_{c_{2}c_{1}^{\prime}}\delta_{c_{1}c_{3}^{\prime}}\delta_{c_{3}c_{2}^{\prime}}\bigg{]}\end{split}$
(32)

Here it is useful to note that for single-channel interactions, where
particles are always in identical spin states, $\hat{P}_{+}^{s}$ forbids
coupling to dimer states with $l$ odd. Finally, to write the AGS equation as a
linear set of algebraic equations we realize that the atom-dimer state can be
written explicitly in terms of the eigenfunctions of the two-body transition
matrix, following from the spectral expansion of Eq. (23),

$\displaystyle\begin{split}\hat{\mathcal{T}}_{\alpha}(z)&\approx\frac{\hat{V}_{\alpha}\ket*{\varphi_{l}^{\nu}pl;c}\bra*{\varphi_{l}^{\nu}pl;c}\hat{V}_{\alpha}}{z-\varepsilon_{\nu}^{l}-\frac{3}{4}\frac{\hbar^{2}p^{2}}{m}-\varepsilon_{c}}\qquad\mathrm{for}\qquad
z\approx\varepsilon_{l}^{\nu}+\frac{3}{4}\frac{\hbar^{2}p^{2}}{m}+\varepsilon_{c}.\end{split}$
(33)

It follows that we can write,

$\displaystyle\begin{split}\ket*{\varphi_{l}^{\nu}pl;c}=\hat{G}_{0}(E)\ket*{\chi_{l}^{n}\left(Z_{c}(p)\right)}\ket*{pl;c},\end{split}$
(34)

connecting each eigenstate $\nu$ of the two-body Hamiltonian with a unique
term $n$ in the spectral expansion of the two-body transition matrix.

We are now finally in a position to write the AGS equation in the partial-wave
representation. For notational convenience we will write the recombination
matrix element as,

$\displaystyle\begin{split}\bar{U}_{lc}^{n}(p)=\matrixelement*{\varphi_{l}^{n}pl;c}{\hat{\bar{U}}_{\alpha
0}}{000;\bar{C}_{\mathrm{in}}c_{\mathrm{in}}},\end{split}$ (35)

which is stored numerically as a vector. The full AGS equation can then be
written as,

$\displaystyle\begin{split}\bar{U}_{lc}^{n}(p)&=\frac{2}{3}\sum_{n^{\prime}}\sum_{c^{\prime}}\mathcal{Z}_{l0;cc^{\prime}}^{nn^{\prime}}\left(p,0\right)\tau_{0}^{n^{\prime}}\left(Z_{c^{\prime}}(0)\right)\sum_{\bar{C}}\innerproduct*{\chi_{0}^{n^{\prime}}\left(Z_{c^{\prime}}(0)\right)}{0;\bar{C}}\matrixelement*{\bar{C}c^{\prime}}{1+2\hat{P}_{+}^{\mathrm{s}}}{\bar{C}_{\mathrm{in}}c_{\mathrm{in}}}\\\
&+8\pi\sum_{l^{\prime}}\sum_{n^{\prime}}\sum_{c^{\prime}}\int_{0}^{\infty}dq\
q^{2}\
\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,q\right)\tau_{l^{\prime}}^{n^{\prime}}\left(Z_{c^{\prime}}(q)\right)\bar{U}_{c^{\prime}}^{n^{\prime}l^{\prime}}(q).\end{split}$
(36)

Here we have defined the kernel function as,

$\displaystyle\begin{split}&\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,q\right)=\frac{1}{2}\sum_{\bar{C}\bar{C}^{\prime}}\int_{-1}^{1}dx\
\frac{\innerproduct*{\chi_{l}^{n}\left(Z_{c}(p)\right)}{\chi_{l}^{n}\left(Z_{c}(p)\right)}{\pi(q,p,x);\bar{C}}\innerproduct*{\pi(p,q,x);\bar{C}^{\prime}}{\chi_{l^{\prime}}^{n^{\prime}}\left(Z_{c^{\prime}}(q)\right)}}{E-\frac{\hbar^{2}}{m}\left(q^{2}+p^{2}+qpx\right)-E_{\bar{C},c}}\mathfrak{P}_{ll^{\prime}}(p,q,x)\matrixelement*{\bar{C}c}{P_{+}^{\mathrm{s}}}{\bar{C}^{\prime}c^{\prime}}.\end{split}$
(37)

In these equations $q$ is an integrated atom-dimer momentum, and
$x=\hat{\mathbf{p}}\cdot\hat{\mathbf{q}}$.

###  II.5 Three-body bound states

The energy of three-body bound states can be obtained from a similar but
homogeneous integral equation, commonly referred to as the Faddeev equation
[52, 41, 36]. It reads,

$\displaystyle\begin{split}\ket*{\Phi_{\alpha}(E)}=2\hat{G}_{0}(E)\hat{\mathcal{T}}_{\alpha}(E)\hat{P}_{+}\ket*{\Phi_{\alpha}(E)},\end{split}$
(38)

for a bound state $\ket*{\Phi(E)}=(1+\hat{P})\ket*{\Phi_{\alpha}(E)}$ at
three-body energy $E<0$. Just like for the AGS equation, the Faddeev equation
can be expressed in the partial-wave three-body basis. For a separable
transition matrix it can then be reduced to a one-dimensional integral
equation by splitting off the dimer degrees of freedom as,

$\displaystyle\begin{split}\innerproduct*{kpl;\bar{C}c}{\Phi_{\alpha}(E)}=\sum_{n}\frac{\tau_{l}^{n}\left(Z_{c}(p)\right)\innerproduct*{k;\bar{C}}{\chi_{l}^{n}(Z_{c}(p;E))}\Phi_{c}^{nl}(p;E)}{E-\frac{\hbar^{2}k^{2}}{m}-\frac{3}{4}\frac{\hbar^{2}p^{2}}{m}-E_{\bar{C},c}}.\end{split}$
(39)

The Faddeev equation becomes,

$\displaystyle\begin{split}&\Phi_{lc}^{n}(p;E)=8\pi\sum_{l^{\prime}}\sum_{c^{\prime}}\sum_{n^{\prime}}\int_{0}^{\infty}dq\
q^{2}\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,q\right)\tau_{l^{\prime}}^{n^{\prime}}\left(Z_{c^{\prime}}(q)\right)\Phi_{l^{\prime}c^{\prime}}^{n^{\prime}}(q;E),\end{split}$
(40)

again dependent on the kernel function defined in Eq. (37).

###  II.6 Poles in the transition matrix

As follows from Eq. (33), the function $\tau_{l}^{n}\left(Z_{c}(q)\right)$ has
a pole whenever $q=\bar{q}_{lc}^{n}$, where,

$\displaystyle\begin{split}\bar{q}_{lc}^{n}=\sqrt{\frac{4m}{3\hbar^{2}}\left(E-\varepsilon_{c}-\varepsilon_{lc}^{n}\right)},\end{split}$
(41)

correspondent to the momenta where the energy of the dimer subsystem is
exactly equal to a two-body bound state. As is typical in scattering theory we
circumvent these poles by adding an infinitesimal imaginary part $i0$ to the
energy [40, 47], after which we can write,

$\displaystyle\begin{split}\tau_{l}^{n}\left(Z_{c}(q\approx\bar{q}_{lc}^{n})\right)\approx\frac{4m}{3\hbar^{2}}\frac{1}{(\bar{q}_{lc}^{n})^{2}-q^{2}+i0}.\end{split}$
(42)

This ensures that the integrals in the AGS and bound-state equations converge
to a finite but complex value. In numerical practice we adjust our integral
equations by subtracting and adding all integrals over the divergent parts of
the transition matrix. The subtracted term then ensures convergence of the
discretized integral, while the added term is evaluated analytically via the
Sokhotski-Plemelj theorem, such that one obtains,

$\displaystyle\begin{split}&8\pi\sum_{l^{\prime}}\sum_{n^{\prime}}\sum_{c^{\prime}}\int
dq\ q^{2}\
\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,q\right)\tau_{l^{\prime}}^{n^{\prime}}\left(Z_{c^{\prime}}(q)\right)\bar{U}_{l^{\prime}c^{\prime}}^{n}(q)=\\\
&\hskip 85.35826pt8\pi\sum_{l^{\prime}}\sum_{n^{\prime}}\sum_{c^{\prime}}\int
dq\ q^{2}\
\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,q\right)\tau_{l^{\prime}}^{n^{\prime}}\left(Z_{c^{\prime}}(q)\right)\bar{U}_{l^{\prime}c^{\prime}}^{n^{\prime}}(q)\\\
&\hskip
85.35826pt-8\pi\frac{4m}{3\hbar^{2}}\sum_{l^{\prime}}\sum_{n^{\prime}}^{\mathrm{poles}}\sum_{c^{\prime}}\int
dq(\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}})^{2}\
\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}}\right)\frac{1}{\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}2}-q^{2}}\bar{U}_{l^{\prime}c^{\prime}}^{n^{\prime}}(\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}})\\\
&\hskip
85.35826pt+8\pi\frac{4m}{6\hbar^{2}}\sum_{l^{\prime}}\sum_{n^{\prime}}^{\mathrm{poles}}\sum_{c^{\prime}}\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}}\
\mathcal{Z}_{ll^{\prime};cc^{\prime}}^{nn^{\prime}}\left(p,\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}}\right)\left[\ln\left(\frac{q_{\mathrm{max}}+\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}}}{q_{\mathrm{max}}-\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}}}\right)-i\pi\right]\bar{U}_{l^{\prime}c^{\prime}}^{n^{\prime}}(\bar{q}_{l^{\prime}c^{\prime}}^{n^{\prime}})\end{split}$
(43)

Here the terms $n\in\mathrm{poles}$ are associated with dimer energies
$\varepsilon_{lc}^{n}$ that fall in the range,

$\displaystyle\begin{split}E-\varepsilon_{c}-\frac{3}{4}\frac{\hbar^{2}q_{\mathrm{max}}^{2}}{m}<\varepsilon_{lc}^{n}<E-\varepsilon_{c}-\frac{3}{4}\frac{\hbar^{2}q_{\mathrm{min}}^{2}}{m}.\end{split}$
(44)

Here $q_{\mathrm{min(max)}}$ is the minimum(maximum) momenta in the numerical
grid used for $q$ (see next section). At the momenta $\bar{q}_{lc}^{n}$ the
three-body recombination process is exactly on-shell. Hence, to obtain the
three-body recombination rate, we sum all
$\absolutevalue{\bar{U}_{cl}^{n}(\bar{q}_{lc}^{n})}^{2}$. For the bound state
energies we use Eq. (43) with
$\bar{U}_{lc}^{n}(p)\rightarrow\Phi_{lc}^{n}(p)$.

###  II.7 Numerical method and parameters

To solve Eqs. (36) and (40) we first discretize the atom-dimer momenta $p$ and
$q$ using Gauss-Legendre quadrature [53], segmented into three sections as
$q_{1}r_{\mathrm{vdW}}\in\left[10^{-5},10^{-2}\right]$,
$q_{2}r_{\mathrm{vdW}}\in\left[10^{-2},1\right]$ and
$q_{3}r_{\mathrm{vdW}}\in\left[1,q_{\mathrm{max}}r_{\mathrm{vdW}}\right]$.
Here $q_{\mathrm{max}}$ acts as an ultraviolet cut-off, which limits the
maximum binding energy of dimers that couple to the three-body complex. After
discretization, both the AGS and Faddeev equations may be written as
straightforward matrix equations,

$\displaystyle\begin{split}\mathbf{U}=\mathbf{Z}_{0}+\underline{Z}\
\underline{\tau}\ \mathbf{U},\end{split}$ (45)

where $\mathbf{U}$ is a vector containing the elements $\bar{U}_{lc}^{n}(p)$
or $\Phi_{lc}^{n}(p)$ for the AGS or Faddeev equation respectively. The vector
$\mathbf{Z}_{0}$ contains the inhomogeneous term, which vanishes in the
Faddeev equation. The matrix $\underline{Z}$ contains the kernel function
$\mathcal{Z}$, and the diagonal matrix $\underline{\tau}$ contains the values
of $\tau$. In numerical practice, we are typically limited by the considerable
memory required to store the matrix $\mathcal{Z}$, which restrains us both in
the ultraviolet cut-off $q_{\mathrm{max}}$ and the maximum orbital quantum
number $l_{\mathrm{max}}$. We can however reduce the dimension of the matrix
equation by separating all terms into two distinct spaces, after which we can
write a matrix equation for the terms in space 1 as,

$\displaystyle\begin{split}\mathbf{U}_{1}=\tilde{\mathbf{Z}}_{0}+\underline{\tilde{Z}}\
\underline{\tau}_{1}\ \mathbf{U}_{1},\end{split}$ (46)

where,

$\displaystyle\begin{split}\tilde{\mathbf{Z}}_{0}&=\mathbf{Z}_{0,1}+\underline{Z}_{12}\left(\underline{\tau}_{2}^{-1}-\underline{Z}_{22}\right)^{-1}\mathbf{Z}_{0,2}\\\
\underline{\tilde{Z}}&=\underline{Z}_{11}+\underline{Z}_{12}\left(\underline{\tau}_{2}^{-1}-\underline{Z}_{22}\right)^{-1}\underline{Z}_{21},\end{split}$
(47)

Such a two-space approach is more generally used in quantum physics when one
desires to eliminate a subspace of states whilst retaining the effects of
their coupling to the states that remain, including for example the formal
theory of scattering resonances [10], and the mapping between the Hubbard and
Heisenberg Hamiltonians in condensed matter physics [54].

Assuming that the terms $\tau_{2}$ are small, we can make the first order
approximation,

$\displaystyle\begin{split}\tilde{\mathbf{Z}}_{0}&\approx\mathbf{Z}_{0,1}+\underline{Z}_{12}\underline{\tau}_{2}\mathbf{Z}_{0,2}\\\
\underline{\tilde{Z}}&\approx\underline{Z}_{11}+\underline{Z}_{12}\underline{\tau}_{2}\underline{Z}_{21},\end{split}$
(48)

which significantly reduces the required memory as we do not need the elements
$\underline{Z}_{22}$. Essentially the method above amounts to another
reduction on the number of terms in the expansion (19) for the transition
matrix, although now the terms are included to first order instead of
neglected entirely.

Although the Faddeev equation looks superficially simpler than the AGS
equation, finding its solution requires some iterative procedure in the three-
body energy $E$ until one obtains a bound state energy. This iteration is
complicated by the fact that the Efimov trimers we look for are embedded in
the atom-dimer continua associated with deep dimer states. Hence, the Efimov
state is not a true bound state but rather a three-body resonance with complex
energy $E_{*}=\mathrm{Re}(E_{*})-i\Gamma/2$, where $\Gamma$ is the energetic
width of the state which defines a characteristic lifetime
$\tau_{\mathrm{trimer}}=\hbar/\Gamma$. To simplify the search for bound
states, we assume the Efimov state to be long lived enough such that
$\mathrm{Re}(E_{*})\gg\Gamma$. Then we scan the absolute logarithm of the
determinant
$D(E)=\mathrm{det}\left(\underline{\tilde{Z}}\underline{\tau}_{1}-\mathbb{I}\right)$
along the real axis of $E$, noting that $D(E_{*})=0$. Upon passing this root,
we can expand the determinant to lowest order in $\Gamma$,

$\displaystyle\begin{split}\mathrm{ln}\absolutevalue{D(E)}\approx\mathrm{ln}\absolutevalue{\derivative{D}{E}}_{E=E_{*}}+\mathrm{ln}\left(\sqrt{(E-\mathrm{Re}(E_{*}))^{2}+\frac{1}{4}\Gamma^{2}}\right).\end{split}$
(49)

It follows that the location of a minimum in $\mathrm{ln}\absolutevalue{D(E)}$
corresponds with the value of $\mathrm{Re}(E_{*})$, which we extract to obtain
the energy of the Efimov trimer. Note that this approach will only work if the
root in the determinant lies close enough to the real axis such that Eq. (49)
is valid. Following experiment [34], we will be primarily interested in the
trimer-dimer energy difference $\Delta
E=\varepsilon_{\mathrm{res}}-\mathrm{Re}(E_{*})$, measuring the three-body
energy relative to the energy of the Feshbach dimer
$\varepsilon_{\mathrm{res}}$. To illustrate our method and the general
behavior of $\Delta E$ near a two-body resonance, we plot in Fig. 3 scans of
$\mathrm{ln}\absolutevalue{D(E)}$ along the appropriately shifted energy axis
$\varepsilon_{\mathrm{res}}-E$. We use here a single-channel van der Waals
potential at the 8’th $s$-wave potential resonance [23], allowing us to
compute the determinant at a large number of points thus illustrating the
appearance of sharp local minima which we identify as Efimov resonances,
fitted well by Eq. (49). Plotting the resulting values of the trimer-dimer
energy difference $\Delta E$ as a function of scattering length then shows the
characteristic behavior of the Efimov trimer, merging with the atom-dimer
continuum at the critical scattering length $a_{*}\approx 1.1a_{+}$ [3, 29].

Figure 3:  Illustration of the numerical solution to the Faddeev equation for
an Efimov resonance in a single-channel van der Waals potential near the 8’th
two-body potential resonance. In the left-hand panel we plot scans of
$\mathrm{ln}\absolutevalue{D(E)}$ as a function of the shifted three-body
energy $\varepsilon_{\mathrm{res}}-E$, for a set of (positive) scattering
lengths. Local minima in $\mathrm{ln}\absolutevalue{D(E)}$ signify the
location of Efimov resonances, as follows from Eq. (49). Solid lines show fits
of this equation to the minima. In the right-hand panel we show the resulting
dependence of the trimer-dimer energy difference $\Delta E$ on the scattering
length, where the scatter points match with the fits in the left-hand panel.

###  II.8 Spin models

An inherent advantage of our method, based on integral equations, is that it
is reasonably simple to exclude specific three-body channels from the
computation if so desired. Up until this point, the method that we have
formulated is quantum mechanically rigorous, except for the omission of non-
additive three-body forces, which we assume to be negligible in the systems of
interest [4]. We will refer to this as the Full Multichannel Spin (FMS) model
[55, 37]. In the literature, FMS calculations are rare, and it is instead
common practice to reduce the spin degrees-of-freedom of the total system by
various approximations, see e.g. [8, 19, 56, 20, 57, 58, 21, 25, 22, 59]. We
will study a specific simplifying approach to the spin structure known as the
Fixed Spectating Spin (FSS) approximation [55], which turns off the two-body
interaction whenever the third particle (the spectator) is not in the incoming
spin state $c_{\mathrm{in}}$. Formally, this is achieved through the
replacement $\hat{V}_{\alpha}\rightarrow\hat{V}_{\alpha}^{\mathrm{FSS}}$,
where,

$\displaystyle\begin{split}\hat{V}_{\alpha}^{\mathrm{FSS}}=\hat{V}_{\alpha}\outerproduct*{c_{\mathrm{in}}}{c_{\mathrm{in}}}.\end{split}$
(50)

It follows that the two-body transition matrix similarly vanishes whenever the
spectating particle has changed spin, meaning that in both the AGS and Faddeev
equations we can simply set $c\rightarrow c_{\mathrm{in}}$ everywhere. This
leads to a significant reduction in the dimension of the integral equations,
such that the resulting numerical computation is much less expensive. The
validity of the FSS approximation relies on the characteristic suppression of
three-body probability in the short-range, which should make it relatively
unlikely for three particles to approach within the spin-exchange regime
simultaneously. Numerical studies have shown however that for some systems,
the effects of three-body spin-exchange where all three particles change their
spin can still be significant [55, 37], and are essential in reproducing the
correct three-body parameter for 7Li specifically [35]. In the next section we
will expand on these results, comparing the FSS and FMS approaches also for
positive scattering lengths.

##  III Results

As an application of our model, we will study the rate of three-body
recombination $K_{3}$ near two high-field Feshbach resonances in 7Li, whose
characteristic properties we collect in table 1. Both resonances are of
intermediate width with $s_{\mathrm{res}}\sim 1$ [11], well into the non-
universal regime for the three-body parameter [23].

State |  $B_{0}$ [G] |  $a_{\mathrm{bg}}$ $[a_{0}]$ |  $\Delta B$ [G] | $s_{\mathrm{res}}$ | $M_{F}$
---|---|---|---|---|---
$aa$ | 738.2(2) | -18.24 | -237.8 | 0.56 | $3$
$bb$ | 893.7(3) | -20.98 | -171.0 | 0.493 | $0$
Table 1: Parameters of the two Feshbach resonances studied in this work [29,
30, 60]. Here $B_{0}$ is the resonant magnetic field, $a_{\mathrm{bg}}$ the
background scattering length, $\Delta B$ the magnetic field width,
$s_{\mathrm{res}}$ the resonance strength parameter[11], and $M_{F}$ the total
magnetic projection quantum number for the three-body calculation.

As we are interested in Efimovian signatures, it is useful to recast our
results for $K_{3}$ using the theoretical prediction
$K_{3}(a)=3C_{\pm}(a)\hbar a^{4}/m$ [3]. Here the function $C_{\pm}(a)$
captures the Efimovian modulation of the universal $a^{4}$ scaling of $K_{3}$
near resonance, for positive ($+$) and negative ($-$) scattering lengths
respectively. Closed form expressions for $C_{\pm}(a)$ can be obtained within
a zero-range effective field theory upon introducing a three-body parameter
[3, 13, 29],

$\displaystyle C_{-}(a)$
$\displaystyle=4590\frac{\mathrm{\sinh}(2\eta_{-})}{\mathrm{sin}^{2}(s_{0}\ln(a/a_{-}))+\mathrm{sinh}^{2}(\eta_{-})},$
(51) $\displaystyle C_{+}(a)$
$\displaystyle=67.1e^{-2\eta_{+}}\left(\mathrm{cos}^{2}(s_{0}\ln(a/a_{+}))+\mathrm{sinh}^{2}(\eta_{+})\right)+16.8(1-e^{-4\eta_{+}}).$
(52)

We will assume a priori that our calculated values of $K_{3}(a)$ obey the
$a^{4}$ power law, and subsequently extract the associated Efimovian functions
$C_{\pm}(a)=mK_{3}(a)/3\hbar a^{4}$. We can then fit the data directly using
Eqs. (51) and (52), thus checking the validity of these universal expressions
in 7Li, and obtaining explicit values for the three-body parameters $a_{\pm}$,
with associated trimer widths $\eta_{\pm}$. The results, for both Feshbach
resonances, are shown in Fig. 4 and table 2. Here we note that the results for
$K_{3}$ for negative scattering lengths have been published previously
elsewhere [35].

Figure 4:  Efimovian enhancement $C_{\pm}(a)$ of the three-body recombination rate near both Feshbach resonances in table 1. FMS results (red) and FSS results (blue) are compared with the experimental data of Ref. [29] (grey). Dash-dotted lines show fits of the data to the predictions of effective field theory in Eqs. (51) and (52), which we use to extract the three-body parameters $a_{\pm}$ and trimer widths $\eta_{\pm}$. State | Model | $a_{-}/r_{\mathrm{vdW}}$ | $\eta_{-}$ | $a_{+}/r_{\mathrm{vdW}}$ | $\eta_{+}$ | $a_{+}/\absolutevalue{a_{-}}$
---|---|---|---|---|---|---
| FSS | -40.5 | 0.63 | 43.1 | 0.91 | 1.06
$M_{F}=3$ | FMS | -8.79 | 0.27 | 8.69 | 0.13 | 0.99
| Experiment[29] | -8.43(37) | 0.253(62) | 8.16(49) | 0.170(41) | 0.97(8)
| FSS | -36.3 | 0.78 | 33.6 | 0.35 | 0.92
$M_{F}=0$ | FMS | -12.7 | 0.38 | 20.0 | 0.37 | 1.57
| Experiment[29] | -8.62(37) | 0.180(48) | 7.33(77) | 0.213(79) | 0.85(11)
Table 2: Three-body parameters $a_{\pm}$ and trimer widths $\eta_{\pm}$
obtained from the fits to universal expressions for $C_{\pm}(a)$, as shown in
Fig. 4, compared with similar fits to the experimental data.

Our results show significant differences between FSS and FMS spin models,
clearly indicating that the FSS approximation is inadequate for 7Li atoms in
the strongly interacting regime. For negative scattering lengths, where no
two-body bound state exists, we find that the FSS approximation produces a
broad Efimov resonance at relatively large scattering length, in accordance
with the general expectation that the Efimov trimer is shifted towards the
three-body threshold for intermediate strength Feshbach resonances [20, 21,
22, 23]. In contrast, upon going to the FMS model we observe a dramatic shift
of the position of the Efimov resonance to lower absolute scattering lengths,
with the three-body parameter for the $M_{F}=3$ case even receeding below the
universal van der Waals value. The resonance itself is also significantly
narrowed, indicating a non-trivial stabilizing effect of the additional spin
channels on the Efimov trimer. Both these effects significantly improve on the
match between the theoretical calculations and the experimental data,
indicating the strong presence of three-body spin-exchange in experiments, and
the accuracy of our numerical model. One will note however that for the case
$M_{F}=0$, a significant discrepancy between the experimental data and FMS
results remains. We attribute this mismatch to a lack of convergence of our
numerical method, originating from the fact that the $M_{\mathrm{tot}}=0$
channel has a much larger number of coupled three-body channels, which limits
us in the choice of numerical parameters. Specifically, for the $M_{F}=3$
channel we use $l_{\mathrm{max}}=12$ and
$q_{\mathrm{max}}r_{\mathrm{vdW}}=40$, which we have confirmed gives adequate
convergence in the value of $a_{-}$. For the $M_{F}=0$ channel numerical
limitations force us to reduce the maximum partial wave to
$l_{\mathrm{max}}=8$, which is not sufficient to fully converge the result,
thus explaining the mismatch with the experiment. For further discussion of
this point we refer the reader to the Supplemental Material of Ref. [35].

For positive scattering lengths, we find once more that the value of $a_{+}$
is decreased significantly in the FMS model compared to an FSS calculation.
For the $M_{F}=3$ channel, this shift again leads to an excellent match with
the experiment of Ref. [29], while the associated FSS results fail to
reproduce the data to any satisfactory degree. In addition, the ratio
$a_{+}/\absolutevalue{a_{-}}$ of our FMS results matches well with both the
experimental data and the theoretical prediction
$a_{+}/\absolutevalue{a_{-}}=0.96(3)$ [3, 29]. A similar shift of the Efimov
resonance is seen in the $M_{F}=0$ channel, although as before the numerical
complexity in this more highly coupled three-body channel does not allow us to
fully converge the calculations. As for negative scattering lengths however,
the general trend of our numerics suggest that a converged calculation will
give a further decrease of $a_{+}$, towards the experimental data.

Next to the comparison with experiment, our results for positive scattering
length additionally reveal that the universal prediction for $C_{+}(a)$, as
given in Eq. (52), can not give a good fit to the recombination data over the
full range of scattering lengths we study, indicating the presence of finite
range effects which are unaccounted for in the universal theory. This effect
is especially prevalent for the FSS results, where the fit of Eq. (52) is very
poor leading also to a large uncertainty in $a_{+}$ and $\eta_{+}$. This
observation however, is fully in line with the experimental data, which is
also quite poorly fitted by the universal expression, and actually finds
better match with our numerical computations. Similarly, in Ref. [30] it was
found that a good fit of the Efimov resonance for $a>0$ could only be obtained
by adding an overall scaling factor to Eq. (52), which is also suggested by
our FMS data in the $M_{F}=3$ channel. These discrepancies indicate that an
extended form of Eq. (52) which accounts for finite range corrections may be
required to match the data. Approaches for adding such corrections to the
universal Efimovian theory have indeed been discussed in the literature, see
for example Ref. [61].

##  IV Conclusion

In this article we have described a state-of-the-art multichannel method for
computing the three-body recombination rate and three-body bound state energy
in a system of alkali-metal atoms near a Feshbach resonance, using realistic
molecular pairwise interaction potentials. We have subsequently applied this
method to study Efimovian signatures near two experimentally relevant Feshbach
resonances in 7Li, where we show strong influence of spin-exchange effects
which are typically neglected in the three-body literature. In particular, we
find that such effects are essential in reproducing experimental data.

Our work poses several opportunities and challenges for further research.
First it is clear that the numerical complexity of simulating interacting
three-body systems at the level of accuracy required to reproduce experiment
is considerable, requiring expensive high-performance computing facilities and
considerable time investments. These points motivate further searches for
accurate simplifying assumptions or improvements on the numerical efficiency
of our method, which may then allow us to fully converge all numerical
computations and thus close the remaining gaps between theory and experiment,
which are highlighted once more in the present article. It may also allow for
an application of the bound state method introduced in this work to realistic
multichannel systems, which poses an even more challenging problem numerically
due to the larger number of three-body energies required.

Looking towards the future, our model can be applied generally in the search
for non-universal features in three-body physics in ultracold atomic gases. We
have shown that our method, once converged, accurately reproduces the
experimentally observed three-body recombination rate, even when this rate
deviates strongly from the universal expectation. We have thus developed a
tool which can be used to guide future experiments, by predicting the three-
body loss rate in general systems to a high degree of accuracy. For example,
quantum gas experiments and technologies are typically limited by the lifetime
of the gas [62, 63], and thus may wish to exploit non-universal minima in the
three-body loss rate, which can be searched for using our method. Finally, it
is interesting to extend our method to also capture elastic three-body
scattering, which is known to significantly affect the phase diagram of Bose-
Einstein condensates in the regime of small scattering lengths, and can
provide a mechanism for the stabilization of quantum droplets [64, 65, 66,
67]. The direct observation and characterization of this scattering remains an
elusive challenge to contemporary experiments, making accurate numerical
models invaluable.

## Acknowledgments

This work is financially supported by the Dutch Ministry of Economic Affairs
and Climate Policy (EZK), as part of the Quantum Delta NL program, and by the
Netherlands Organisation for Scientific research (NWO) under Grant No.
680.92.18.05. The results presented in this work were obtained on the Dutch
national supercomputer Snellius, with support from the Eindhoven
Supercomputing Center and SURF.

## Statements and Declarations

The authors have no competing interests to declare that are relevant to the
content of this article.

## References

  * Efimov [1970] V. Efimov, Phys. Lett. B 33, 563 (1970).
  * Efimov [1971] V. Efimov, Sov. J. Nucl. Phys. 12, 589 (1971).
  * Braaten and Hammer [2006] E. Braaten and H.-W. Hammer, Physics Reports 428, 259 (2006).
  * Naidon and Endo [2017] P. Naidon and S. Endo, Rep. Prog. Phys. 80, 056001 (2017).
  * Greene  _et al._ [2017] C. H. Greene, P. Giannakeas, and J. Pérez-Ríos, Rev. Mod. Phys. 89, 035006 (2017).
  * D’Incao [2018] J. P. D’Incao, J. Phys. B: At. Mol. Opt. Phys. 51, 043001 (2018).
  * D'Incao  _et al._ [2009] J. P. D'Incao, C. H. Greene, and B. D. Esry, J. Phys. B: At. Mol. Opt. Phys. 42, 044016 (2009).
  * Petrov [2004] D. S. Petrov, Phys. Rev. Lett. 93, 143201 (2004).
  * Feshbach [1958] H. Feshbach, Ann. Phys. 5, 357 (1958).
  * Feshbach [1962] H. Feshbach, Ann. Phys. 19, 287 (1962).
  * Chin  _et al._ [2010] C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Esry  _et al._ [1999] B. D. Esry, C. H. Greene, and J. P. Burke, Phys. Rev. Lett. 83, 1751 (1999).
  * 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, and R. Grimm, 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).
  * 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).
  * Wang  _et al._ [2012] J. Wang, J. P. D’Incao, B. D. Esry, and C. H. Greene, Phys. Rev. Lett. 108, 263001 (2012).
  * Naidon  _et al._ [2014a] P. Naidon, S. Endo, and M. Ueda, Phys. Rev. A 90, 022106 (2014a).
  * Naidon  _et al._ [2014b] P. Naidon, S. Endo, and M. Ueda, Phys. Rev. Lett. 112, 105301 (2014b).
  * Gogolin  _et al._ [2008] A. O. Gogolin, C. Mora, and R. Egger, Phys. Rev. Lett. 100, 140404 (2008).
  * Schmidt  _et al._ [2012] R. Schmidt, R. Rath, and W. Zwerger, Eur. Phys. J. B 85, 386 (2012).
  * Langmack  _et al._ [2018] C. Langmack, R. Schmidt, and W. Zwerger, Phys. Rev. A 97, 033623 (2018).
  * Secker  _et al._ [2021a] T. Secker, D. J. M. Ahmed-Braun, P. M. A. Mestrom, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 103, 052805 (2021a).
  * van de Kraats  _et al._ [2023] J. van de Kraats, D. J. M. Ahmed-Braun, J.-L. Li, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 107, 023301 (2023).
  * Roy  _et al._ [2013] S. Roy, M. Landini, A. Trenkwalder, G. Semeghini, G. Spagnolli, A. Simoni, M. Fattori, M. Inguscio, and G. Modugno, Phys. Rev. Lett. 111, 053202 (2013).
  * Chapurin  _et al._ [2019] R. Chapurin, X. Xie, M. J. Van de Graaff, J. S. Popowski, J. P. D’Incao, P. S. Julienne, J. Ye, and E. A. Cornell, Phys. Rev. Lett. 123, 233402 (2019).
  * Etrych  _et al._ [2023] J. c. v. Etrych, G. Martirosyan, A. Cao, J. A. P. Glidden, L. H. Dogra, J. M. Hutson, Z. Hadzibabic, and C. Eigen, Phys. Rev. Res. 5, 013174 (2023).
  * 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).
  * Gross  _et al._ [2011] N. Gross, Z. Shotan, O. Machtey, S. Kokkelmans, and L. Khaykovich, C. R. Phys. 12, 4 (2011).
  * Dyke  _et al._ [2013] P. Dyke, S. E. Pollack, and R. G. Hulet, Phys. Rev. A 88, 023625 (2013).
  * Yudkin and Khaykovich [2021] Y. Yudkin and L. Khaykovich, Phys. Rev. A 103, 063303 (2021).
  * D’Incao  _et al._ [2022] J. D’Incao, Y. Yudkin, P. Julienne, and L. Khaykovich, Bulletin of the American Physical Society (2022).
  * Yudkin  _et al._ [2023] Y. Yudkin, P. S. Julienne, and L. Khaykovich, Phys. Rev. A 107, 053303 (2023).
  * Yudkin  _et al._ [2024] Y. Yudkin, R. Elbaz, J. P. D’Incao, P. S. Julienne, and L. Khaykovich, Nat. Commun. 15, 2127 (2024).
  * van de Kraats  _et al._ [2024] J. van de Kraats, D. J. M. Ahmed-Braun, J.-L. Li, and S. J. J. M. F. Kokkelmans, Phys. Rev. Lett. 132, 133402 (2024).
  * Sitenko [1991] A. Sitenko,  _Scattering theory_ , Springer Series in Nuclear and Particle Physics (Springer-Verlag, Berlin Heidelberg, 1991).
  * Li  _et al._ [2022] J.-L. Li, T. Secker, P. M. A. Mestrom, and S. J. J. M. F. Kokkelmans, Phys. Rev. Research 4, 023103 (2022).
  * Stoof  _et al._ [1988] H. T. C. Stoof, J. M. V. A. Koelman, and B. J. Verhaar, Phys. Rev. B 38, 4688 (1988).
  * Julienne and Hutson [2014] P. S. Julienne and J. M. Hutson, Phys. Rev. A 89, 052715 (2014).
  * Taylor [2006] J. Taylor,  _Scattering Theory: The Quantum Theory of Nonrelativistic Collisions_ , 1st ed. (Dover Publications, Inc., New York, 2006).
  * Glöckle [1983] W. Glöckle,  _The quantum mechanical few-body problem_, Texts and monographs in physics (Springer, Berlin, 1983).
  * Esry  _et al._ [2001] B. D. Esry, C. H. Greene, and H. Suno, Phys. Rev. A 65, 010705(R) (2001).
  * Light  _et al._ [1985] J. C. Light, I. P. Hamilton, and J. V. Lill, J. Chem. Phys. 82, 1400 (1985).
  * Fattal  _et al._ [1996] E. Fattal, R. Baer, and R. Kosloff, Phys. Rev. E 53, 1217 (1996).
  * Willner  _et al._ [2004] K. Willner, O. Dulieu, and F. Masnou-Seeuws, J. Chem. Phys. 120, 548 (2004).
  * Secker  _et al._ [2021b] T. Secker, J.-L. Li, P. M. A. Mestrom, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 103, 032817 (2021b).
  * Mestrom  _et al._ [2019a] P. M. A. Mestrom, T. Secker, R. M. Kroeze, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 99, 012702 (2019a).
  * Lee  _et al._ [2007] M. D. Lee, T. Köhler, and P. S. Julienne, Phys. Rev. A 76, 012720 (2007).
  * Moerdijk  _et al._ [1996] A. J. Moerdijk, H. M. J. M. Boesten, and B. J. Verhaar, Phys. Rev. A 53, 916 (1996).
  * Smirne  _et al._ [2007] G. Smirne, R. M. Godun, D. Cassettari, V. Boyer, C. J. Foot, T. Volz, N. Syassen, S. Dürr, G. Rempe, M. D. Lee, K. Góral, and T. Köhler, Phys. Rev. A 75, 020702(R) (2007).
  * Alt  _et al._ [1967] E. Alt, P. Grassberger, and W. Sandhas, Nuclear Physics B 2, 167 (1967).
  * Faddeev [1961] L. Faddeev, Sov. Phys. JETP 12, 1014 (1961).
  * Press and Teukolsky [2007] W. H. Press and S. A. Teukolsky,  _Numerical Recipes 3rd Edition: The Art of Scientific Computing_ (Cambridge University Press, 40 W. 20 St. New York, NY, United States, 2007).
  * Mila and Schmidt [2011] F. Mila and K. P. Schmidt, in  _Introduction to Frustrated Magnetism_, Springer Series in Solid-State Sciences, Vol. 164, edited by C. Lacroix, P. Mendels, and F. Mila (Springer, Berlin, Heidelberg, 2011).
  * Secker  _et al._ [2021c] T. Secker, J.-L. Li, P. M. A. Mestrom, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 103, 022825 (2021c).
  * Massignan and Stoof [2008] P. Massignan and H. T. C. Stoof, Phys. Rev. A 78, 030701(R) (2008).
  * Wang and Julienne [2014] Y. Wang and P. Julienne, Nat. Phys. 10, 768–773 (2014).
  * Wolf  _et al._ [2017] J. Wolf, M. Deiß, A. Krükow, E. Tiemann, B. P. Ruzic, Y. Wang, J. P. D’Incao, P. S. Julienne, and J. H. Denschlag, Science 358, 921 (2017).
  * Haze  _et al._ [2022] S. Haze, J. P. D’Incao, D. Dorer, M. Deiß, E. Tiemann, P. S. Julienne, and J. H. Denschlag, Phys. Rev. Lett. 128, 133401 (2022).
  * Jachymski and Julienne [2013] K. Jachymski and P. S. Julienne, Phys. Rev. A 88, 052701 (2013).
  * Gattobigio  _et al._ [2019] M. Gattobigio, M. Göbel, H.-W. Hammer, and A. Kievsky, Few-Body Syst. 60, 40 (2019).
  * Kunimi and Danshita [2019] M. Kunimi and I. Danshita, Phys. Rev. A 99, 043613 (2019).
  * Rätzel and Schützhold [2021] D. Rätzel and R. Schützhold, Phys. Rev. A 103, 063321 (2021).
  * Bulgac [2002] A. Bulgac, Phys. Rev. Lett. 89, 050402 (2002).
  * Tan [2008] S. Tan, Phys. Rev. A 78, 013636 (2008).
  * Mestrom  _et al._ [2019b] P. M. A. Mestrom, V. E. Colussi, T. Secker, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 100, 050702 (2019b).
  * Mestrom  _et al._ [2020] P. M. A. Mestrom, V. E. Colussi, T. Secker, G. P. Groeneveld, and S. J. J. M. F. Kokkelmans, Phys. Rev. Lett. 124, 143401 (2020).