#  _Ab initio_ calculation of the spectrum of Feshbach resonances in NaLi + Na
collisions

Tijs Karman  Institute for Molecules and Materials, Radboud University,
Heijendaalseweg 135, 6525 AJ Nijmegen, the Netherlands  Marcin Gronowski
Faculty of Physics, University of Warsaw, Pasteura 5, 02-093 Warsaw, Poland
Michał Tomza  Faculty of Physics, University of Warsaw, Pasteura 5, 02-093
Warsaw, Poland  Juliana J. Park  Research Laboratory of Electronics, MIT-
Harvard Center for Ultracold Atoms, Department of Physics, Massachusetts
Institute of Technology, Cambridge, Massachusetts 02139, USA  Hyungmok Son
Research Laboratory of Electronics, MIT-Harvard Center for Ultracold Atoms,
Department of Physics, Massachusetts Institute of Technology, Cambridge,
Massachusetts 02139, USA  Department of Physics, Harvard University,
Cambridge, Massachusetts 02138, USA  Yu-Kun Lu  Research Laboratory of
Electronics, MIT-Harvard Center for Ultracold Atoms, Department of Physics,
Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA
Alan O. Jamison  Institute for Quantum Computing and Department of Physics &
Astronomy, University of Waterloo, Waterloo, Ontario N2L 3G1, Canada  Wolfgang
Ketterle  Research Laboratory of Electronics, MIT-Harvard Center for Ultracold
Atoms, Department of Physics, Massachusetts Institute of Technology,
Cambridge, Massachusetts 02139, USA

###### Abstract

We present a combined experimental and theoretical study of the spectrum of
magnetically tunable Feshbach resonances in NaLi $(a^{3}\Sigma^{+})$ $+$ Na
collisions. In the accompanying paper, we observe experimentally 8 and 17
resonances occur between $B=0$ and $1400$ G in upper and lower spin-stretched
states, respectively. Here, we perform _ab initio_ calculations of the NaLi
$+$ Na interaction potential and describe in detail the coupled-channel
scattering calculations of the Feshbach resonance spectrum. The positions of
the resonances cannot be predicted with realistic uncertainty in the state-of-
the-art _ab initio_ potential, but our calculations yield a _typical_ number
of resonances that is in near-quantitative agreement with experiment. We show
that the main coupling mechanism results from spin-rotation and spin-spin
couplings in combination with the anisotropic atom-molecule interaction. The
calculations furthermore explain the qualitative difference between the
numbers of resonances in either spin state.

##  I Introduction

Magnetically tunable resonances have become an indispensable tool for
controlling the interactions between ultracold atoms Chin  _et al._ (2010).
Control over the interactions enables using atoms as quantum simulators to
study for example condensed matter physics Bloch  _et al._ (2012); Bakr  _et
al._ (2010); Struck  _et al._ (2011); Simon  _et al._ (2011); Chiu  _et al._
(2018). Feshbach resonances are also used to associate pairs of ultracold
atoms into weakly bound molecules Köhler  _et al._ (2006). These weakly bound
molecules can be transferred coherently to their absolute ground state using a
stimulated Raman adiabatic passage (STIRAP) Bergmann  _et al._ (1998), which
has become a common tool for producing ultracold ground-state molecules Ni
_et al._ (2008); Takekoshi  _et al._ (2014); Molony  _et al._ (2014); Guo  _et
al._ (2016). Resonances in molecular collisions could similarly enable the
control of contact interactions between molecules, or between molecules and
atoms. This could enable sympathetic cooling of molecules Son  _et al._
(2020), and the bottom up construction of ultracold polyatomic molecules Yang
_et al._ (2022). Furthermore, scattering resonances are incredibly sensitive
to the details of the interaction potential and provide an interesting testing
ground for theory Gronowski  _et al._ (2020). In particular, there is an
active debate about the nature of so-called “sticky collisions” between
ultracold molecules Mayle  _et al._ (2012, 2013); Christianen  _et al._
(2019a), their role in collisional losses Christianen  _et al._ (2019b), and
what simplified statistical models can be used to describe these Mayle  _et
al._ (2012, 2013); Croft  _et al._ (2021); Christianen  _et al._ (2021).
Observation of atom-molecule and molecule-molecule resonances could directly
test some of these models by measuring the “effective” density of states
Christianen  _et al._ (2019a); Man  _et al._ (2022), the probability of short-
range loss Croft  _et al._ (2020); Christianen  _et al._ (2021), and by
probing the mechanism of this loss.

Collisional loss also represents a major hurdle in the observation of Feshbach
resonances. The resonance states correspond to quasi-bound states of the atom-
molecule or molecule-molecule collision complex. In the presence of
collisional loss, however, these states experience decay and the resonances
broaden Christianen  _et al._ (2021). In the case of so-called universal loss
Idziaszek and Julienne (2010), which is consistent with many experimental
observations of molecule-molecule Ni  _et al._ (2008); Danzl  _et al._ (2010);
Takekoshi  _et al._ (2014); Molony  _et al._ (2014); Park  _et al._ (2015);
Guo  _et al._ (2016); Rvachov  _et al._ (2017); Seeßelberg  _et al._ (2018);
Yang  _et al._ (2019) loss rates, the short-range loss becomes complete and
the resonances disappear entirely. Universal loss is also observed for some
but not all ultracold molecule-atom collisions Gregory  _et al._ (2021); Wang
_et al._ (2021); Nichols  _et al._ (2022); Voges  _et al._ (2022). For these
reasons it is an open question for which systems scattering resonances exist.

Recently, scattering resonances have been observed in NaK+K Yang  _et al._
(2019), NaLi+Na Son  _et al._ (2022), and even NaLi+NaLi Park  _et al._ (2022)
collisions, but the three systems are very different. The resonances observed
in NaK+K could be interpreted as long-range states that are further split into
several resonances by hyperfine interactions, and which could tentatively be
assigned free-molecule and free-atom quantum numbers Wang  _et al._ (2021). In
the case of collisions between triplet NaLi molecules and Na atoms in the spin
stretched state, the background loss is much less than universal. This
situation, however, could be unique to triplet NaLi, which is the lightest
bialkali realized that furthermore has relatively weak interactions with Na in
the spin-stretched state. By contrast, NaLi bimolecular collisions exhibit
universal background loss. Nevertheless, a single resonance has been observed
in NaLi+NaLi collisions, though assignment of the observed resonance in terms
of a stable resonance state is complicated Park  _et al._ (2022).

Here, we present a combined experimental and theoretical study of the spectrum
of magnetically tunable Feshbach resonances in NaLi $(a^{3}\Sigma^{+})$ \+ Na
collisions. In the accompanying paper we show 8 and 17 resonances occur
between $B=0$ and $1400$ G in the upper and lower spin-stretched states,
respectively com . We perform _ab initio_ calculations of the NaLi+Na
interaction potential, and use these to perform coupled-channel scattering
calculations of the Feshbach resonance spectrum. We show that the positions of
the resonances cannot be predicted with realistic uncertainty in the state-of-
the-art _ab initio_ potential, but that our calculations can nevertheless
confirm the expected number of resonances in this magnetic field range, and
furthermore explain the qualitative difference between the numbers of
resonances in either spin state, and the nature of the coupling mechanism. In
this paper, we describe in detail how these calculations were performed, and
also present additional supporting calculations of the density of states and
the effect of hyperfine interactions on the spectrum of resonances.

##  II Interaction potential

Intermolecular interactions in the triatomic NaLi+Na system can be formally
decomposed into a pairwise additive two-body part and a pairwise nonadditive
three-body part. The two-body part is the sum of interatomic interaction
between all atomic pairs obtained neglecting the presence of other atoms. The
nonadditive three-body interaction results from the electron density
modification and correlation extending beyond atomic pairs and can also be
understood as the change of pairwise interactions due to the presence of other
atoms.

###  II.1 Two-body interactions

We first consider the two-body interactions. The interaction between two
alkali-metal atoms $A$ and $B$ depends on the distance between the atoms and
whether the spins are singlet or triplet coupled, which can be represented as

$\displaystyle\hat{V}_{AB}(r_{AB})=\bar{V}_{AB}(r_{AB})+\hat{s}_{A}\cdot\hat{s}_{B}\Delta
V_{AB}(r_{AB}),$ (1)

where $\bar{V}_{AB}(r)=[3V_{AB}^{(T)}(r)+V_{AB}^{(S)}(r)]/4$ and $\Delta
V_{AB}(r)=V_{AB}^{(T)}(r)-V_{AB}^{(S)}(r)$. The singlet $V_{AB}^{(S)}$ and
triplet $V_{AB}^{(T)}$ potentials for NaNa and NaLi are taken from Refs. Knoop
_et al._ (2011); Steinke  _et al._ (2012); Rvachov  _et al._ (2018). These are
empirical potentials that reproduce the known molecular spectroscopy and atom-
atom scattering lengths accurately. The total additive electronic interaction
is then computed for fixed molecular bond length, as a function of the
molecule-atom center-of-mass separation $R$ and the Jacobi angle $\theta$ as

$\hat{V}^{2b}(R,\theta)=\sum_{A<B}\hat{V}_{AB}(r_{AB}),$ (2)

where, the three interatomic distances $r_{AB}$ are obtained as a function of
$R$ and $\theta$ and evaluated with Eq. (1) for each pair. This total additive
interaction is represented in a Legendre expansion convenient for scattering
calculations

$\displaystyle\hat{V}^{2b}(R,\theta)$
$\displaystyle=\sum_{L}P_{L}(\cos\theta)\times\Big{[}V^{(0)}_{L}(R)+\hat{s}_{1}\cdot\hat{s}_{2}V^{(12)}_{L}(R)$
$\displaystyle+\hat{s}_{1}\cdot\hat{s}_{3}V^{(13)}_{L}(R)+\hat{s}_{2}\cdot\hat{s}_{3}V^{(23)}_{L}(R)\Big{]},$
(3)

computed by Gauss-Legendre quadrature including terms up to $L=70$. The
required matrix elements in the channel basis are given in the Appendix.

###  II.2 Nonadditive three-body potential

The nonadditive three-body contribution to the interaction energy in the
triatomic $ABC$ system is computed as

$V^{3b}(R,\theta)=E_{ABC}-\sum_{A<B}E_{AB}+\sum_{A}E_{A}\,,$ (4)

where the total energies of the trimer $E_{ABC}$, dimers $E_{AB}$, and
monomers $E_{A}$ at the trimer geometry $(R,\theta)$ are computed using a
trimer basis set.

First, the total energies and resulting three-body term are computed using the
coupled cluster method Bartlett and Musiał (2007) with the full treatment of
single and double excitations and estimation of the connected triples
contribution non-iteratively by many-body perturbation theory, CCSD(T),
Knowles  _et al._ (1993). Next, the correction to account for the remaining
triple excitations in the coupled-cluster expansion, CCSDT Kállay and Surján
(2001), is added. Large Gaussian basis sets are employed. Thus, all total
energies in Eq. (4) are calculated using

$E=E^{\text{HF}}_{\text{apCV5Z}}+\delta
E^{\text{CCSD(T)}}_{\text{CBS(Q,5)}}+\delta
E^{\text{CCSDT}}_{\text{apCVTZ}}\,,$ (5)

where the Hartree-Fock energy, $E^{\text{HF}}_{\text{apCV5Z}}$, is calculated
in the Douglas-Kroll correlation-consistent polarized core-valence
quintuple-$\zeta$ quality basis sets, aug-cc-pCV5Z-DK Prascher  _et al._
(2011). Next, the correlation energy at the CCSD(T) level, $\delta
E^{\text{CCSD(T)}}_{\text{CBS(Q,5)}}$, is extrapolated to the complete basis
set (CBS) limit with the two-point formula Helgaker  _et al._ (1997) and the
aug-cc-pCVQZ-DK and aug-cc-pCV5Z-DK basis sets. Relativistic effects are
included in those calculations with eXact-2-Component (X2C) Hamiltonian Peng
and Reiher (2012). Finally, the full triples correction $\delta
E^{\text{CCSDT}}_{\text{apCVTZ}}$ defined as a difference between CCSDT and
CCSD(T) results is calculated using the aug-cc-pCVTZ basis sets. We correlate
only the three valence electrons in the CCSDT calculations, which in the case
of the NaLi+Na system is equivalent to reaching the full configuration
interaction (FCI) quality for valence electrons.

The CCSD(T) calculations are performed with the MOLPRO package of ab initio
programs Werner  _et al._ (2012, ), while the CCSDT results are obtained with
the MRCC2019 software Kállay  _et al._ (2020, ). The Legendre expansion of the
nonadditive three-body interaction potential is available in the Supplementary
Material.

The nonadditive three-body interaction potential is computed for the spin-
stretched quartet state where the coupled-cluster method can be applied. We
subsequently assume the nonadditive part of the interaction is spin-
independent and use it for the doublet states, too. This approximation may be
justified by much higher importance of the three-body interaction for the
spin-stretched state where it is larger than the two-body contribution, while
deeply bound doublet states are strongly dominated by two-body interactions.
Additionally, the decomposition into additive and nonadditive parts is
ambiguous for the relevant low-spin states. Finally, later we will show that
the details of doublet states are less important for resonance prediction. We
note that the spin-dependence of the two-body interactions is fully accounted
for.

###  II.3 Electron spin-spin interaction

The electron spin-spin interaction can originate either from direct spin–spin
magnetic dipole–dipole interaction or indirect interaction in the second-order
of perturbation theory mediated by spin-orbit coupling. The total effective
Hamiltonian for the electron spin-spin interaction reads as

$\hat{H}_{SS}=\bm{\hat{S}}^{T}\bm{D}\bm{\hat{S}}$ (6)

where $\bm{\hat{S}}$ represents spin of the system and $\bm{D}$ is $3\times 3$
matrix. In the magnetic axis frame the $\bm{D}$ become diagonal, and the
effective Hamiltonian can be parameterized by $D$ and $E$ constants Jaszuński
_et al._ (2017), as

$\hat{H}_{SS}=D\left[\hat{S}_{z}^{2}-\frac{1}{3}S(S+1)\right]\
+E\left[\hat{S}_{x}^{2}-\hat{S}_{y}^{2}\right]$ (7)

where $\hat{S}_{x}$ is the $x$ component of the spin operator $\hat{\bm{S}}$
with respect to the magnetic axis system, and similar for $y$ and $z$. The
magnetic axes depend on the geometry of the NaLi+Na complex. In the case of
linear configurations, the $z$ axis coincides with the NaLi axis. In the
global minima of the potential energy surface, the $z$ is perpendicular to the
plane defined by atoms, whereas the angle between $y$ and the NaLi axis is
about 30 degrees. Preliminary computations with multi-reference configuration
interaction methods show that the spin-orbit interaction plays a minor role
here, and $E$ is orders of magnitude smaller than $D$. Thus, we neglect the
spin-orbit interaction, include only two-electron spin–spin direct magnetic
interaction. We use the multi-reference averaged quadratic coupled-cluster
(MR-AQCC) Szalay and Bartlett (1993) electronic wave function as implemented
Ganyushin and Neese (2006) in Orca Neese (2012, 2018). We describe scalar
relativistic effects by the Douglass-Kroll-Hess (DKH) Hamiltonian Dyall (2000)
and include the picture-change effects Sandhoefer and Neese (2012). The aug-
cc-pCVTZ-DK basis sets are used. The components of the $\bm{D}$ matrix, the
coupling coefficient $D$, and $E$/$D$ ratio are available in the Supplementary
Material.

We have found that the order of magnitude of the isotropic coupling
coefficient $D$ is conserved for most intermonomer configurations available
during low-energy NaLi+Na collision. In the scattering calculations reported
below, we neglected the interaction-induced variation of the spin-spin
coupling. This means the modification of the spin-spin coupling in a NaLi
molecule due to its interaction with a Na atom is omitted, but the geometry
dependence of the spin-spin coupling in a NaLi molecule on its orientation
within a NaLi+Na complex is included.

###  II.4 Geometry and interpolation

We use Jacobi coordinates to parameterize the geometry of NaLi+Na complex. The
interatomic separation in NaLi is fixed at the vibrationally averaged
interatomic distance of 9.139$~{}a_{0}$ for the ground state of NaLi
(${}^{3}\Sigma^{+}$) Gronowski  _et al._ (2020). We perform the computations
for more than 20 values of the atom-molecule separation $R$. The cosine of
angle ($x_{i}=\cos(\theta_{i}),i\in\\{0,1,2,..,13,14\\}$) are selected from
the roots of the 15-th order Legendre polynomial, and are ordered as
$x_{i}>x_{i+1}$. We independently calculated the nonadditive three-body
contribution [$V^{3b}(R,\theta)$] and electron spin-spin interaction coupling
[$D(R,\theta)$] for each combination of geometric parameters. For each
distance, we interpolate the three-body interaction term and the spin-spin
interaction coupling constants by Legendre polynomials Ogura (1920). We
determine Legendre expansion coefficients using the Gauss–Legendre quadrature
as

$V_{L}^{(3b)}(R)=(2L+1)\sum_{i=0}^{n-1}w_{i}P_{L}(x_{i})V^{3b}(R,\theta_{i}),$
(8)

and similar for the spin-spin coupling, and $w_{i}$ are the quadrature
weights. For use in scattering calculations, the radial dependence of the
expansion coefficients is interpolated using the reproducing kernel Hilbert
space method Ho and Rabitz (1996).

###  II.5 Accuracy

The global minimum of the molecule-atom interaction potential for the quartet
state of NaLi+Na with the internuclear distance in NaLi fixed at
9.139$\,a_{0}$ occurs at $R_{e}=7.48\,a_{0}$ and
$\theta_{e}=117.3\,^{\text{o}}$ with the well depth of $D_{e}=812.5\,$cm-1
(see Fig. 1). At this equilibrium geometry two-body and three-body
interactions contribute 230.7 cm-1 and 581.8 cm-1 to the binding energy,
respectively. The CCSD(T) method reproduced 555.0 cm-1 (95%) of the
nonadditive part and the CCSDT correction accounts for 26.8 cm-1. The final
uncertainty of the calculated nonadditive three-body potential is about 17
cm-1 (3$\,\%$) at the global minimum of the potential energy surface and
results mainly from the incompleteness of the basis set in the CCSD(T)
($\sim$5 cm-1) and CCSDT ($\sim$8 cm-1) computations, and the lack of core-
electron correlation at the CCSDT level ($\sim$3 cm-1). We estimate the
magnitude of the neglected contribution of core and core-valence correlation
beyond CCSD(T) level by running all-electron CCSDT calculations with the
virtual one-electron space truncated by 20 % Rolik and Kállay (2011).

Our calculations show that the three-body term is larger than the two-body one
and predominantly attractive. Consequently, the intermonomer NaLi distance is
shorter than the intramolecular one. The interatomic distances in the complex
are closer to the classical turning point of diatomic potentials than to their
minima. The repulsive part of the two-body potential is less accurate than its
long-range part. This introduces additional uncertainty to the overall
potential, which can be as high as 16 cm-1 (an additional uncertainty of
$\sim$3 cm-1 of the three-body interaction), estimated from the difference
between an experimental two-body potential for Na2 Bauer and Toennies (2019)
and a highly accurate theoretical potential for NaLi Gronowski  _et al._
(2020).

The native uncertainty of our state-of-the-art electronic calculation of the
three-body term is relatively small, although still too large to predict
molecular scattering lengths. Additional uncertainty is introduced by using
the rigid rotor approximation to describe NaLi within NaLi+Na. This
approximation works very well for van der Waals complexes of deeply-bound
molecules, but for high-spin weakly bound alkali-metal systems may result in
non-negligible errors Tomza  _et al._ (2013). Our calculations show that the
global equilibrium geometry for NaLi+Na with relaxed distance in NaLi has
isosceles triangular symmetry and internuclear distances smaller by around 20%
than in the diatomic molecules due to large three-body forces, the magnitude
of which monotonically increases with decreasing atom-molecule separation. To
reflect the errors related to the difference in three-body interaction at the
global minimum and the minimum within the rigid rotor approximation used, as
well as the increased uncertainty of the two-body interactions at these
shorter interatomic distances, we used a more conservative estimate of the
uncertainty of the three-body term in the remainder of this paper.

##  III Scattering calculations

###  III.1 Basis set and Hamiltonian

In our coupled-channels calculations the scattering wavefunction is expanded
in a basis of fully coupled channel functions of the form

$\displaystyle|(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}\rangle=\sum_{M_{J},M_{S}}\langle
JM_{J}SM_{S}|\mathcal{J}\mathcal{M}\rangle$ $\displaystyle\times$
$\displaystyle|(NL)JM_{J}\rangle|(s~{}s_{3})SM_{S}\rangle$ (9)

where $\langle j_{1}m_{1}j_{2}m_{2}|jm\rangle$ is a Clebsch-Gordan
coefficient. The quantum number $N$ represents the rotational angular momentum
of the NaLi molecule, and $L$ the angular momentum associated with the end-
over-end rotation of the atom and molecule about one another, which are
Clebsch-Gordan coupled to a total mechanical angular momentum $J$ with
$B$-field projection $M_{J}$. Similarly, $s$ denotes the molecular electronic
spin resultant from coupling atomic spins $s_{1}$ and $s_{2}$ within NaLi,
whereas $s_{3}$ is the electronic spin of Na, and $S$ the total electronic
spin with $B$-field projection $M_{S}$. In the coupled basis, $J$ and $S$ are
subsequently coupled to a total angular momentum $\mathcal{J}$ with magnetic-
field projection $\mathcal{M}=M_{J}+M_{S}$. Nuclear spin is initially not
taken into account, see Sec. III.7. We assume the molecular bond length fixed
at the triplet equilibrium position.

The total angular momentum projection along the magnetic field axis
$\mathcal{M}=M_{J}+M_{S}$ is strictly conserved. For large enough magnetic
field $M_{S}$ becomes a good quantum number, and therefore also
$M_{J}=\mathcal{M}-M_{S}$ is good. The Na atomic spin is $s_{3}=1/2$
throughout. Due to the large singlet-triplet splitting in the NaLi molecule,
$s=0$ or 1 is also a good quantum number. For a separated atom and molecule,
$m_{s}$ and $m_{s_{3}}$ would separately become good quantum numbers, but at
chemically relevant distances the exchange splitting between the doublet and
quartet interaction potentials is dominant such that $S=1/2$ and $3/2$ are
good quantum numbers. Hence, we can effectively consider each $|SM_{S}\rangle$
state separately, with only perturbatively weak couplings between them. For
each of these spin channels, there are strong and anisotropic interactions
that couple different $N$ and $L$ channels, but conserve $J$ and $M_{J}$. The
initial channel corresponds to $s$-wave collisions in the spin-stretched
rotational ground state, $|(NL)JM_{J}\rangle|SM_{S}\rangle=|(00)00\rangle|3/2\
3/2\rangle$ for the low-field seeking upper spin-stretched state, and
$|(00)00\rangle|3/2,\ -3/2\rangle$ for the high-field seeking lower spin-
stretched state, respectively.

In our coupled-channels calculations we include the electronic interaction
described in Sec. II,

$\displaystyle\hat{V}(R,\theta)$
$\displaystyle=\sum_{L}P_{L}(\cos\theta)\times\Big{[}V^{(0)}_{L}(R)+\hat{s}_{1}\cdot\hat{s}_{2}V^{(12)}_{L}(R)$
$\displaystyle+$
$\displaystyle\hat{s}_{1}\cdot\hat{s}_{3}V^{(13)}_{L}(R)+\hat{s}_{2}\cdot\hat{s}_{3}V^{(23)}_{L}(R)+V_{L}^{(3b)}(R)\Big{]},$
(10)

the Zeeman interaction with the magnetic field,

$\displaystyle\hat{H}_{\mathrm{Zeeman}}=\mu_{B}g_{e}B(\hat{s}_{1,z}+\hat{s}_{2,z}+\hat{s}_{3,z}),$
(11)

the magnetic dipole-dipole interaction

$\displaystyle\hat{V}_{\mathrm{magn.dip}}=$
$\displaystyle-\sqrt{30}\frac{(\mu_{B}g_{e}\alpha)^{2}}{R^{3}}$
$\displaystyle\times\left[\left[\hat{s}\otimes\hat{s}_{3}\right]^{(2)}\otimes
C^{(2)}(\hat{R})\right]^{(0)}_{0},$ (12)

and the spin-rotation

$\displaystyle\hat{H}_{\mathrm{spin-rotation}}=\gamma_{s}\hat{N}\cdot\hat{s},$
(13)

and spin-spin couplings

$\displaystyle\hat{H}_{\mathrm{spin-
spin}}=\lambda_{s}\sqrt{30}/3\left[\left[\hat{s}\otimes\hat{s}\right]^{(2)}\otimes
C^{(2)}(\hat{r}_{\mathrm{NaLi}})\right]^{(0)}_{0},$ (14)

where $[\hat{A}\otimes\hat{B}]^{(k)}_{q}$ indicates a tensor product, see the
Appendix. The spin-spin coupling parameter $\lambda_{s}=0.0189$ cm-1 is
computed here, whereas for the spin-rotation coupling we use the
$\gamma_{s}=0.005$ cm-1 upper limit estimated in Ref. Gronowski  _et al._
(2020). The reader is referred to the Appendix for a full description of the
Hamiltonian and for matrix elements in the channel basis.

The spin-rotation coupling $\gamma_{s}\hat{N}\cdot\hat{s}$ and spin-spin
coupling
$\lambda_{s}\left[\left(\hat{s}\cdot\hat{r}_{\mathrm{NaLi}}\right)\left(\hat{s}\cdot\hat{r}_{\mathrm{NaLi}}\right)-\frac{1}{3}\hat{s}^{2}\right]$
play an important role here as these are spin-dependent couplings that are not
diagonal in $S$ and $M_{S}$. Losses by Zeeman relaxation or transitions to the
chemically reactive doublet potential must involve these couplings. It is
worth noting that if the electronic interaction were isotropic, _i.e._ ,
independent of the relative orientation of the atom and molecule, these spin-
dependent couplings would not be enough to lead to Zeeman relaxation. However,
after accounting for the spin-rotation and spin-spin coupling, molecular
eigenstates in different Zeeman levels have different rotational-state
decompositions, and hence are coupled by the anisotropic part of the
electronic interaction. That is, physically, the strong anisotropic electronic
interaction can reorient the molecule which effectively flips the electronic
spin because the spin is coupled to the molecular axis by spin-rotation and
spin-spin coupling. The role of anisotropic interactions implies that
transitions occur at short range and cannot be described by simpler models
based on isotropic long-range $R^{-6}$ interactions alone Wang  _et al._
(2021).

After having described the atom-molecule interactions, approximately good
quantum numbers, and the critical role of various coupling mechanisms, we
continue by performing full quantum mechanical coupled-channels calculations.
For computational tractability we will start out ignoring hyperfine and
vibrational degrees of freedom. By ignoring the vibrational coordinate of the
molecule, and fixing the molecular bond length to the triplet NaLi equilibrium
bond length, our model cannot directly describe chemical reactions. The only
energetically accessible products are to form singlet NaLi or Na2 molecules.
Chemical reactions will occur on the low-spin potential. In our coupled-
channels calculations, we model these by imposing an absorbing boundary
condition at $R_{\mathrm{min}}=4.5~{}a_{0}$, which can be reached on the low-
spin potential, but not on the high-spin potential which is highly repulsive
at these short distances, see Fig. 1. This choice seems arbitrary but it does
not affect the results as long as the boundary condition is imposed in a
region where the high-spin potential is highly repulsive and simultaneously
the low-spin potential is strongly attractive. Any flux that reaches this
region in the low-spin state will continue classically to smaller $R$, such
that it does not matter where exactly in this region one matches to the
absorbing boundary condition. We have confirmed this numerically for
$R_{\mathrm{min}}$ between 4 and 4.5 $a_{0}$.

Figure 1:  Interaction potentials for the NaLi+Na collision complex. The
potentials are shown as a function of $R$, the center of mass distance between
the atom and molecule, and $\theta$ the Jacobi angle between the orientation
of the molecule and the approach of the atom, with $\theta=0$ corresponding to
Na-NaLi and $\theta=180^{\circ}$ to Na-LiNa. The molecular bond length is
fixed at the triplet equilibrium bond length $r=8.8~{}a_{0}$. (a) high-spin
quartet $S=3/2$ potential. (b) low-spin doublet $S=1/2$ potential. In
principle there are two doublet potentials that can have an avoided crossing.
What is shown here is the doublet potential for a pure triplet NaLi molecule,
$s=1$.

###  III.2 Cross sections and rate coefficients

We solve the coupled-channels equations numerically using the renormalized
Numerov propagator. Using the method of Janssen  _et al._ (2013), and
described in more detail in Ref. Karman , we match to reactive boundary
conditions at short range and the usual scattering boundary conditions at long
range. Again, the short-range boundary condition is imposed at
$R_{\mathrm{min}}=4.5~{}a_{0}$, a separation that can be reached only on the
low-spin potential and effectively models chemical reactions. This procedure
yields an “inelastic” $S$-matrix, $\bm{S}^{(\mathrm{LR})}$ and a “reactive”
$S$-matrix $\bm{S}^{(\mathrm{SR})}$. The elements of the inelastic $S$-matrix,
${S}^{(\mathrm{LR})}_{f,L^{\prime};i,L}$, describe the amplitudes for
scattering from an initial state $i$ and partial wave $L$, to a final state
$f$ and partial wave $L^{\prime}$. The elements of the reactive $S$-matrix,
${S}^{(\mathrm{SR})}_{r;i,L}$ describe scattering from an initial state $i$
and partial wave $L$ into a reactive channel $r$ at short range. The reactive
channels are determined by diagonalizing the Hamiltonian excluding radial
kinetic energy at the short-range matching point. From the $S$-matrices one
can determine the elastic cross section

$\displaystyle\sigma^{\mathrm{elastic}}=\frac{\pi}{k^{2}}\sum_{L,L^{\prime}}\left|\delta_{L,L^{\prime}}-S^{(\mathrm{LR})}_{i,L^{\prime};i,L}\right|^{2},$
(15)

where $i$ is the initial state and $k=\hbar^{-1}\sqrt{2\mu E}$ is the intial
wavenumber. From the $S$ matrices one can also determine the inelastic rate
coefficient

$\displaystyle k^{(\mathrm{inelastic})}=\frac{\pi}{\mu
k}\sum_{L,L^{\prime},f\neq
i}\left|S^{(\mathrm{LR})}_{f,L^{\prime};i,L}\right|^{2},$ (16)

the rate coefficient for loss at short range

$\displaystyle k^{(\mathrm{short})}=\frac{\pi}{\mu
k}\sum_{L,r}\left|S^{(\mathrm{SR})}_{r;i,L}\right|^{2},$ (17)

and we define a total loss rate coefficient,
$k^{(\mathrm{loss})}=k^{(\mathrm{inelastic})}+k^{(\mathrm{short})}$. Both the
elastic cross section and the rate coefficients are energy independent for
energies well below the van der Waals energy which is in the order of 500
$\mu$K for NaLi+Na. We computed the cross sections and loss rate coefficients
at a single collision energy of 4 $\mu$K.

###  III.3 Convergence

First we consider convergence of the calculation with $J_{\mathrm{max}}$, the
highest value of $J$ included in the basis set. As explained above, we can
essentially consider each spin channel $|SM_{S}\rangle$ independently with
only perturbative couplings between them, and within each spin level, we can
consider $J$ and $M_{J}$ to be good quantum numbers. Spin-rotation coupling
has the selection rules $\Delta J=0\rightarrow 1$ and $|\Delta M_{S}|\leq 1$.
Spin-spin coupling has the selection rules $\Delta J=0\rightarrow 2$ and
$|\Delta M_{S}|\leq 2$. Hence, if these spin-dependent couplings act only
perturbatively we expect cross sections do not change by including functions
with $J=3$ or higher, and that cross sections scale quadratically in the spin-
rotation and spin-spin coupling strengths. This is exactly what we observe in
Fig. 2(a). Inclusion of functions with $J=3$ does not affect the cross
sections (dots), and the cross sections scale perturbatively with the spin-
rotation and spin-spin coupling strengths (crosses). The contribution of the
magnetic dipolar interaction between the atomic and molecular magnetic moment
is far smaller and we have discussed in less detail for this reason, although
it is included in the calculation. Finally we note that the mechanism
responsible for the loss rates requires anisotropic electronic interactions,
but the anisotropy is not perturbatively weak.

The convergence of the calculation with $N_{\mathrm{max}}$ for fixed
$J_{\mathrm{max}}=1$ is shown in Fig. 2(b). The expectation is that the
calculation converges when all locally open channels are included, _i.e._ ,
when the excitation energy of channels that are excluded are all higher than
the depth of the interaction potential. Morita  _et al._ (2018, 2019, 2020) If
we estimate the required $N_{\mathrm{max}}$ based on the depth of the spin-
stretch potential of 800 cm-1, one would expect the calculation to converge
with $N_{\mathrm{max}}\approx 70$. Instead, we see the calculation converges
with much higher $N_{\mathrm{max}}=350$. The reason for this appears to be
that these higher rotational channels continue to contribute locally open
channels in the low-spin state that affect the background scattering length.

As we have seen, converging the scattering calculation requires including also
functions with $J=2$. Scattering calculations for basis sets truncated at
$N_{\mathrm{max}}=350$ and $J_{\mathrm{max}}=2$ become prohibitively computer
intensive. Fortunately, the typical number of resonances appears to converge
more rapidly with $N_{\mathrm{max}}$, and can be predicted with much lower
truncation of $N$. Figure 2(c) shows loss rates computed with various
$N_{\mathrm{max}}$ up to 30 for fixed $J_{\mathrm{max}}=2$. We will come back
to convergence of the typical density of resonances after discussing the
dependence on the interaction potential.

Figure 2:  Convergence of the coupled-channels calculations with truncation of
the basis set. (a) Calculation with fixed $N_{\mathrm{max}}=20$. Cross
sections converge with $J_{\mathrm{max}}=2$. When the basis is truncated at
$J_{\mathrm{max}}=1$, the cross sections scale quadratically with the spin-
rotation coupling constant $\gamma_{s}$. This is demonstrated by agreement
with the crosses which show one-quarter of the cross section obtained by
scaling $\gamma_{s}$ up by a factor of two. When the basis is truncated at
$J_{\mathrm{max}}=2$, spin-spin coupling contributes and the cross section
scales perturbatively with both coupling constants (crosses). Spin-spin
coupling is typically dominant, but not by a large factor so both mechanisms
contribute. Magnetic dipole-dipole coupling does not play an important role.
(b) Calculation with fixed $J_{\mathrm{max}}=1$ The scattering cross section
can be converged with $N_{\mathrm{max}}$, but requires an impractically large
basis set. (c) Calculation with fixed $J_{\mathrm{max}}=2$ The scattering
cross section is not converged with $N_{\mathrm{max}}$, but the density and
typical width of resonances are not strongly dependent on $N_{\mathrm{max}}$.
Figure 3:  Dependence of loss rates on scaling of three-body interactions.
Loss rates as a function of $B$-field and scaling of the three-body
interaction by a factor $1+\lambda$ for $N_{\mathrm{max}}=20$, shown for the
upper stretched state. This shows that due to the uncertainty of $\lambda$ of
at least several percent, the background loss and resonance positions are
undetermined. Several $B$-independent resonances are observed, where the
$\lambda$-scaling tunes the initial spin-stretched potential such that it
supports a resonance, i.e., a bound state near zero energy. Hence, prediction
of resonance positions requires knowledge of the interaction potentials to an
accuracy that cannot realistically be achieved by _ab initio_ calculations.

Figure 4:  Spin-dependence of three-body interactions. Loss rates are shown as
a function of $B$-field and scalings of the interactions potential.
Calculations are done with $N_{\mathrm{max}}=6$, which artificially reduced
the number of resonances compared to Fig. 3. (a) The spin-independent three-
body interaction is scaled by $1+\lambda$. In this case, all resonance
positions depend on the scaling of the potential, the background loss rate
varies strongly, and new $B$-independent resonances appear by scaling the
initial-state potential such that supports a bound state near zero energy. (b)
The spin splitting is scaled by $1+\lambda$, while the spin-stretched
potential is kept fixed. Therefore, only the doublet potential is varied. The
resonances near 150, 350, 1150, and 1250 G are independent of the scaling of
the low-spin potentials, i.e., they are completely supported by the non-
reactive spin-stretched potential. The feature just above 1000 G has a weak
dependence on the scaling of the doublet potentials, and several broader
features such as that around 500 G have a stronger dependence.

Figure 3 shows the sensitivity of the scattering rates to the interaction
potential, parameterized by $\lambda$. Here, we scale by a factor $1+\lambda$
the non-additive three-body part of the interaction potential, _i.e._ , the
part that is computed _ab initio_ and is uncertain up to an estimated 3 %
within the rigid rotor approximation for NaLi, and considerably larger when
this approximation is relaxed, see Sec. II.5. Figure 3 is computed with
$J_{\mathrm{max}}=2$ and $N_{\mathrm{max}}=20$. By modifying the three-body
interaction by less than $10^{-4}$, we find that the resonances are
unaffected. By modifying the potential by about $10^{-3}$, we find that the
resonances start to shift such that realistically their positions are
completely undetermined. When the scaling of the potential is at the level of
several percents, we find that horizontal $B$-field independent resonances
appear. This occurs when the potential is scaled to support a resonance near
zero collision energy for the _initial_ Zeeman level, which is therefore not
tuned by the magnetic field. Therefore, with realistic uncertainties in the
three-body part of the interaction and the rigid rotor approximation, both the
position of magnetically tunable resonances and the background scattering
length are completely undetermined but we can still draw probabilistic
conclusions about a typical range of values for the scattering length.

Figure 4 shows a similar “$\lambda$-scan”, but now for $N_{\mathrm{max}}=6$
which artificially reduces the density of resonances somewhat and produces a
more sparse figure. In panel 4(a) we have scaled the non-additive three-body
part of the interaction potential for both the high-spin and the low-spin
potential, as before. In panel 4(b), by contrast, we have assumed this
uncertainty is entirely in the low-spin potential, and we have kept the high-
spin potential constant. In this case, we no longer observe horizontal
$B$-field resonances as these are supported by the initial top Zeeman energy
level, which is a high-spin state. We find that several of the resonances are
now completely independent of the scaling of the low-spin potential up to
$\lambda=0.1$. This indicates the resonances are supported by the non-reactive
high-spin potential. When the scaling of the low-spin potential reaches the
percent level, coupling to the low-spin state starts to affect the collision
rate, especially for the broader features.

The analysis above shows that the resonances are supported by the high-spin
potential, and that their positions are sensitive to 0.1 % uncertainty in the
non-additive three-body part of the interaction potential. This means that the
_ab initio_ prediction of the resonance positions is beyond the capability of
state-of-the-art theory. In addition to highly accurate atom-molecule
interactions, the quantitative prediction of resonance positions would require
converged coupled-channels calculations that also fully account for hyperfine,
molecular vibrations, and chemical reactions. This is not attempted here, and
instead we interpret only the typical number of resonances, their widths, and
coupling mechanisms.

###  III.4 Background elastic and inelastic scattering

Next we consider interaction potentials obtained by scaling the non-additive
three-body interactions by $1+\lambda$ with $\lambda$ between $-0.1$ and
$+0.1$, on the order of the uncertainty in the interaction potential. We
consider each of these Hamiltonians statistically independent realizations of
the physical NaLi+Na system. By performing scattering calculations for these
different realizations, we gather statistics about the number of resonances
and their widths.

Figure 5:  Cross sections and rates as a function of $\lambda$ scaling. (a)
The elastic collision cross section obtained by scaling the non-additive
three-body interaction by $1+\lambda$ at a fixed magnetic field $B=1500$ G.
Different colors correspond to different spin states,
$|sm_{s}\rangle|s_{3}m_{s_{3}}\rangle$, as indicated. (b) Loss rate
coefficients as a function $\lambda$ scaling for fixed $B=1500~{}$G.

Figure 6:  Probability distribution of collision rates. Cumulative probability
distribution obtained by sorting cross sections and rate coefficients for
different $\lambda$ scaling at fixed $B=1500$ G Morita  _et al._ (2019).
Different colors correspond to different states,
$|sm_{s}\rangle|s_{3}m_{s_{3}}\rangle$, as indicated. (a) Distribution of the
elastic cross section, (b) loss rate coefficient, and (c) the ratio of
elastic-to-inelastic collisions. Solid (dashed) lines correspond to
$N_{\mathrm{max}}=20$ ($10$). In panel(a) the black markers indicated the
expected distribution of elastic cross sections for a van der Waals potential
Gribakin and Flambaum (1993).

We first consider the nonresonant background. Figure 5 shows the dependence of
the elastic cross section and loss rate coefficient on scaling of the non-
additive three-body interaction by $1+\lambda$ at a fixed magnetic field
$B=1500$ G, $N_{\mathrm{max}}=20$. These calculations are performed for three
initial states; the top stretched state where the molecular electronic spin
projection $m_{s}=1$ and the atomic electronic spin projection $m_{3}=1/2$,
the bottom stretched state where $m_{s}=-1$ and $m_{3}=-1/2$, and the non-
stretched state $m_{s}=1$, $m_{3}=-1/2$. The elastic cross sections are nearly
identical in the two stretched states, and different but of the same order in
the non-stretched initial state. For the loss rate coefficient, the _typical_
behavior is also that the cross sections are similar between the two stretched
states, while the differences at fixed $\lambda$ can be as large as an order
of magnitude. For the non-stretched initial state the loss rate coefficient is
significantly larger. This can be understood as follows: For non-stretched
states the collision can directly proceed on the reactive potential and lead
to loss. For stretched states, the initial potential is nonreactive such that
loss processes require a spin flip which is unlikely since it is perturbative
in the weak spin-spin and spin-rotation coupling, as we have seen above.

In Figure 6 we show again the effect of scaling of the non-additive three-body
interaction by $1+\lambda$, but ordered by increasing cross section and loss
rate. We interpret the resulting horizontal axis as the cumulative probability
to obtain a specific cross section or rate coefficient Morita  _et al._
(2019). Panel 6(a) shows that elastic cross sections, for which all cumulative
probability distributions are similar and agree closely with the expected
distribution for a $R^{-6}$ potential Gribakin and Flambaum (1993), shown as
the black markers. This distribution is characterized by a mean scattering
length $\bar{a}=0.47799(2\mu C_{6}/\hbar)^{1/4}\approx 52~{}a_{0}$. The
experimentally measured scattering length is 260 $a_{0}$, the corresponding
elastic cross section is larger than one might expect.

Panel 6(b) shows the cumulative probability distribution of the loss rate
coefficient. The two spin-stretched states show a similar distribution,
whereas the typical loss rate is much higher in the non-stretched state, as
observed before. Similar probability distributions are obtained with
$N_{\mathrm{max}}=10$ and $20$. In the spin-stetched states, the background
loss rate is likely to lie between $10^{-13}$ and a few $10^{-12}$ cm3/s.
Experimentally, a typical background loss rate of $5\times 10^{-12}$ cm3/s is
observed, which on the higher end of this distribution function, but not in
disagreement with it.

Finally, in panel 6(c), we consider the background ratio of elastic-to-
inelastic collisions, $\gamma$, We compute the ratio of elastic-to-inelastic
collisions as

$\displaystyle\gamma=\frac{\sigma^{(\mathrm{elastic})}\langle
v\rangle}{k^{(\mathrm{loss})}},$ (18)

where $\langle v\rangle$ is the thermally averaged velocity at $T=2~{}\mu$K
Son  _et al._ (2020), and we take for the elastic cross section $4\pi a^{2}$
where $a=260~{}a_{0}$ is the measured elastic cross section. The resulting
ratio of elastic-to-inelastic collisions is likely between 100 and 1000 in the
spin-stretched states, in agreement with the experimental value of 300 Son
_et al._ (2020).

Figure 7:  Representative magnetic field scans for different basis set
truncation. The figure shows representative magnetic field scans of the
collisional loss rate coefficient. Different panels correspond to truncation
of the basis set at $N_{\mathrm{max}}=2,$ 10, 30, as indicated. The left
(right) hand column shows results for the upper (lower) spin-stretched state.

###  III.5 Density of resonances

Next, we compare the typical number of resonances to the density of states of
the spin-stretched NaLi+Na collision complex. Figure 7 shows typical magnetic
field scans for $N_{\mathrm{max}}=2$, 10, and 30 for both the top and bottom
spin-stretched states. The density of resonances increases with
$N_{\mathrm{max}}$ and is higher for the lower spin-stretched state than it is
for the upper one.

In order to explain the observed number of resonances, we calculate their
density of states. Again, the resonances are supported by the spin-stretched
interaction potential, and the magnetically tunable resonances correspond to
Zeeman sub-states that are different from the initial state. This means each
resonance can be assigned total electron spin $S$, $M_{S}$ quantum numbers.
Spin is coupled to the spatial degrees of freedom perturbatively through spin-
spin and spin-rotation coupling. This means we can assign each state total
mechanical angular momentum $J$, $M_{J}$, and the spin-rotation and spin-spin
selection rules tell us which values of these nearly-good quantum numbers
contribute. Spin-rotation coupling is rank-1 in the spin and spatial degrees
of freedom and couples to states with $J=1$ and $|\Delta M_{S}|\leq 1$. The
density of such states is determined by the $J=1$ density of states on the
spin-stretched potential. The total number of bound states below energy $E$
can be computed quasiclassically using a phase-space integral Christianen  _et
al._ (2019a). We determine the total number of resonances expected as the
total number of bound states less than the Zeeman energy shift for $\Delta
M_{S}=1$ below threshold. Spin-spin coupling is rank-2 in the spin and spatial
degrees of freedom and couples to $J=2$ and $|\Delta M_{S}|\leq 2$ states. We
similarly determine the number of resonances by computing the number of $J=2$
states bound by less than the Zeeman shifts for $\Delta M_{S}=1$ and $2$.

We computed the total numbers of bound states using the phase-space integrals
of Ref. Christianen  _et al._ (2019a),

$\displaystyle
N^{(\mathrm{3D})}=(2J+1)\frac{g_{\mathrm{parity}}8\pi\sqrt{2}m_{\mathrm{Na}}^{2}m_{\mathrm{Li}}}{3h^{3}(2m_{\mathrm{Na}}+m_{\mathrm{Li}})}$
$\displaystyle\times\iiint\frac{rR}{\sqrt{\mu
R^{2}+\mu_{\mathrm{NaLi}}r^{2}}}\left[E-V(\bm{q})\right]^{3/2}\ dr\ dR\
d\theta$ (19)

as well as the number of bound states for the NaLi vibrational coordinate
$r=r_{e}$ fixed at the equilibrium distance

$\displaystyle N^{(\mathrm{2D})}=(2J+1)\frac{g_{\mathrm{parity}}2\pi
m_{\mathrm{Na}}^{2}m_{\mathrm{Li}}}{h^{2}(2m_{\mathrm{Na}}+m_{\mathrm{Li}})\sqrt{\mu_{\mathrm{NaLi}}}}$
$\displaystyle\times\iint\frac{rR}{\sqrt{\mu
R^{2}+\mu_{\mathrm{NaLi}}r^{2}}}\left[E-V(\bm{q})\right]\ dR\ d\theta,$ (20)

where $\bm{q}=\\{R,r,\theta\\}$ are the Jacobi coordinates, $V(\bm{q})$ is the
interaction potential, $m_{x}$ are the atomic masses, and
$g_{\mathrm{parity}}=1/2$ is a factor that accounts for parity conservation.
From this we find that we should expect to encounter approximately 14 Feshbach
resonances when excluding the vibrational coordinate. This is a useful
reference for the scattering calculations, where the vibrational coordinate is
fixed. When vibrations are accounted for, the number of resonances increases
from 14 to 21. Hence, the NaLi molecular vibrations contribute significantly
to the density of states of the complex, but not by orders of magnitude. Hence
a large fraction of the resonances corresponds to NaLi in the vibrational
ground state, justifying freezing the vibrational coordinate in our scattering
calculations. Figure 8(a) shows the quasiclassical density of resonances as a
function of $R_{\mathrm{max}}$, the upper integration limit used in evaluating
Eqs. (19) and (20). This shows that the resonances are supported by atom-
molecule distances up to approximately 40 $a_{0}$, which is considerably
shorter than the range of van der Waals interactions $R_{6}=(2\mu
C_{6})^{1/4}\approx 109~{}a_{0}$ in the system, suggesting that in some sense
the resonances are of short-range nature. Nevertheless, at atom-molecule
distances around 40 $a_{0}$ the electronic interaction is close to its
$R^{-6}$ asymptotic form, and the vast majority of the density of states is
hence supported by the van der Waals potential. This was previously also
argued for NaK+K collisions Wang  _et al._ (2021), and may more generally be
true of atom+molecule collisions Frye and Hutson (2022). Finally, we will see
that the coupling mechanism for these resonances involves the anisotropy of
the atom+molecule interaction, and hence this cannot be understood in terms of
the long-range interaction which is isotropic for a molecule in its rotational
ground state.

Because the density of states increases with $2J+1$, and because for $J=2$
resonances can occur with both $|\Delta M_{S}|=1$ and $2$, the typical density
of resonances for $J=2$ is a factor five larger than for $J=1$. Hence we can
conclude that most of the resonances observed, approximately 5 out of 6
resonances, can be assigned $J=2$ and are due to spin-spin coupling.

Figure 8:  Number of resonances.  The plot shows the expected number of
resonances between 0 and 1500 G from ($\mathbf{a}$) quasiclassical phase-
integrals Eqs. (19) and (20) as a function of the maximum molecule-atom
center-of-mass separation $R_{\mathrm{max}}$ up to which the integrals are
evaluated and ($\mathbf{b}$) quantum mechanical calculations of bound states
on each adiabatic potential energy curve, as a function of $N_{\mathrm{max}}$
that truncates the channel basis set.

We also compute the density of states quantum mechanically using the same
channel basis as used in the scattering calculations. This is useful for a
direct comparison to the scattering calculations. To this end, we first
compute adiabatic potential energy curves by diagonalizing the Hamiltonian
excluding radial kinetic energy at every grid point. On each adiabatic
potential curve, we compute bound state wavefunctions using sinc-function
discrete variable representation Colbert and Miller (1992). We record the
number of bound states below the initial state threshold. We repeat this at
both ends of the magnetic field range, and the difference in the number of
bound states provides an estimate for the number of resonances. The results
are shown in Fig. 8(b). As can be seen, the total number of resonances
increases with $N_{\mathrm{max}}$ truncation of the basis set, and converges
around $N_{\mathrm{max}}=30$. At $N_{\mathrm{max}}=30$ the highest adiabatic
potentials for the spin-stretched state no longer support bound states, thus
the typical number of resonances converges, although the background scattering
length does not converge until much higher $N_{\mathrm{max}}$ due to the
stronger interactions in the low-spin states. The total number of resonances
is close to the quasiclassical estimate, but not in perfect agreement with it
due to the light masses and relatively weak interactions in the spin-stretched
state. The total number of resonances matches observations for a typical
$B$-field scan in the lower spin-stretched state. However, in the upper spin-
stretched state we find a similar density of states, although typically a much
lower number of resonances is observed.

Figure 9:  Representative magnetic field scans with non-initial Zeeman states
removed. Shown are representative magnetic field scans obtained by excluding
the non-initial Zeeman states for $N=0$ from the calculation. Different panels
correspond to truncation of the basis set at $N_{\mathrm{max}}=2,$ 10, 30, as
indicated. The left (right) hand column shows results for the upper (lower)
spin-stretched state. Excluding the lower-lying Zeeman levels for the upper
spin-stretched state increases the typical number of resonances, whereas
excluding excited Zeeman levels for the lower spin-stretched state does not
reduce the typical number of resonances.

###  III.6 Missing resonances in the upper spin-stretched state

One might expect that the higher number of resonances in the lower spin-
stretched state comes from quasi-bound states in the Zeeman-excited $N=0$
states. However, this would be reflected in a higher density of states for the
lower stretched state in the calculation above. What happens is not that the
total number of states is higher for the lower state, but rather that we do
not observe every state as a resonance in the upper spin-stretched state. This
is caused by fast decay of some of the resonances to the lower-lying Zeeman
levels. This is illustrated in Fig. 9, which shows typical $B$-field scans
from a calculation that excludes non-initial Zeeman states in $N=0$. In the
case of the lower stretched state, the typical number of resonances is not
reduced by omitting the excited Zeeman levels, _i.e._ , these do not cause the
higher number of resonances. In the case of the upper stretched state, the
typical number of resonances is increased by _omitting_ the lower-lying Zeeman
states. Some of these resonances were previously not visible due to fast loss
to the lower-lying Zeeman states.

Loss to lower-lying Zeeman states WITH N=0 can occur only via spin-spin
coupling. Spin-rotation coupling cannot simultaneously fulfil the parity
selection rule and the selection rule $\Delta J=1$. The argument is as
follows. In the rotational ground state $N=0$, the orbital and total
mechanical angular momentum are equal, $J=L$. For $J=1$ this means that all
levels in the rotational ground state have $L=1$, and hence odd parity. Since
the parity of $N+L$ is conserved, these inelastic exit channels are
inaccessible. Note that spin-rotation coupling can lead to resonances and
chemical loss since in these cases, N is changing.

Figure 10:  Representative magnetic field scans with the channel basis
truncated with $J_{\mathrm{max}}=1$. Different panels correspond to truncation
of the basis set at $N_{\mathrm{max}}=2$, 10, 30, as indicated. The left
(right) hand column shows results for the upper (lower) spin-stretched state.
Truncation of the basis set with $J_{\mathrm{max}}=1$ has reduced the number
of resonances, and the qualitative difference in the number of visible
resonances between the two spin states has disappeared.

Figure 10 shows once more typical $B$-field scans but now including $J=0,1$
states only. This has artificially removed the contribution of the spin-spin
interaction, which couples to $J=2$ states. Therefore, no decay to lower
Zeeman states WITH N=0 is possible. This removes the effect discussed
previously where fast inelastic scattering to lower Zeeman states renders
fewer resonances observable in the upper spin-stretched state. The qualitative
differences between the upper and lower spin-stretched states occur only for
$J=2$ resonances, are attributed to spin-spin coupling and would therefore not
occur for collisions between atoms and spin-doublet molecules.

###  III.7 Hyperfine interactions

Figure 11:  Including hyperfine interactions in the calcuations.
Representative magnetic field scans are obtained with and without hyperfine
interactions. Panel (a) shows results for the upper spin-stretched state and
panel (b) the lower spin-stretched state. Including hyperfine interactions
increases the number of resonances somewhat, but not significantly, and it
does not lead to clearly identifiable multiplet splittings.

Next, we investigate the effect hyperfine interactions by including nuclear
spin in the coupled-channels scattering calculations. To this end, the basis
set of Eq. (9) is extended to functions of the form

$\displaystyle|(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}\rangle|i_{1}m_{i_{1}}\rangle|i_{2}m_{i_{2}}\rangle|i_{3}m_{i_{3}}\rangle.$
(21)

The Na nuclear spin is $i_{1}=i_{3}=3/2$ and the Li nuclear spin is $i_{2}=1$.
Only functions with $\mathcal{M}+m_{i_{1}}+m_{i_{2}}+m_{i_{3}}$ equal to
$+11/2$ are included for the upper spin-stretched state, and $-11/2$ for the
lower spin-stretched state, respectively. The hyperfine couplings included
take the form $a\ \hat{i}\cdot\hat{s}$ for each atom. We write this dot
product as
${\hat{i}_{z}\hat{s}_{z}-\hat{i}_{+1}\hat{s}_{-1}-\hat{i}_{-1}\hat{s}_{+1}}$.
Matrix elements of the spherical components of the electron spin operators are
given in the Appendix. The action of the nuclear spin operators is

$\displaystyle\hat{i}_{z}|im\rangle$ $\displaystyle=m|im\rangle,$ (22)
$\displaystyle\hat{i}_{+1}|im\rangle$
$\displaystyle=-\sqrt{\frac{1}{2}(i-m)(i+m+1)}\ |i,m+1\rangle,$
$\displaystyle\hat{i}_{-1}|im\rangle$
$\displaystyle=\sqrt{\frac{1}{2}(i+m)(i-m+1)}\ |i,m-1\rangle.$ (23)

Inclusion of hyperfine substantially increases the dimension of the basis set
and we performed calculations for modest $N_{\mathrm{max}}=10$, shown in Fig.
11. As can be seen, the inclusion of hyperfine structure modifies the spectrum
of resonances, but it does not substantially increase the number of resonances
nor does it lead to clearly identifiable multiplet structure on existing
resonances.

Figure 12:  Collisional loss rates at higher magnetic fields. Plotted are
representative magnetic field scans up to higher magnetic fields than are
probed experimentally. In the upper spin-stretched state, a strong feature
appears which is caused by a resonance between the Zeeman relaxation energy
and a rotational excitation. In the lower spin-stretched state, more
resonances are typically visible, but the strong feature near 5000 G is
missing as Zeeman relaxation is not possible in this spin state.

###  III.8 Higher magnetic field strengths

Finally, we consider a wider $B$-field scan of the collisional loss rate that
is not accessible in the experiment. In Fig. 12 we see that a strong resonance
feature occurs around 5000 G for the upper spin-stretched state. This feature
occurs consistently at the same magnetic field strength and is insensitive to
uncertainty of the interaction potential. This is not a Feshbach resonance,
but rather results from resonant energy transfer via spin-spin coupling where
the energy released by Zeeman relaxation matches the rotational energy
associated with excitation from $N=0$ to $N=2$. This occurs only in the upper
spin-stretched state, as Zeeman relaxation is not possible in the lower state.

A similar feature may be expected around 1650 G and 3300 G, where the Zeeman
energy released by a double and a single spin flip, respectively, is resonant
with the $N=0$ to $N=1$ rotational transition. This, however, is not observed
as these rotational states are not coupled by spin-spin nor by spin-rotation
coupling.

###  III.9 Discussion

A summary of our findings is represented in Fig. 13 which can also be found in
the accompanying paper com . This figure compares the spectra of Feshbach
resonances for the upper and lower spin-stretched states. Markers show four
times the loss rate obtained with the spin-spin and spin-rotation couplings
halved. Agreement with the solid lines indicates these spin-dependent
couplings act perturbatively, and the loss rates scale as the coupling
constants squared. Note that this does not apply to the sharp resonances,
which are narrower for smaller coupling. The dashed-dotted line indicates the
loss rate obtained with the interaction anisotropy turned off. Turning off the
interaction anisotropy effectively turns off the dominant loss mechanism,
which requires the combination of spin-spin or spin-rotation coupling and the
anisotropic atom-molecule interaction. The resulting, much smaller loss rate
is entirely due to the magnetic dipole-dipole interaction, and vanishes within
our model if the dipole-dipole interaction is also turned off. The role of the
anisotropic atom-molecule interaction, however, is not perturbative. This can
be seen from the disagreement between the solid line, which resulted from the
full calculation, and the dotted line, which is obtained from a calculation
where we halved the strength of anisotropic interactions and multiplying the
resulting loss rate by a factor four. We also see that changing the strength
of the anisotropic short-range interaction affects the resonance positions.
This means that the pattern of resonances cannot be explained in terms of the
long-range _isotropic_ van der Waals interaction alone, as is argued in Ref.
Frye and Hutson (2022), despite the fact that vast majority of resonances is
supported by the long-range $R^{-6}$ interaction, as argued in Ref. Frye and
Hutson (2022) and confirmed in our density of states calculations.

Figure 13:  Perturbative analysis of calculated collision rates. Resonance
spectra are displayed for the upper and lower spin-stretched states in orange
and blue solid lines, respectively, calculated for $\lambda=-0.02$ and
$N_{\mathrm{max}}=30$, also shown as Fig. 4(a) in the accompanying paper com .
Markers show four times the loss rate obtained with the spin-spin and spin-
rotation couplings halved, and agreement with the solid lines indicates these
spin-dependent couplings act perturbatively. For the upper stretched state,
the dashed-dotted line indicates the loss rate obtained with the interaction
anisotropy turned off, which is entirely due to the magnetic dipole-dipole
interaction. The much smaller loss rate in this case identifies a combination
of the anisotropic interaction and the spin-spin and spin-rotation coupling as
the dominant loss mechanism. The dotted line is obtained by halving the
strength of anisotropic interactions and multiplying the resulting loss rate
by a factor four. This completely changes the shape and resonance positions.
The disagreement between the dotted and the full line indicates that the
anisotropic interactions do not act perturbatively.

In Figure 13 we also witness the qualitative difference in the number of
resonances observed in the two maximally stretched spin states. By performing
calculations for various $\lambda$ scalings of the potential, we _typically_
observe around 5 and 10 resonances for the upper and lower spin-stretched
state, respectively. This is in qualitative agreement with experimental
observation of 17 and 8 resonances in the upper and lower stretched state,
respectively. The difference between theory and experiment is partially
explained by the neglect of the vibrational degree of freedom and hyperfine,
as discussed above. This difference between the two spin states in observed
resonance density is not attributed to a difference in density of states, but
rather to the decay of resonances to Zeeman-relaxed channels from the upper
spin-stretched state, as discussed in Sec. III.6.

Figure 13 also suggests substantial differences in the background loss rate
between the upper and lower spin-stretched states, in contrast to the
experimental observations where very similar loss rates were found. Indeed,
this demonstrates that for a particular realization of the calculation, for
specific $\lambda$ scalings of the potential, such differences can occur, and
there is no guarantee that the two state exhibit the same background loss rate
despite the losses being determined by the same mechanisms. However, this is
not necessarily the _expected_ behavior. In the $\lambda$-scans of the ratio
of elastic-to-inelastic collisions, shown in Fig. 5, we see that the two
stretched spin-states exhibit qualitatively the same behavior, and that for
most values of $\lambda$ the loss rates are similar between the two states,
but for specific values of $\lambda$ there can be large differences. Such
large differences can occur where a resonance occurs in one of the two spin
states, but they can also reflect differences in the background scattering
rates. In Figs. 7, 9, 10 we see several examples where the background loss
rate between the two spin states can be either similar, different by a small
factor, or different about an order of magnitude. Where the background
scattering rates are different – for a specific $\lambda$ scaling and basis-
set truncation, $N_{\mathrm{max}}$ – these differences are not systematic; it
can either be the upper or lower spin-stretched state that experiences the
higher background scattering rate. Again, as shown most clearly in Fig. 5, the
_expected_ background behavior is similar for the two spin states.

##  IV Conclusions

We have performed coupled-channels scattering calculations of Feshbach
resonances in spin-polarized NaLi $(a^{3}\Sigma^{+})$ \+ Na collisions based
on _ab initio_ interaction potentials calculated in this work. Quantitatively
predicting the background scattering length or the positions of resonances is
beyond the reach of state-of-the-art theory. However, the calculations do
explain experimental observations qualitatively. The background loss is fast
in non-stretched spin states, whereas in stretched states the ratio of
elastic-to-inelastic collisions can be around 100, in agreement with previous
observations of sympathetic cooling. When comparing the upper and lower
stretched states we find the _expected_ background loss rate to be similar,
also in agreement with experimental observations.

The calculations furthermore capture a series of Feshbach resonances. We show
that these resonance states are supported by relatively short atom-molecule
separations up to $40~{}a_{0}$, We show that the dominant coupling mechanism
is a combination of the anisotropic atom-molecule interaction and the spin-
spin coupling and to a lesser extent spin-rotation coupling. The resonance
states are supported by atom-molecule separations up to $40~{}a_{0}$, where
the interaction can be described by its asymptotic $R^{-6}$ form, as was
previously argued for alkali-metal atom-molecule collisions Wang  _et al._
(2021); Frye and Hutson (2022). However, due to the critical role of the
_anisotropic_ atom-molecule interaction in the coupling mechanism, the
resonance positions depend sensitively on the anisotropic short-range
interactions, and the pattern of resonances cannot be described in terms of
the long-range $R^{-6}$ interaction alone.

In the lower spin-stretched state we observe approximately 10 resonances up to
1500 G. In the upper spin-stretched state only around 5 resonances are visible
due to fast decay to lower-lying Zeeman states in $N=0$. This qualitative
difference between the upper and lower spin-stretched states has also been
observed experimentally, where the two states support 8 and 17 resonances,
respectively. Molecular vibrations and hyperfine interactions, which were
excluded in most of the scattering calculations that we performed, are
expected to further increase the number of observable resonances. Calculations
of the density of states suggests molecular vibrations increase the number of
resonances by 50 %, and a scattering calculation including hyperfine
interaction in a small basis suggests that this too can increase the number of
observed resonances somewhat. Hence, the calculations are in semi-quantitative
agreement with the experimental observations.

Our combined experimental and theoretical study has shown that Feshbach
resonances and collisional complexes can be well understood on the basis of
state-of-the-art first principles calculations. Due to the light elements,
highly accurate electronic structure calculations can be performed on this
triatomic collision complex, with an uncertainty smaller than 3% of the
pairwise nonadditive three-body part of the potential. To make fully
quantitative predictions, the accuracy of the electronic structure
calculations should be improved by more than one order of magnitude, and the
coupled-channels calculations should be converged, including vibration and
hyperfine structure. While it is not feasible with present-day methods and
computational power, predicting the exact positions of measured Feshbach
resonances in NaLi+Na collisions constitutes a perfect testbed and playground
for near-future developments of both electronic structure and scattering
theory. In the mean time, this work emphasizes we can still draw statistical
conclusions regarding the density and typical width of resonances, and the
typical background loss rate, that are in nearly quantitative agreement with
experimental observations. Our work showcases how scattering calculations can
be used as a “ _numerical experiment_ ” in which the interaction can be scaled
at will, specific couplings can be turned off, or exit channels removed, as a
versatile tool to identify the dominant coupling mechanisms, the nature of the
resonance states, and their relevant decay pathways.

###### Acknowledgements.

M. G. and M. T. were supported by the National Science Centre Poland (Sonata
Bis Grant no. 2020/38/E/ST2/00564) and the PL-Grid Infrastructure
(computational grant no. PLG/2020/014342). The MIT work was supported by the
NSF through the Center for Ultracold Atoms and Grant No. 1506369 and from the
Air Force Office of Scientific Research (MURI, Grant No. FA9550-21-1-0069). J.
P. and H. S. acknowledge additional support from the Samsung Scholarship.

##  Appendix A Matrix elements

In this Appendix we give all required matrix elements in the coupled basis
functions of Eq. (9). First we consider the matrix elements of the electronic
spin operators required for the Zeeman interaction
$\hat{H}_{\mathrm{Zeeman}}=\mu_{B}g_{e}B(\hat{s}_{1,z}+\hat{s}_{2,z}+\hat{s}_{3,z})$.
These are

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|\hat{s}_{1,q}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}^{\prime}\mathcal{M}\rangle=$
$\displaystyle\delta_{N,N^{\prime}}\delta_{L,L^{\prime}}\delta_{J,J^{\prime}}(-1)^{2\mathcal{J}-\mathcal{M}+J+2S^{\prime}+s+s^{\prime}+s_{1}+s_{2}+s_{3}+1}\left[\mathcal{J},\mathcal{J}^{\prime},S,S^{\prime},s,s^{\prime}\right]^{1/2}$
$\displaystyle\times\left(\begin{array}[]{ccc}\mathcal{J}&1&\mathcal{J}^{\prime}\\\
\mathcal{M}&q&\mathcal{M}^{\prime}\end{array}\right)\left\\{\begin{array}[]{ccc}\mathcal{J}&1&\mathcal{J}\\\
S^{\prime}&J&S\end{array}\right\\}\left\\{\begin{array}[]{ccc}S&1&S^{\prime}\\\
s^{\prime}&s_{3}&s\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&1&s^{\prime}\\\
s_{1}&s_{2}&s_{1}\end{array}\right\\}\sqrt{s_{1}(s_{1}+1)(2s_{1}+1)},$ (32)

and similarly

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|\hat{s}_{2,q}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}^{\prime}\mathcal{M}\rangle=$
$\displaystyle\delta_{N,N^{\prime}}\delta_{L,L^{\prime}}\delta_{J,J^{\prime}}(-1)^{2\mathcal{J}-\mathcal{M}+J+2S^{\prime}+2s^{\prime}+s_{1}+s_{2}+s_{3}+1}\left[\mathcal{J},\mathcal{J}^{\prime},S,S^{\prime},s,s^{\prime}\right]^{1/2}$
$\displaystyle\times\left(\begin{array}[]{ccc}\mathcal{J}&1&\mathcal{J}^{\prime}\\\
\mathcal{M}&q&\mathcal{M}^{\prime}\end{array}\right)\left\\{\begin{array}[]{ccc}\mathcal{J}&1&\mathcal{J}\\\
S^{\prime}&J&S\end{array}\right\\}\left\\{\begin{array}[]{ccc}S&1&S^{\prime}\\\
s^{\prime}&s_{3}&s\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&1&s^{\prime}\\\
s_{2}&s_{1}&s_{2}\end{array}\right\\}\sqrt{s_{2}(s_{2}+1)(2s_{2}+1)},$ (41)

and

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|\hat{s}_{3,q}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}^{\prime}\mathcal{M}\rangle=$
$\displaystyle\delta_{N,N^{\prime}}\delta_{L,L^{\prime}}\delta_{s,s^{\prime}}\delta_{J,J^{\prime}}(-1)^{2\mathcal{J}-\mathcal{M}+J+2S^{\prime}+s+s_{3}}\left[\mathcal{J},\mathcal{J}^{\prime},S,S^{\prime}\right]^{1/2}$
$\displaystyle\times\left(\begin{array}[]{ccc}\mathcal{J}&1&\mathcal{J}^{\prime}\\\
\mathcal{M}&q&\mathcal{M}^{\prime}\end{array}\right)\left\\{\begin{array}[]{ccc}\mathcal{J}&1&\mathcal{J}\\\
S^{\prime}&J&S\end{array}\right\\}\left\\{\begin{array}[]{ccc}S&1&S^{\prime}\\\
s^{\prime}&s_{3}&s\end{array}\right\\}\sqrt{s_{3}(s_{3}+1)(2s_{3}+1)},$ (48)

where the $q=-1,0,1$ spherical components of the spin operators are given by
$\hat{s}_{0}=\hat{s}_{z}$ and $\hat{s}_{\pm 1}=\mp(\hat{s}_{x}\pm
i\hat{s}_{y})/\sqrt{2}$.

The remaining operators we consider are all scalar, and diagonal in
$\mathcal{J}$ and $\mathcal{M}$. For the spin-rotation coupling
$\gamma_{s}\hat{N}\cdot\hat{s}$ with $\gamma_{s}=0.005~{}$cm-1 we need the
matrix elements

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|\hat{N}\cdot\hat{s}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}\mathcal{M}\rangle=$
$\displaystyle\delta_{N,N^{\prime}}\delta_{L,L^{\prime}}\delta_{S,S^{\prime}}(-1)^{2J^{\prime}+S+\mathcal{J}+N+L+s+s_{3}+S^{\prime}}\left[J,J^{\prime},S,S^{\prime}\right]^{1/2}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}J&1&J^{\prime}\\\
N^{\prime}&L&N\end{array}\right\\}\left\\{\begin{array}[]{ccc}S&1&S^{\prime}\\\
s^{\prime}&s_{3}&s\end{array}\right\\}\sqrt{N(N+1)(2N+1)s(s+1)(2s+1)}.$ (53)

The spin-spin coupling is given by
$\lambda_{s}\sqrt{30}/3\left[\left[\hat{s}\otimes\hat{s}\right]^{(2)}\otimes
C^{(2)}(\hat{r}_{\mathrm{NaLi}})\right]^{(0)}_{0}$ with
$\lambda_{s}=-0.0189~{}$cm-1. Here,

$\displaystyle\left[A^{(k_{1})}\otimes
B^{(k_{2})}\right]^{(k)}_{q}=\sum_{q_{1},q_{2}}A^{(k_{1})}_{q_{1}}B^{(k_{2})}_{q_{2}}\langle
k_{1}q_{1}k_{2}q_{2}|kq\rangle$ (54)

is the rank-$k$ spherical tensor product and $C^{(2)}(\hat{R})$ is a tensor of
Racah-normalized spherical harmonics depending on the Euler angles of the
molecular axis, $\bm{r}_{\mathrm{NaLi}}$. For the matrix elements we have

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|\left[\left[\hat{s}\otimes\hat{s}\right]^{(2)}\otimes
C^{(2)}(\hat{r}_{\mathrm{NaLi}})\right]^{(0)}_{0}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}\mathcal{M}\rangle=$
$\displaystyle\delta_{L,L^{\prime}}\delta_{s,s^{\prime}}(-1)^{2J^{\prime}+S+\mathcal{J}+N+L+s+s_{3}+S^{\prime}}\left[J,J^{\prime},S,S^{\prime}\right]^{1/2}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}J&2&J^{\prime}\\\
N^{\prime}&L&N\end{array}\right\\}\left\\{\begin{array}[]{ccc}S&2&S^{\prime}\\\
s^{\prime}&s_{3}&s\end{array}\right\\}\left(\begin{array}[]{ccc}N&2&N^{\prime}\\\
0&0&0\end{array}\right)\frac{\sqrt{s(2s-1)}}{\sqrt{5}\left(\begin{array}[]{ccc}s&2&s\\\
-s&2&s-2\end{array}\right)}.$ (63)

In order to evaluate the magnetic dipole-dipole interaction
$\hat{V}_{\mathrm{magn.dip}}=-\sqrt{30}(\mu_{B}g_{e}\alpha)^{2}R^{-3}\left[\left[\hat{s}\otimes\hat{s}_{3}\right]^{(2)}\otimes
C^{(2)}(\hat{R})\right]^{(0)}_{0}$, we use the following

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|\left[\left[\hat{s}\otimes\hat{s}_{3}\right]^{(2)}\otimes
C^{(2)}(\hat{R})\right]^{(0)}_{0}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}\mathcal{M}\rangle=$
$\displaystyle\delta_{N,N^{\prime}}\delta_{s,s^{\prime}}(-1)^{J^{\prime}+S+\mathcal{J}+N+L^{\prime}+J+L}\left[L,L^{\prime},J,J^{\prime},S,S^{\prime}\right]^{1/2}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}J&2&J^{\prime}\\\
L^{\prime}&N&L\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&s^{\prime}&1\\\
s_{3}&s_{3}&1\\\ S&S^{\prime}&2\\\
\end{array}\right\\}\left(\begin{array}[]{ccc}L&2&L^{\prime}\\\
0&0&0\end{array}\right)\sqrt{s(s+1)(2s+1)s_{3}(s_{3}+1)(2s_{3}+1)}.$ (71)

Finally, for the electronic interaction we have the following expressions

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|P_{\ell}(\cos\theta)|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}\mathcal{M}\rangle=$
$\displaystyle\delta_{J,J^{\prime}}\delta_{s,s^{\prime}}\delta_{S,S^{\prime}}(-1)^{J+S+\mathcal{J}+N+L}\left[\mathcal{J},J,J^{\prime},\ell,N,N^{\prime},L,L^{\prime}\right]^{1/2}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}N&N^{\prime}&\ell\\\
L&L^{\prime}&\ell\\\ J&J^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}\mathcal{J}&0&\mathcal{J}^{\prime}\\\
J^{\prime}&S&J\end{array}\right\\}\left(\begin{array}[]{ccc}N&\ell&N^{\prime}\\\
0&0&0\end{array}\right)\left(\begin{array}[]{ccc}L&\ell&L^{\prime}\\\
0&0&0\end{array}\right),$ (81) $\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|P_{\ell}(\cos\theta)\hat{s}_{1}\cdot\hat{s}_{3}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}\mathcal{M}\rangle=$
$\displaystyle\delta_{J,J^{\prime}}\delta_{S,S^{\prime}}(-1)^{N+L+s_{1}+s_{2}+s^{\prime}}\left[\mathcal{J},J,J^{\prime},\ell,N,N^{\prime},L,L^{\prime},S,S^{\prime},s,s^{\prime},1\right]^{1/2}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}N&N^{\prime}&\ell\\\
L&L^{\prime}&\ell\\\ J&J^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}J&J^{\prime}&0\\\
S&S^{\prime}&0\\\ \mathcal{J}&\mathcal{J}^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&s^{\prime}&1\\\
s_{3}&s_{3}&1\\\ S&S^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&1&s^{\prime}\\\
s_{1}&s_{2}&s_{1}\end{array}\right\\}\left(\begin{array}[]{ccc}N&\ell&N^{\prime}\\\
0&0&0\end{array}\right)\left(\begin{array}[]{ccc}L&\ell&L^{\prime}\\\
0&0&0\end{array}\right)$ (97)
$\displaystyle\times\sqrt{s_{1}(s_{1}+1)(2s_{1}+1)s_{3}(s_{3}+1)(2s_{3}+1)}$
(98)

and

$\displaystyle\langle(NL)J(s~{}s_{3})S;$
$\displaystyle\mathcal{J}\mathcal{M}|P_{\ell}(\cos\theta)\hat{s}_{2}\cdot\hat{s}_{3}|(N^{\prime}L^{\prime})J^{\prime}(s^{\prime}~{}s_{3})S^{\prime};\mathcal{J}\mathcal{M}\rangle=$
$\displaystyle\delta_{J,J^{\prime}}\delta_{S,S^{\prime}}(-1)^{N+L+s_{1}+s_{2}+s}\left[\mathcal{J},J,J^{\prime},\ell,N,N^{\prime},L,L^{\prime},S,S^{\prime},s,s^{\prime},1\right]^{1/2}$
$\displaystyle\times\left\\{\begin{array}[]{ccc}N&N^{\prime}&\ell\\\
L&L^{\prime}&\ell\\\ J&J^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}J&J^{\prime}&0\\\
S&S^{\prime}&0\\\ \mathcal{J}&\mathcal{J}^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&s^{\prime}&1\\\
s_{3}&s_{3}&1\\\ S&S^{\prime}&0\\\
\end{array}\right\\}\left\\{\begin{array}[]{ccc}s&1&s^{\prime}\\\
s_{2}&s_{1}&s_{2}\end{array}\right\\}\left(\begin{array}[]{ccc}N&\ell&N^{\prime}\\\
0&0&0\end{array}\right)\left(\begin{array}[]{ccc}L&\ell&L^{\prime}\\\
0&0&0\end{array}\right)$ (114)
$\displaystyle\times\sqrt{s_{2}(s_{2}+1)(2s_{2}+1)s_{3}(s_{3}+1)(2s_{3}+1)}.$
(115)

The full expansion of the interaction is given in Eq. (3), and the expansion
coefficients are determined as explained in the main text.

## References

  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Bloch  _et al._ (2012) I. Bloch, J. Dalibard, and S. Nascimbene, Nature Phys. 8, 267 (2012).
  * Bakr  _et al._ (2010) W. S. Bakr, A. Peng, M. E. Tai, R. Ma, J. Simon, J. I. Gillen, S. Foelling, L. Pollet, and M. Greiner, Science 329, 547 (2010).
  * Struck  _et al._ (2011) J. Struck, C. Ölschläger, R. Le Targat, P. Soltan-Panahi, A. Eckardt, M. Lewenstein, P. Windpassinger, and K. Sengstock, Science 333, 996 (2011).
  * Simon  _et al._ (2011) J. Simon, W. S. Bakr, R. Ma, M. E. Tai, P. M. Preiss, and M. Greiner, Nature 472, 307 (2011).
  * Chiu  _et al._ (2018) C. S. Chiu, G. Ji, A. Mazurenko, D. Greif, and M. Greiner, Phys. Rev. Lett. 120, 243201 (2018).
  * Köhler  _et al._ (2006) T. Köhler, K. Góral, and P. S. Julienne, Rev. Mod. Phys. 78, 1311 (2006).
  * Bergmann  _et al._ (1998) K. Bergmann, H. Theuer, and B. W. Shore, Rev. Mod. Phys. 70, 1003 (1998).
  * Ni  _et al._ (2008) K. K. Ni, S. Ospelkaus, M. H. G. de Miranda, A. Pe’er, B. Neyenhuis, J. J. Zirbel, S. Kotochigova, P. S. Julienne, D. S. Jin, and J. Ye, Science 322, 231 (2008).
  * Takekoshi  _et al._ (2014) T. Takekoshi, L. Reichsöllner, A. Schindewolf, J. M. Hutson, C. R. Le Sueur, O. Dulieu, F. Ferlaino, R. Grimm, and H. C. Nägerl, Phys. Rev. Lett. 113, 205301 (2014).
  * Molony  _et al._ (2014) P. K. Molony, P. D. Gregory, Z. Ji, B. Lu, M. P. Köppinger, C. R. Le Sueur, C. L. Blackley, J. M. Hutson, and S. L. Cornish, Phys. Rev. Lett. 113, 255301 (2014).
  * Guo  _et al._ (2016) M. Guo, B. Zhu, B. Lu, X. Ye, F. Wang, R. Vexiau, N. Bouloufa-Maafa, G. Quéméner, O. Dulieu, and D. Wang, Phys. Rev. Lett. 116, 205303 (2016).
  * Son  _et al._ (2020) H. Son, J. J. Park, W. Ketterle, and A. O. Jamison, Nature 580, 197 (2020).
  * Yang  _et al._ (2022) H. Yang, X.-Y. Wang, Z. Su, J. Cao, D.-C. Zhang, J. Rui, B. Zhao, C.-L. Bai, and J.-W. Pan, Nature 602, 229 (2022).
  * Gronowski  _et al._ (2020) M. Gronowski, A. M. Koza, and M. Tomza, Phys. Rev. A 102, 020801 (2020).
  * Mayle  _et al._ (2012) M. Mayle, B. P. Ruzic, and J. L. Bohn, Phys. Rev. A 85, 062712 (2012).
  * Mayle  _et al._ (2013) M. Mayle, G. Quéméner, B. P. Ruzic, and J. L. Bohn, Phys. Rev. A 87, 012709 (2013).
  * Christianen  _et al._ (2019a) A. Christianen, T. Karman, and G. C. Groenenboom, Phys. Rev. A 100, 032708 (2019a).
  * Christianen  _et al._ (2019b) A. Christianen, M. W. Zwierlein, G. C. Groenenboom, and T. Karman, Phys. Rev. Lett. 123, 123402 (2019b).
  * Croft  _et al._ (2021) J. F. E. Croft, J. L. Bohn, and G. Quéméner, (2021), arXiv:2111.09956 .
  * Christianen  _et al._ (2021) A. Christianen, G. C. Groenenboom, and T. Karman, Phys. Rev. A 104, 043327 (2021), 2108.02724 .
  * Man  _et al._ (2022) M. P. Man, T. Karman, and G. C. Groenenboom, New J. Phys. 24, 055001 (2022), 2108.02511 .
  * Croft  _et al._ (2020) J. F. E. Croft, J. L. Bohn, and G. Quéméner, Phys. Rev. A 102, 033306 (2020).
  * Idziaszek and Julienne (2010) Z. Idziaszek and P. S. Julienne, Phys. Rev. Lett. 104, 113202 (2010).
  * Danzl  _et al._ (2010) J. G. Danzl, M. J. Mark, E. Haller, M. Gustavsson, R. Hart, J. Aldegunde, J. M. Hutson, and H. C. Nägerl, Nature Phys. 6, 265 (2010).
  * Park  _et al._ (2015) J. W. Park, S. A. Will, and M. W. Zwierlein, Phys. Rev. Lett. 114, 205302 (2015).
  * Rvachov  _et al._ (2017) T. M. Rvachov, H. Son, A. T. Sommer, S. Ebadi, J. J. Park, M. W. Zwierlein, W. Ketterle, and A. O. Jamison, Phys. Rev. Lett. 119, 143001 (2017).
  * Seeßelberg  _et al._ (2018) F. Seeßelberg, N. Buchheim, Z.-K. Lu, T. Schneider, X.-Y. Luo, E. Tiemann, I. Bloch, and C. Gohle, Phys. Rev. A 97, 013405 (2018).
  * Yang  _et al._ (2019) H. Yang, D.-C. Zhang, L. Liu, Y.-X. Liu, J. Nan, B. Zhao, and J.-W. Pan, Science 363, 261 (2019).
  * Gregory  _et al._ (2021) P. D. Gregory, J. A. Blackmore, L. M. Fernley, S. L. Bromley, J. M. Hutson, S. L. Cornish, _et al._ , New J. Phys. 23, 125004 (2021).
  * Wang  _et al._ (2021) X.-Y. Wang, M. D. Frye, Z. Su, J. Cao, L. Liu, D.-C. Zhang, H. Yang, J. M. Hutson, B. Zhao, C.-L. Bai, and J.-W. Pan, New J. Phys. 23, 115010 (2021).
  * Nichols  _et al._ (2022) M. A. Nichols, Y.-X. Liu, L. Zhu, M.-G. Hu, Y. Liu, and K.-K. Ni, Phys. Rev. X 12, 011049 (2022).
  * Voges  _et al._ (2022) K. K. Voges, P. Gersema, T. Hartmann, S. Ospelkaus, and A. Zenesini, Phys. Rev. Res. 4, 023184 (2022).
  * Son  _et al._ (2022) H. Son, J. J. Park, Y.-K. Lu, A. O. Jamison, T. Karman, and W. Ketterle, Science 375, 1006 (2022).
  * Park  _et al._ (2022) J. J. Park, Y.-K. Lu, A. O. Jamison, T. Tscherbul, and W. Ketterle, “A feshbach resonance in collisions between ultracold ground state molecules,”  (2022).
  * (36) “Accompanying paper,” .
  * Knoop  _et al._ (2011) S. Knoop, T. Schuster, R. Scelle, A. Trautmann, J. Appmeier, M. K. Oberthaler, E. Tiesinga, and E. Tiemann, Phys. Rev. A 83, 042704 (2011).
  * Steinke  _et al._ (2012) M. Steinke, H. Knöckel, and E. Tiemann, Phys. Rev. A 85, 042720 (2012).
  * Rvachov  _et al._ (2018) T. M. Rvachov, H. Son, J. J. Park, S. Ebadi, M. W. Zwierlein, W. Ketterle, and A. O. Jamison, Phys. Chem. Chem. Phys. 20, 4739 (2018).
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
  * Knowles  _et al._ (1993) P. J. Knowles, C. Hampel, and H.-J. Werner, J. Chem. Phys. 99, 5219 (1993).
  * Kállay and Surján (2001) M. Kállay and P. R. Surján, J. Chem. Phys. 115, 2945 (2001).
  * Prascher  _et al._ (2011) B. P. Prascher, D. E. Woon, K. A. Peterson, T. H. Dunning, and A. K. Wilson, Theor. Chem. Acc. 128, 69–82 (2011).
  * Helgaker  _et al._ (1997) T. Helgaker, W. Klopper, H. Koch, and J. Noga, J. Chem. Phys. 106, 9639 (1997).
  * Peng and Reiher (2012) D. Peng and M. Reiher, Theor. Chem. Acc. 131, 1081 (2012).
  * Werner  _et al._ (2012) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, and M. Schütz, WIREs Comput. Mol. Sci. 2, 242 (2012).
  * (47) H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, P. Celani, W. Györffy, D. Kats, T. Korona, R. Lindh, A. Mitrushenkov, G. Rauhut, K. R. Shamasundar, T. B. Adler, R. D. Amos, S. J. Bennie, A. Bernhardsson, A. Berning, D. L. Cooper, M. J. O. Deegan, A. J. Dobbyn, F. Eckert, E. Goll, C. Hampel, A. Hesselmann, G. Hetzer, T. Hrenar, G. Jansen, C. Köppl, S. J. R. Lee, Y. Liu, A. W. Lloyd, Q. Ma, R. A. Mata, A. J. May, S. J. McNicholas, W. Meyer, T. F. Miller III, M. E. Mura, A. Nicklass, D. P. O’Neill, P. Palmieri, D. Peng, K. Pflüger, R. Pitzer, M. Reiher, T. Shiozaki, H. Stoll, A. J. Stone, R. Tarroni, T. Thorsteinsson, M. Wang, and M. Welborn, “Molpro, version 2019, a package of ab initio programs,” See https://www.molpro.net.
  * Kállay  _et al._ (2020) M. Kállay, P. R. Nagy, D. Mester, Z. Rolik, G. Samu, J. Csontos, J. Csóka, P. B. Szabó, L. Gyevi-Nagy, B. Hégely, I. Ladjánszki, L. Szegedy, B. Ladóczki, K. Petrov, M. Farkas, P. D. Mezei, and A. Ganyecz, J. Chem. Phys. 152, 074107 (2020).
  * (49) M. Kállay, P. R. Nagy, D. Mester, Z. Rolik, G. Samu, J. Csontos, J. Csóka, P. B. Szabó, L. Gyevi-Nagy, B. Hégely, I. Ladjánszki, L. Szegedy, B. Ladóczki, K. Petrov, M. Farkas, P. D. Mezei, and A. Ganyecz, “Mrcc, a quantum chemical program,” See www.mrcc.hu.
  * Jaszuński  _et al._ (2017) M. Jaszuński, A. Rizzo, and K. Ruud, “Molecular electric, magnetic, and optical properties,” in  _Handbook of Computational Chemistry_, edited by J. Leszczynski, A. Kaczmarek-Kedziera, T. Puzyn, M. G. Papadopoulos, H. Reis, and M. K. Shukla (Springer International Publishing, Cham, 2017) pp. 497–592.
  * Szalay and Bartlett (1993) P. G. Szalay and R. J. Bartlett, Chem. Phys. Lett. 214, 481 (1993).
  * Ganyushin and Neese (2006) D. Ganyushin and F. Neese, J. Chem. Phys. 125, 024103 (2006).
  * Neese (2012) F. Neese, WIREs Comput. Mol. Sci. 2, 73 (2012).
  * Neese (2018) F. Neese, WIREs Comput. Mol. Sci. 8, e1327 (2018).
  * Dyall (2000) K. G. Dyall, Int. J. Quantum Chem. 78, 412 (2000).
  * Sandhoefer and Neese (2012) B. Sandhoefer and F. Neese, J. Chem. Phys. 137, 094102 (2012).
  * Ogura (1920) K. Ogura, Tohoku Math. J. 18, 61 (1920).
  * Ho and Rabitz (1996) T.-S. Ho and H. Rabitz, J. Chem. Phys. 104, 2584 (1996).
  * Rolik and Kállay (2011) Z. Rolik and M. Kállay, J. Chem. Phys. 134, 124111 (2011).
  * Bauer and Toennies (2019) J. B. Bauer and J. P. Toennies, J. Chem. Phys. 150, 144310 (2019).
  * Tomza  _et al._ (2013) M. Tomza, K. W. Madison, R. Moszynski, and R. V. Krems, Phys. Rev. A 88, 050701 (2013).
  * Janssen  _et al._ (2013) L. M. C. Janssen, A. van der Avoird, and G. C. Groenenboom, Phys. Rev. Lett. 110, 063201 (2013), arXiv:1301.1931 .
  * Karman (0) T. Karman, J. Phys. Chem. A 0, 0 (0), arXiv:2212.03065 .
  * Morita  _et al._ (2018) M. Morita, M. B. Kosicki, P. S. Żuchowski, and T. V. Tscherbul, Phys. Rev. A 98, 042702 (2018).
  * Morita  _et al._ (2019) M. Morita, R. V. Krems, and T. V. Tscherbul, Phys. Rev. Lett. 123, 013401 (2019).
  * Morita  _et al._ (2020) M. Morita, J. Kłos, and T. V. Tscherbul, Phys. Rev. Research 2, 043294 (2020).
  * Gribakin and Flambaum (1993) G. F. Gribakin and V. V. Flambaum, Phys. Rev. A 48, 546 (1993).
  * Frye and Hutson (2022) M. D. Frye and J. M. Hutson, (2022), arXiv:2212.08030 .
  * Colbert and Miller (1992) D. T. Colbert and W. H. Miller, J. Chem. Phys. 96, 1982 (1992).