# Photoassociation of ultracold long-range polyatomic molecules

Marko Gacesa  Bay Area Environmental Research Institute, Moffett Field, CA
94035-0001, USA  Space Science Division, NASA Ames Research Center, MS 245-3,
Moffett Field, CA 94035, USA  Department of Physics, Khalifa University, Abu
Dhabi 127788, United Arab Emirates<EMAIL_ADDRESS>Jason N. Byrd
ENSCO, Inc., 4849 North Wickham Road, Melbourne, FL 32940, USA  Jonathan
Smucker  Department of Physics, University of Connecticut, Storrs, CT
06269-3046, USA  John A. Montgomery, Jr  Department of Physics, University of
Connecticut, Storrs, CT 06269-3046, USA  Robin Côté  Department of Physics,
University of Connecticut, Storrs, CT 06269-3046, USA  Physics Department,
University of Massachusetts Boston, Boston, Massachusetts 02125, USA
<EMAIL_ADDRESS>
###### Abstract

We explore the feasibility of optically forming long-range tetratomic and
larger polyatomic molecules in their ground electronic state from ultracold
pairs of polar molecules aligned by external fields. Depending on the relative
orientation of the interacting diatomic molecules, we find that a tetratomic
can be formed either as a weakly bound complex in a very extended halo state
or as a pure long-range molecule composed of collinear or nearly-collinear
diatomic molecules. The latter is a novel type of tetratomic molecule
comprised of two diatomic molecules bound at long intermolecular range and
predicted to be stable in cold and ultracold regimes. Our numerical studies
were conducted for ultracold KRb and RbCs, resulting in production of (KRb)2
and (RbCs)2 complexes, respectively. Based on universal properties of long-
range interactions between polar molecules, we identify triatomic and
tetratomic linear polar molecules with favorable ratio of dipole and
quadrupole moments for which the apporach could be generalized to form larger
polyatomic molecules.

###### pacs:

##  I Introduction

Ultracold polar molecules have been proposed as an ideal model system to
explore novel physical phenomena at the intersection of molecular physics with
few- and many-body quantum physics. For example, they could be used to
engineer and study lattice spin models in strongly interacting many-body
hamiltoniansYan  _et al._ (2013); Gorshkov  _et al._ (2013); Eisert  _et al._
(2013); Richerme  _et al._ (2014); Doçaj  _et al._ (2016); Kruckenhauser  _et
al._ (2020), supersolidityCapogrosso-Sansone (2011), unconventional superfluid
phases and quantum magnetismHazzard  _et al._ (2013); Wall  _et al._ (2015),
to name a few. Moreover, the tunability of interactions combined with the
extraordinary degree of control available in current experimentsChin  _et al._
(2010) make such systems suitable for engineering quantum simulatorsMicheli
_et al._ (2006); Bloch  _et al._ (2008), as well as quantum entanglement and
information processing DeMille (2002); Krems  _et al._ (2009); Kuznetsova  _et
al._ (2011); Côté (2014); Micheli  _et al._ (2006); Yelin  _et al._ (2006);
Kuznetsova  _et al._ (2010). Unlike atomic gases, where inter-particle
interactions are isotropic and short-range, gases of ultracold polar molecules
exhibit much richer dynamics and macroscopic properties due to the long-range
anisotropic electric dipole-dipole interactions between their constituents.
While dipole-dipole interactions can be realized using atomic magnetic
dipolesStuhler  _et al._ (2005) or between excited Rydberg atoms Boisseau  _et
al._ (2002); Singer  _et al._ (2005), they will be much weaker (ground state
atoms) or be short-lived (Rydberg atoms) than those between permanent electric
dipoles of polar molecules, and allow less tunability via external
electromagnetic fieldsJaksch  _et al._ (2000); Lukin  _et al._ (2001); Cote
(2005); Stanojevic  _et al._ (2006); Hazzard  _et al._ (2014).

Producing ultracold polar molecules was achieved for only a handful of
diatomic species and up to this day remains challenging. Fermionic 40K87RbNi
_et al._ (2008); Moses  _et al._ (2015) were the first heteronuclear molecules
to be produced in deeply bound molecular states, with NaK Park  _et al._
(2015); Seeßelberg  _et al._ (2018), RbCsTakekoshi  _et al._ (2014); Molony
_et al._ (2014); Reichsöllner  _et al._ (2017), NaRbGuo  _et al._ (2016), and
LiNaHeo  _et al._ (2012); Son  _et al._ (2020) being added to the list more
recently. In these experiments, pairs of atoms in an ultracold atomic gas were
coherently associated into loosely bound molecules by a magnetic field sweep
across a Fano-Feshbach resonance, and stimulated Raman adiabatic passage
(STIRAP)Gaubatz  _et al._ (1990) was used to transfer the population into the
molecular ground state. Alternative approaches to produce ultracold molecules
is the photoassociation (PA) of ultracold atomic pairsThorsheim  _et al._
(1987); Jones  _et al._ (2006) or Feshbach Optimized PA (FOPA) Pellegrini  _et
al._ (2008); Pellegrini and Côté (2009); Krzyzewski  _et al._ (2015); Hai  _et
al._ (2020).

Thus far, those techniques were limited to diatomic molecules, even though the
theory suggests that they could be extended to larger molecules Dulieu and
Gabbanini (2009); Balakrishnan (2016). In practice, in addition to a
significantly increased complexity of the experiments, one of the greatest
obstacles to production of larger molecules is the rapid loss of cold
molecules from the trap, caused by chemical reactions (e.g., KRb+KRb
$\rightarrow$ K2+Rb2 is exothermic) or not well-understood “complex nature of
scattering processes” found to occur even for endothermic intermolecular
reactions Reichsöllner  _et al._ (2017). In order to explain the loss due to
bi-molecular exothermic reactions, and how to prevent it, detailed studies of
ultracold reactions between KRb molecules including energetics, were conducted
Byrd  _et al._ (2010); Buchachenko  _et al._ (2012). Byrd et al. Byrd  _et
al._ (2012a) have carried out calculations of structure, energetics, and
reactions of alkali metal tetramer molecules and analyzed conditions for their
controllable binding in a trap with attractive dipoles Byrd  _et al._ (2012b,
c). These studies can be interpreted as initial investigations of electronic
structure of ground-state tetramers composed of two alkali-metal polar
molecules.

In addition, the knowledge of electronic potential energy surfaces (PESs) for
ground and excited electronic states over a wide range of internuclear
separations is critical for experiments’ design and selection of suitable
atomic and molecular systems. Thus, preliminary theoretical studies have been
focused on understanding long-range atom-diatom Lepers  _et al._ (2010);
Pérez-Ríos  _et al._ (2015) and diatom-diatom interactions Byrd  _et al._
(2012c, c, b, a); Vexiau  _et al._ (2015); Lepers  _et al._ (2016); Quéméner
(2017) for experimentally accessible species including Cs-Cs2, KRb-KRb, and
other alkali metal diatomic molecules. Notably, Perez-Rios et al. Pérez-Ríos
_et al._ (2015) computed the rate coefficient for the PA of Cs3 and found that
it would be comparable to the formation rates observed for diatomic cases,
and, more recently, Schnabel et al. Schnabel  _et al._ (2021) theoretically
investigated the prospects for photoassociation (PA) of Rb3.

In this paper, we explore the possibility that ultracold long-range tetratomic
and larger polyatomic molecular complexes can be formed by photoassociation of
pairs of ultracold polar molecules. The properties of the long-range complexes
can be predicted from the ratio of dipole and quadrupole electric moments of
the pairs, while their orientation determines if such long-range states exist.
The concept is robust and extends beyond known diatomic pairs to larger
polyatomic complexes.

This article is organized as follows. In Section II we outline the working
concepts used in the photoassociation of pairs of polar molecules. In Section
III we describe the theoretical methods employed on two benchmark ultracold
molecules, KRb and RbCs, both of which have been experimentally produced. In
Section IV, the numerical results and their implications for the production of
larger polyatomic molecules are presented. We summarize the main findings of
this study and conclude in Section V.

##  II Concept: Photoassociation of diatomic into tetratomic molecules

###  II.1 Background

The photoassociation (PA) of a pair of ultracold atoms was initially proposed
by Thorsheim et al. Thorsheim  _et al._ (1987) in the context of molecular
spectroscopy. The PA is particularly effective at ultralow temperatures, where
two atoms colliding with extremely low relative kinetic energy ($k_{B}T<1$ mK
$\simeq$ 21 MHz) can efficiently absorb a photon from a laser field tuned to a
quasi-resonant free-bound electronic dipolar transition and form a molecule in
an excited electronic state, typically in a weakly bound ro-vibrational level
close to the dissociation energy. The PA has found numerous applications in
low-temperature atomic, molecular, and optical physics, as described in a
number of excellent review articlesLett  _et al._ (1995); Stwalley and Wang
(1999); Jones  _et al._ (2006); Krems  _et al._ (2009); Dulieu and Gabbanini
(2009); Balakrishnan (2016); Côté (2010).

The PA can be used to form ultracold molecules in deeply bound ro-vibrational
levels of the electronic ground state. A demonstrated approach consists of two
steps: the PA of an ultracold atomic pair into a molecule in an excited
electronic state, followed by the relaxation into the ground state via
spontaneous emissionLett  _et al._ (1995); Fioretti  _et al._ (1998);
Vanhaecke  _et al._ (2002, 2004). The main difficulty with the PA formation
can be traced to the Franck-Condon principle, according to which the PA is
almost always most effective at large internuclear separations where the
spontaneous emission has very low probability to populate deeply bound levels.
Early proposals to remedy this obstacle suggested to use re-pumping lasers to
reach excited molecular levels that overlap better with deeply molecular
ground state levels Côté and Dalgarno (1997, 1999). Another approach to
enhancing photoassociative production efficiency of stable ultracold molecules
can be achieved via Feshbach resonances Pellegrini  _et al._ (2008); Hu  _et
al._ (2014); Gacesa and Côté (2014); Hai  _et al._ (2020); Sun  _et al._
(2020), short laser pulses Koch  _et al._ (2006); Ghosal  _et al._ (2009); Hu
_et al._ (2014); Carini  _et al._ (2015); Wang  _et al._ (2017), or both
Gacesa  _et al._ (2013); Hu  _et al._ (2015).

###  II.2 Photoassociation in specific geometries

Figure 1: (Color online) Schematic representation of photoassociative
formation of a ground-state molecule (XY)2 from two aligned ultracold polar
dimers XY whose relative orientation is given by the angle $\theta$. The
excited tetramer in the ro-vib. state $|b\rangle$ (wave function $\psi_{b}$),
detuned by $\Delta$ from its dissociation energy, is created by absorbing a
photon of energy $\hbar\omega$ from the laser field of intensity $I$. The
excited state spontaneously decays into the ground state, ro-vib. level
$|b^{\prime}\rangle$ with binding energy $E_{b^{\prime}}$ (wave function
$\psi_{b^{\prime}}$). Inset: The long-range region of the ground state. Long-
range potential well (shown for $\theta$=0°) capable of supporting several
bound states exists for $\theta<\theta_{c}$, where $\theta_{c}\approx 20$° for
alkali metal pairs.

In the two-step PA formation of tetratomic molecules, and larger polyatomic
molecules, the relative orientation of the colliding pair will determine the
Franck-Condon factors and whether deeply bound molecules in the ground state
can be formed in significant numbers. The basic assumption is that ultracold
polar molecules aligned by external fields will approach each other very
slowly with a specific relative angle. The otherwise very complex potential
energy surface can be replaced at large separation between two interacting
polar molecules by a simpler set of curves parametrized by the separation $R$
between their center-of-mass, and the angle $\theta$ between their orientation
and the intermolecular axis (see Fig. 1). Since PA into the attractive excited
electronic state occurs at very large separation, the massive molecules do not
accelerate significantly towards each other before spontaneous decay takes
place red, with the molecular partners therefore still well separated and at
large distances. While one of the partners remains in its initial state and
aligned by the external field, the other absorbed a photon and gained angular
momentum leading to its rotation and misalignment of the orientation between
molecular partners. However, since the Franck-Condon principle dictates that
the transition will be most probable where the wave functions overlap best,
decay into geometries with long-range wells in the vicinity of the location of
the ground-excited molecular pair will be favored. As we explain in the next
section, this decay at large separation will lead to the formation of
polyatomic molecules in their electronic ground state, albeit in high excited
ro-vibrational levels. These assumptions allow a drastic simplification of the
full problem, by reducing the full multidimensional problem to a much simpler
one that can described by the intermolecular separation $R$ and the alignment
angle $\theta$.

Byrd et al. Byrd  _et al._ (2012b) found that strong quadrupole interactions
between the molecules create a repulsive barrier in potential energy curves
for parallel and collinear alignment of the diatomic molecules. The barrier,
found to be several Kelvin high, gives rise to outer wells (for collinear
alignments) deep enough to support bound states between two KRb molecules (see
Fig. 1). Byrd et al. Byrd  _et al._ (2012b) also determined that the first two
terms in the long-range potential expansion determine the existence and height
of the barrier for KRb-KRb pairs as well as for several other alkali metal
dimers. For a more general geometry, other terms of the long-range expansion
have to be includedMulder  _et al._ (1979); van der Avoird  _et al._ (1980);
Byrd  _et al._ (2011). In case of KRb-KRb, the barrier was found to survive
deviations up to 20° from the collinear alignment, as well as similar
deviations from a purely planar geometry.

The existence of the barrier in the ground state is serendipitous for the PA
formation of the complex: the long-range potential well in the ground state
will, in general, have a significantly higher probability to be populated by
the spontaneous emission from an electronic excited state accessible to the
PA. Again, this can be understood from the Franck-Condon principleKrems  _et
al._ (2009). In Fig. 1, we depict a sketch for the PA formation of the (XY)2
complex from aligned two polar diatomics XY, with relative orientation
described by the angle $\theta$. Here, the spontaneous emission will populate
deeply bound ro-vibrational levels in the long-range well for nearly-aligned
dimers ($\theta<20$), while only weakly bound and very extended levels (that
would not be stable against collisons and other dissociation processes in the
trap) would be populated in the absence of the barrier.

###  II.3 Approximations and Simplifications

In principle, treating the problem exactly would require the calculation of
multi-dimensional potential energy surfaces for both ground and electronically
excited states for all possible orientations of the approaching molecules, as
well as the corresponding transition dipole moments. In addition, the
molecules need to be aligned using external fields, such as electric fields,
in order to take advantage of the long-range angle-dependent barrier. The
coupling of this external field with the internal degrees of freedom also adds
to the level of complexity to be dealt with.

To render this problem tractable, a first assumption is made that the
individual molecules can be aligned via external fields. In that case, each
pair of molecules is always in the same plane and the angle $\theta$ between
the molecular axis and the line joining the molecules’ center of mass is the
same for both molecules. Relaxing this condition would modify the long-range
interaction Byrd  _et al._ (2012b) and complicate the analysis. However, we
assume a strong enough field so that the alignment is possible, and long-range
barriers and wells persist even if the alignment is not perfect.

Second, in order to further simplify the treatment, we take advantage of the
conditions offered by ultracold temperatures, and the relatively low densities
encountered in experimental settings. The low density allows us to consider
only binary event, i.e. the possibility of having a third molecule in the
vicinity of a pair of molecules can be neglected. The ultralow temperatures
imply very slow moving particles, which is amplified by the considerable mass
of the polar molecules considered, anywhere between roughly 50-150 $u$ or more
(atomic mass unit: $u=1$ g/mol = 1822.8885 $m_{e}$). For example, assuming
$m\sim 100$ $u$, $E\sim k_{B}T\sim\frac{1}{2}mv^{2}$ for temperatures $T\sim
1$ nK - 1 $\mu$K lead to velocities $v\sim 0.4-12.9$ mm/s. Neglecting the
small momentum kick of the photon absorbed during the PA process, a pair of
molecule will travel a very short distance while excited. In fact, the
lifetime of the molecular electronic excited states being in the range of a
few nsec (say 10), and during that time, the molecules cover roughly
$4.08\times 10^{-12}-1.29\times 10^{-10}$ m, i.e. a fraction to a few Bohr
radii 0.08 - 2.4 $a_{0}$ (1 $a_{0}=5.29\times 10^{-11}$ m) . During that time
and distance, the molecules do not have significant opportunity to exert a
torque on each other and their alignment should remain intact.

In addition, as we will describe in the following sections, the PA excitation
takes place in the vicinity of the classical outer turning point of the
excited curve, which is selected to be at large separation by the laser
detuning $\Delta$. At those large distances of several tens of Bohr radii,
each “compact” molecule is well separated from the colliding partner, and does
not change its orientation significantly. Since both molecules are originally
aligned in their electronic ground state, we thus can assume that they
maintain their alignment during the PA process.

However, before the second step of the scheme takes place, i.e. the de-
excitation into the long-range ground state intermolecular well, we need to
account for the molecular rotation induced by the PA process in the first
step. Indeed, once a molecule in the pair of approaching molecules is excited,
it will start rotating, and during the time it is excited (roughly its
lifetime of tens of nsec or so), it will execute many full rotations
(classically speaking), and hence the ground and excited molecules will become
misaligned.

We recall the Franck-Condon principle, which states that the transition will
be most probable where the wave function overlap best. In fact, the Einstein
$A$-coefficient for spontaneous emission of a photon of energy $\hbar\omega$
is proportional to $\omega^{3}|\langle\psi_{g}|D|\psi_{e}\rangle|^{2}$, i.e.
the overlap between the ground ($\psi_{g}$) and excited ($\psi_{e}$) wave
functions will dictate the most probable decay (assuming the transition dipole
moment $D\sim{\rm const.}$ in the appropriate range of separation). Because
the ultracold molecules move a very short distance before decaying to the
ground state (a few $a_{0}$ as described above), the de-excitation will take
place while the molecules are still at large separations. For example,
assuming that the initial PA takes place at the outer turning point of the
barrier (roughly 125 $a_{0}$ for KRb or 45 $a_{0}$ for RbCs: see next
sections), during the roughly 10 nsec in the existed state, the pair should
still be at large separations after de-excitation.

So, while the ground state molecule would remain oriented in the laboratory
frame (e.g., by the electric field), we can thus assume $\theta_{1}=0$ for
convenience. The excited molecule will execute many rotation before decaying,
and the transition itself being vertical, the decay will preserve the angles
when it takes place. Hence the decay will be most probable in tetramer ground
state curves that provide the largest overlap within the few bohr radii of the
barrier of the initial PA step. This should happen when there is a turning
point nearby (i.e. for $\theta_{1}=0$, and any $\theta_{2}$ and $\phi$ leading
to a barrier nearby). For any other orientation of the pair of molecules, the
tetramer surface will be repulsive or deeply attractive in that region,
leading to small wave function probability density and overlap in that range,
and thus small $A$-coefficients.

Therefore, the most probable decay will take place into long-range wells, with
a probability given by an average $A$-coefficient over all angles. A more
explicit discussion is given in the next sections.

##  III Methods

###  III.1 Electronic structure calculation

Figure 2: (Color online) Top: Calculated potential energy curves for (KRb)2
(left panels) and (RbCs)2 (right panels) for $\theta=0^{\mathrm{\circ}}$.
Middle: Zoom on the excited electronic state b${}^{1}\mathrm{A}^{\prime}$ for
selected orientations $\theta$ (legend). Bottom: Same as above for the ground
state. The potential barrier and long-range well are present for both
complexes and diminish with increasing angle $\theta$.

To keep our study relevant to ongoing experiments, we carried out detailed
numerical analyses of long-range tetratomic molecules’ formation from
ultracold pairs of KRb and RbCs molecules. KRb was the first heteronuclear
molecule to be cooled to ultracold temperaturesNi  _et al._ (2008) and remains
one of the most studied system in the context of novel physical regimes and
phenomena. Although the exothermal reaction KRb+KRb$\rightarrow$ K2 \+ Rb2
occurs Byrd  _et al._ (2010), restricted geometries allowed to trap and study
those molecules De Miranda  _et al._ (2011); Moses  _et al._ (2015); Bohn  _et
al._ (2017). Recently, a degenerate Fermi gas of KRb molecules was producedDe
Marco  _et al._ (2019). On the other hand, RbCs is interesting because it is
an experimentally accessible non-reactive bosonic molecule. Moreover, RbCs is
the most extended and least polar of bi-alkali molecules Aymar and Dulieu
(2005): its long-range attractive dipole-dipole interaction is less
significant and the repulsive exchange interaction is more significant than
for other bi-alkali molecules, allowing us to analyze a practical limiting
case.

We constructed ab initio electronic potential energy curves (PECs)
$V_{Y}(R,\theta)$ for (KRb)2 and (RbCs)2 complexes in the planar geometry,
where $R$ and $\theta$ are the distance between the center-of-mass of two
diatomics and their relative orientation angle, respectively (see Fig. 1).

Figure 3: (Color online) Transition dipole moments for
X${}^{1}A^{\prime}$-b${}^{1}A^{\prime}$ transitions in (KRb)2 and (RbCs)2
(solid curves for $\theta=0^{\mathrm{o}}$, dashed curves for
$\theta=20^{\mathrm{o}}$) . In case of (RbCs)2, differences between relative
orientations given by $\theta=0^{\mathrm{o}}\dots 20^{\mathrm{o}}$ are less
than 2%.

For (KRb)2, the ground (X${}^{1}A^{\prime}$) and two lowest excited singlet
PECs (b${}^{1}A^{\prime}$, and c${}^{1}A^{\prime}$) were constructed. For
(RbCs)2, we constructed the ground (X${}^{1}A^{\prime}$) and lowest six
excited PECs (a${}^{3}A^{\prime}$, b${}^{3}A^{\prime}$,
a${}^{3}A^{\prime\prime}$, b${}^{3}A^{\prime\prime}$, b${}^{1}A^{\prime}$,
c${}^{1}A^{\prime}$) in order to correctly account for the spin-orbit
coupling. The KRb and RbCs bond lengths were set to 2.2014 a.u. and 4.4271
a.u., respectively. The ab initio points were calculated for 40 values of $R$
and 7 angles, $\theta$={0°,5°,10°,15°,18°,20°,45°,90°}, at CCSD(T) level of
theory Purvis III and Bartlett (1982); Raghavachari  _et al._ (1989); Bartlett
and Musiał (2007) with the Karlsruhe def2-TZVPP basis set Weigend and Ahlrichs
(2005). The inner valence electrons of Rb and Cs were replaced with Stuttgart
small-core relativistic ECP-28 and ECP-46 effective core potentials (ECPs),
respectively Leininger  _et al._ (1996). The excited state energies were
computed using EE-EOM-CCSD (no frozen core) theoryStanton and Bartlett (1993);
Bartlett and Musiał (2007) with the same basis and ECPs as for the ground
state. All calculations were performed in MOLPRO 2010.1Werner  _et al._
(2010). In Fig. 2 (top panels) we show the ab initio PECs for ground and
excited states for both molecules. The long-range van der Waals tails of the
interaction potentials were constructed from the theory presented in Ref. Byrd
_et al._ (2012c) (Eqs. (9-13), pp. 4). The first- and second-order interaction
energies (electrostatic, dispersion, and induction) were evaluated up to 8th
ordercod . Final PECs were constructed separately for each angle $\theta$ by
fitting a spline to the ab initio points and ensuring a smooth transition to
the long-range form. Ground state potentials were joined before the barrier,
typically at internuclear distances between 20 and 30 Bohr radii, to ensure
the correct long-range form. The most important characteristic values of the
potential energy curves are given in Table 1.

Table 1: Characteristic points of potential energy curves for (KRb)2 and (RbCs)2 for selected orientation angles $\theta$. Inner well depth ($V^{\mathrm{in}}_{\mathrm{min}}$), equilibrium distance $r_{\mathrm{eq}}$, outer well depth ($V^{\mathrm{out}}_{\mathrm{min}}$), barrier height ($h_{\mathrm{b}}$), and the inner turning point of the outer well ($R_{Q}$) are given. All values are in atomic units, and $[n]$ stand for $\times 10^{n}$.  (KRb)2 | $V^{\mathrm{in}}_{\mathrm{min}}$ | $r_{\mathrm{eq}}$ | $h_{\mathrm{b}}$ | $V^{\mathrm{out}}_{\mathrm{min}}$ | $R_{Q}$
---|---|---|---|---|---
$\theta=0^{\circ}$ | -1.74 [-3] | 17.21 | 1.95 [-5] | -2.13 [-8] | 102.07
$\theta=20^{\circ}$ | -2.37 [-3] | 16.11 | 1.76 [-7] | -5.68 [-8] | 62.48
$\theta=25^{\circ}$ | -2.37 [-3] | 16.11 |  |  |
$\theta=90^{\circ}$ | -2.63 [-3] | 10.68 |  |  |
(RbCs)2 |  |  |  |  |
$\theta=0^{\circ}$ | -8.45 [-4] | 18.81 | 3.41 [-5] | -8.16 [-8] | 41.59
$\theta=10^{\circ}$ | -9.65 [-4] | 18.55 | 2.03 [-5] | -1.57 [-8] | 38.32
$\theta=20^{\circ}$ | -1.32 [-3] | 17.74 | -2.41 [-6] | -5.94 [-6] | 31.87
$\theta=90^{\circ}$ | -2.96 [-3] | 11.07 |  |  |

Similarly, transition dipole moments for different orientations were
constructed from ab initio values and smoothly interpolated to the asymptotic
values (Fig. 3).

###  III.2 Photoassociation

Using the PECs computed above, we estimate the production rates of (KRb)2 and
(RbCs)2 in the ground electronic state. We assume a two-step process: the
first step is using PA to form tetratomic molecules in the first excited
singlet electronic state, b${}^{1}A^{\prime}$, which is followed by the second
step, the spontaneous relaxation that transfers the population to the ground
electronic state (Fig. 1). For both molecules, the PA step consists of a
single-photon excitation from a pair of diatomics in their ro-vibrational
ground state of the ground electronic X${}^{1}\Sigma^{+}$ state, red-detuned
by $\Delta$ from a pair in which one of the diatomic is in its ground ro-
vibrational level of its excited $2^{1}\Sigma^{+}$ electronic state. In case
of RbCs molecules, as the two diatomics approach each other, the
b${}^{1}A^{\prime}$ undergoes a crossing with the a${}^{3}A^{\prime\prime}$
and b${}^{3}A^{\prime\prime}$ states, which for RbCs corresponds to the single
excitation to the $1^{3}\Pi$ state. Consequently, for RbCs, it is necessary to
perform a $3\times 3$ spin-orbit calculation mixing the b${}^{1}A^{\prime}$,
a${}^{3}A^{\prime\prime}$ and b${}^{3}A^{\prime\prime}$ states. This was
performed by computing the $3\times 3$ spin-orbit matrix using a three-state
multi-configuration self consistent field (MCSCF) Werner and Knowles (1985);
Knowles and Werner (1985) calculation at every ab initio point. As depicted in
Fig. 2, the spin-orbit coupling impacts the b${}^{1}A^{\prime}$ electronic
state for separation $R$ smaller than 24 Bohr radii, much smaller than the
long-range barrier considered to increase the formation rate of tetratomic
molecules, and hence will not affect our results (discussed in the following
sections).

The PA rate coefficient can be expressed asThorsheim  _et al._ (1987)

$K^{\mathrm{PA}}_{b}=\langle
v_{\mathrm{rel}}\sigma^{\mathrm{PA}}_{b}\rangle\;,$ (1)

where $v_{\mathrm{rel}}$ is the relative velocity of the colliding dimers, and
$\sigma^{\mathrm{PA}}_{b}$ is the PA cross section for forming the complex in
the ro-vibrational state $|b\rangle=|v,\ell\rangle$ of the selected excited
electronic state. Here, $v$ and $\ell$ are the vibrational and total orbital
quantum number, respectively, and the brackets $\langle\cdots\rangle$ indicate
averaging over a Maxwellian distribution of initial velocities at the mean gas
temperature $T$. At ultracold temperatures ($s$-wave regime, $\ell=0$) and
assuming a low PA laser intensity $I$, the optimal PA rate coefficient
$K_{b}^{\mathrm{PA}}(T,I)$ is given byNapolitano  _et al._ (1994); Juarros
_et al._ (2006); Gacesa  _et al._ (2016); Juarros  _et al._ (2006):

$K_{b}^{\mathrm{PA}}(T,I)=\frac{8\pi^{3}}{h^{2}c}\frac{I}{Q_{T}}e^{-1/2}|\langle
b|D(R)|\varepsilon,\ell=0\rangle|^{2},$ (2)

where $Q_{T}=(2\pi\mu k_{B}T/h^{2})^{3/2}$ is the translational partition
function, $\mu$ is the reduced mass of the pair, $\varepsilon$ is the
asymptotic relative collision energy, and $D(R)$ is the transition dipole
moment between the initial continuum state $|\varepsilon,\ell=0\rangle$ of the
colliding pair and the target bound state $|b\rangle=|v,\ell=1\rangle$ of the
complex in an electronic excited state. Here, $k_{B}$ and $h$ are Boltzmann
and Planck constants, and $c$ is the speed of light in vacuum.

###  III.3 Formation rate of ground state molecules

As described in Section II.3, we assume that the most probable decay takes
place in the same geometry as the initial PA step. We account for other final
geometries by a coefficient $\eta$ discussed at the end of this section.

The rate coefficient $K_{b,b^{\prime}}(T,I)$, for the formation of tetratomic
molecules in the ro-vibrational state
$|b^{\prime}\rangle=|v^{\prime},\ell^{\prime}\rangle$ of their electronic
ground state by spontaneous emission from the state $|b\rangle$ of an excited
electronic state, is given byJuarros  _et al._ (2006); Gacesa  _et al._ (2016)

$K_{b,b^{\prime}}(T,I)=K_{b}^{\mathrm{PA}}(T,I)r_{b,b^{\prime}}^{(\alpha)}\;,$
(3)

where $r_{b,b^{\prime}}^{(\alpha)}$ is the branching ratio for the spontaneous
radiative emission for the branch $\alpha=\\{R,Q,P\\}$. The branching ratio
can be expressed in terms of Einstein $A$ coefficients weighted by Hönl-London
factors Herzberg (1950):

$r_{b,b^{\prime}}^{(\alpha)}={A_{b,b^{\prime}}^{(\alpha)}}/\left[{\sum_{b^{\prime}}A_{b,b^{\prime}}^{(\alpha)}+\int
A_{b}^{(\alpha)}(\varepsilon^{\prime})\mathrm{d}\varepsilon^{\prime}}\right]\,.$
(4)

The Einstein $A$ coefficients for bound-bound and bound-free transitions,
$A_{b,b^{\prime}}^{(\alpha)}$ and $A_{b}^{(\alpha)}(\mathrm{\varepsilon})$,
respectively, are

$\displaystyle A_{b,b^{\prime}}^{(\alpha)}$ $\displaystyle=$
$\displaystyle\frac{4e^{2}[\omega_{b,b^{\prime}}^{(\alpha)}]^{3}}{3\hbar
c^{3}}W_{J_{b^{\prime}}}^{(\alpha)}|\langle b^{\prime}|D(R)|b\rangle|^{2}$
$\displaystyle A_{b}^{(\alpha)}(\mathrm{\varepsilon})$ $\displaystyle=$
$\displaystyle\frac{4e^{2}\omega_{b}^{(\alpha)}(\varepsilon)^{3}}{3\hbar
c^{3}}W_{J_{b^{\prime}}}^{(\alpha)}|\langle\varepsilon\ell|D(R)|b\rangle|^{2},$
(5)

where $\hbar\omega_{b,b^{\prime}}^{(\alpha)}=|E_{b}-E_{b^{\prime}}|$ and
$\hbar\omega_{b}^{(\alpha)}(\varepsilon)=|E_{b}-E_{\varepsilon\ell}|$, are the
frequencies for the bound-bound and bound-free transition, respectively, while
$D(R)$ is the dipole transition moment between the two electronic states and
$W_{J_{b^{\prime}}}^{(\alpha)}$ are the Hönl-London factors for the branch
$\alpha$Herzberg (1950). The largest contribution to the PA rate coefficients
comes from small values of $\varepsilon$, validating the use of the dipole
transition moment in free-bound transitions in ultracold regime.

In evaluating the above expressions, we included all bound-bound transitions
and bound-free transitions for energies up to $\varepsilon/k_{B}=1$ mK, as in
Refs. Juarros  _et al._ (2006); Gacesa  _et al._ (2016). We assumed that all
optical transitions are $Q$-branch transitions for $\ell=\ell^{\prime}=0$.
This approximation has minimal impact on the production rate coefficients
because $|v\rangle\approx|v,\ell\rangle$ for both $v$ and $v^{\prime}$ states.
Thus, we simplify the notation by dropping the index $\alpha$, i.e.,
$A_{b,b^{\prime}}\equiv A_{b,b^{\prime}}^{(\alpha)}$.

The comparison between different pair alignments and different molecules is
simplified if the rate coefficients are expressed in terms of binding energies
$E_{b}$ and $E_{b^{\prime}}$ of the ro-vibrational levels $b$ and
$b^{\prime}$. In places where the energy notation is used, we replaced the
indices $b$ and $b^{\prime}$ in the rate coefficients with $E_{b}$ and
$E_{b^{\prime}}$, respectively. Moreover, the binding energy $E_{b}$ can
expressed in terms of the PA laser detuning $\Delta$, selected such that
$E_{b}=h\Delta$ (see Fig. 1). In this notation, the PA rate coefficient and
formation rate coefficient are expressed as $K_{\Delta}^{\mathrm{PA}}(T,I)$
and $K_{\Delta,E_{b^{\prime}}}(T,I)$, and the Einstein $A$ coefficients as
$A_{\Delta,E_{b^{\prime}}}$ and $A_{\Delta}(\varepsilon)$, respectively.

The matrix elements in Eqs. (2) and (5) were evaluated numerically by
diagonalizing the radial Schrödinger equation for electronic motion in
$\theta$-dependent interaction potentials (Fig. 2) and using the corresponding
transition dipole moments (Fig. 3). The wave functions and bound state
energies were calculated numerically using mapped Fourier grid method
(MFGR)Kokoouline  _et al._ (1999), to simultaneously obtain bound and quasi-
discretized continuum spectrum. The MFGR calculation was performed
independently for different relative orientations $\theta$ of the pairs (no
couplings between different potential curves), assuming a variable grid step
size determined by the total box size ($R_{\mathrm{max}}=5000$ $a_{0}$), and
the mapping potential numerically evaluated from the local momentum. The
accuracy of the wave functions in the highly oscillatory short-range region
was ensured by requiring at least 20 points per a single period of the wave
function. The continuum wave functions were found to be in excellent agreement
with a calculation performed using the renormalized Numerov methodJohnson
(1978) for the energies greater than 500 nK.

Finally, we account for the change of geometry due to the rotation of the
excited molecule arising from the photon absorption. Omitting the index
$\alpha$, we can rewrite the branching ratio $r^{(\alpha)}_{b,b^{\prime}}$ in
Eq.(4) as $r_{b,b^{\prime}}=A_{b,b^{\prime}}\tau_{b}$ where $\tau_{b}$ is
simply the lifetime of the state $b$ due to decay. In principle, this lifetime
includes all possible final geometries. As before, we can use the detuning
$\Delta$ to label the bound level $b$, and write
$r_{\Delta,b^{\prime}}=A_{\Delta,b^{\prime}}\tau_{\Delta}$. We now account for
the rotation of the excited molecule using an averaged $A$-coefficient
$\bar{A}_{\Delta,b^{\prime}}$, so that we can introduce the averaged branching
ratio $\bar{r}_{\Delta,b^{\prime}}=\bar{A}_{\Delta,b^{\prime}}\tau_{\Delta}$.

We define a coefficient $\eta$ as follows

$\bar{A}_{\Delta,b^{\prime}}=\frac{1}{4\pi}\int d\Omega
A_{\Delta,b^{\prime}}(\Omega)\equiv\eta A_{\Delta,b^{\prime}}(0),$ (6)

where the solid angle $\Omega$ stands for all possible orientations, and
$A_{\Delta,b^{\prime}}(0)$ is the Einstein $A$-coefficient for the fully
aligned case with $\theta_{1}=\theta_{2}=\phi=0$. As we will discuss in the
next section on results, $A_{\Delta,b^{\prime}}(\Omega)$ is sizable only when
a long-range barrier together with a long-range well exist. All other
configurations lead to very poor Franck-Condon overlaps, and much smaller
$A$-coefficients. We can therefore rewrite $\bar{r}_{\Delta,b^{\prime}}=\eta
r_{\Delta,b^{\prime}}(0)$, where $r_{\Delta,b^{\prime}}(0)$ is the ratio
corresponding to fully aligned case with $A_{\Delta,b^{\prime}}(0)$, so that
$K_{b,b^{\prime}}(T,I)$ in Eq.(3) becomes

$\displaystyle\bar{K}_{\Delta,E_{b^{\prime}}}(T,I)$ $\displaystyle=$
$\displaystyle\eta K_{\Delta}^{\mathrm{PA}}(T,I)r_{\Delta,b^{\prime}}(0)\;,$
(7) $\displaystyle=$ $\displaystyle\eta K_{\Delta,E_{b^{\prime}}}(T,I,0)\;,$
(8)

where $K_{\Delta,E_{b^{\prime}}}(T,I,0)$ is the rate coefficient for the fully
aligned case. Thus $\bar{K}_{\Delta,E_{b^{\prime}}}(T,I)$ is the averaged rate
coefficient for the formation of ultracold bound tetratomic (or polyatomic)
molecules.

##  IV Results and discussion

###  IV.1 (KRb)2 and (RbCs)2 tetratomic molecules

####  IV.1.1 Photoassociation rates and transition probablities

Figure 4: (Color online) Top panels: PA rate coefficients
$K_{\Delta}^{\mathrm{PA}}(T_{0},I_{0})$ for the production of (KRb)2 (left)
and (RbCs)2 (right) complexes, evaluated for $\theta=0$° and 20°.
Middle & bottom panels: Einstein $A$ coefficients,
$A_{\Delta,E_{b^{\prime}}}$, for the spontaneous emission (bound-bound) from
the excited $b^{1}\mathrm{A}^{\prime}$ state to the ground state, shown in
energy-notation. Both quantities are expressed as base ten logarithms. Note
that binding energies $E_{b}=h\Delta$ and
$E_{b^{\prime}}=h\nu{{}_{b}^{\prime}}$ given in terms of frequency units.

Using Eqs. (2) and (5), we calculated the PA rates and Einstein $A$
coefficients for angles $\theta=0\dots 45$° for both molecules. The most
interesting relative alignments that characterize two different physical
regimes based on the barrier height are obtained for $\theta=0^{\circ}$
(highest barrier) and $\theta=\theta_{c}\approx 20^{\circ}$, where
$\theta_{c}$ is the angle for which the barrier dips below the asymptotic
threshold (see Fig. 2).

For the two relative orientations, in Fig. 4 (top panels) we show the PA rate
coefficients $K_{\Delta}^{\mathrm{PA}}(T_{0},I_{0})$, calculated assuming the
PA laser intensity $I_{0}=1$ kW/cm2 and average gas temperature $T_{0}=1$
$\mu$K. For (KRb)2 (Fig. 4, top left), the PA rates vary between about 10-10
cm3/s and $10^{-14}$ cm3/s for the detuning $\Delta$ set to about 10 MHz and
$10^{2}$ to $10^{3}$ GHz, respectively. For larger detunings, the PA rate
rapidly drops to negligible values. The relative alignment of the dimers does
not strongly affect the magnitude of the PA rates even though the oscillation
period is extended for larger angles $\theta$ as the binding energies of the
target states are shifted. Notably, for $\theta=20^{\circ}$, the PA rates are
as high as $10^{-14}$ cm3/s for the largest detunings considered,
$\Delta>2.7\times 10^{3}$ GHz, suggesting that deeper ro-vibrational levels in
the excited state could be populated for $\theta\approx\theta_{c}$, in
contrast to $\theta=0^{\circ}$. For (RbCs)2 (Fig. 4, top right), we computed
the PA rates to be at least an order of magnitude smaller, varying between
about 10-11 and $10^{-18}$ cm3/s, as well as to decrease more rapidly with the
detuning. Moreover, the PA rates are about an order of magnitude larger for
the collinear alignment, $\theta=0^{\circ}$, than for $\theta=\theta_{c}$,
suggesting that the potential barrier has a more significant role in this
system.

The results suggest that the long-range excited tetratomic molecules, with
intermolecular distances $R>50$ Bohr, could be photoassociated regardless of
the orientation $\theta$. As expected, the PA rate coefficients are the
largest for the formation of weakly bound ($\Delta<10$ GHz) excited complexes.
The cutoff value of the detuning $\Delta$, beyond which the PA rates are
negligible, increases with the relative orientation $\theta$, even though
smaller angles yield larger overall PA rates.

In order to determine the production rates of ground state tetramers, we
calculated the Einstein $A$ coefficients, $A_{b,b^{\prime}}^{(\alpha)}$, for
all pairs of ro-vibrational levels $|b\rangle$ and $|b^{\prime}\rangle$ in the
ground and excited electronic state, for all considered relative orientations
$\theta$. Again, the results for two relative alignments $\theta=0^{\circ}$
and $\theta\approx\theta_{c}=20^{\circ}$ are given in Fig. 4 (middle and
bottom panels). As expected, the transition probability is the largest on the
diagonal, where the wave functions of the pair of ro-vibrational states
involved in the transition have similar nodal structure and extent, leading to
large Franck-Condon overlaps.

In the case of (KRb)2 (Fig. 4, left column), for the collinear alignment
($\theta=0^{\circ}$), the effect of the long-range barrier in the ground state
manifests as an extended “diffuse crest” on top of the last oscillation on the
diagonal, around $\Delta=10^{2}$ GHz, where the Einstein $A$ coefficients
remain significant for the transitions to the ground-state ro-vibrational
levels with binding energies as large as $E_{b}^{\prime}=10^{4}$ MHz. The
enhanced transition probabilities are due to the increased Franck-Condon
overlap between the excited state and the outer turning point of the long-
range well in the ground state. For $\theta=20^{\circ}$, the Einstein $A$
coefficients remain significant for even larger detunings, up to
$\Delta=10^{3}$ GHz, corresponding to the ro-vibrational levels as deep as
$E_{b}^{\prime}=10^{3}$ MHz. The difference is due to the fact that the outer
potential well becomes deeper and more extended as the barrier height is
reduced (see Fig. 2, bottom panels). Finally, for $\Delta<10$ GHz, the
Einstein $A$ coefficients for the two alignments become very similar.

Figure 5: (Color online) Production rate of ground state molecules as a
function of PA laser detuning $\Delta$ and binding energy $E_{b^{\prime}}$
shown for $\theta=0^{\circ}$ (left) and $\theta=20^{\circ}$ (right). All axes
are in log scale.

####  IV.1.2 Ground-state production rate coefficients

The rate coefficients $K_{\Delta,E_{b^{\prime}}}(T_{0},I_{0})$, computed using
Eq. (3) for the production of (KRb)2 and (RbCs)2 tetratomic molecules in the
electronic ground state, are shown in Fig. 5 for the PA laser detuning
$\Delta$ and the binding energy $E_{b^{\prime}}$. The figure illustrates the
main result of the study and warrants a more detailed discussion. As above,
for both molecules, we analyzed the two relative alignments of the diatomics,
given by $\theta=0^{\circ}$ (top panels) and $\theta=20^{\circ}$ (bottom
panels). Note that the critical angle $\theta_{c}$, for which the potential
barrier dips below the kinetic interaction energy of the pair, is
$\theta_{c}=20.6^{\circ}$ for (KRb)2 vs $\theta_{c}=19^{\circ}$ for (RbCs)2.

The common feature in all cases is the diagonal distribution of distinct
vertically-elongated lobes, where $K_{\Delta,E_{b^{\prime}}}>10^{-16}$
cm-3s-1. The lobes correspond to the bound ro-vibrational levels in the
electronic ground state of the particular case close to the dissociation limit
that have the binding energy $E_{b^{\prime}}$ between 0.5 and $10^{5}$ MHz.
The binding energies corresponding to the bound ro-vibrational states in the
outer well are marked as the purple shaded area; the energies below it
correspond to the weakly bound near-dissociation levels present regardless of
the barrier, and the energies above it to the deeply bound levels in the inner
well.

For (KRb)2 with $\theta=0^{\circ}$ (Fig. 5, top left), we count a total of
nine lobes, four of which correspond to the weakly-bound levels, three at
least partly overlap with the outer well (and likely correspond to the outer
well bound levels), and two correspond to the deeply bound states in the inner
potential well. For $\theta=20^{\circ}$ (Fig. 5, bottom left), the outer well
is deeper and more extended and supports up to seven bound levels. However, no
inner well levels can be populated with a significant probability. In Fig. 6,
we illustrate the probability density $|\psi_{v^{\prime}}(R)|^{2}$ of the ro-
vibrational wave functions $\psi_{v^{\prime}}(R)$ in the outer well
($v^{\prime}=34\dots 38$) for collinear (KRb)2 $(\theta=0^{\circ})$.

In the case of (RbCs)2, the situation is somewhat less favorable for molecule
formation: only the four highest-energy bound levels in the outer well are
accessible for $\theta=0^{\circ}$ (Fig. 5, top right) and none for
$\theta=20^{\circ}$ (Fig. 5, bottom right). This is the case because the outer
well in (RbCs)2 is much deeper (see Fig. 2), and bound levels more localized
but more difficult to populate, according to Franck-Condon principle. For both
molecules, in case of collinear diatomics ($\theta=0^{\circ}$), the production
rates are about hundred times greater than for $\theta=20^{\circ}$. The
dependence on the orientation is the most pronounced for small binding
energies, $E_{b^{\prime}}<50$ MHz in (KRb)2 as well as for
$10<E_{b^{\prime}}<100$ MHz in (RbCs)2, where the rates are negligible for
$\theta=20^{\circ}$.

Figure 6: (Color online) Probability density of the five lowest vibrational
levels ($v^{\prime}=34\dots 38$, indicated on the wave functions) in the outer
potential well of (KRb)2 molecule at collinear alignment $(\theta=0^{\circ})$.
The colors correspond to the production rates given in Fig. 5.

Once the tetratomic molecules are formed in the ground electronic state, they
will spontaneously relax to the absolute ground state via a radiative cascade,
providing they are not destroyed by collisions or lost from the trap through
other mechanisms.

In order to estimate the experimental feasibility of the proposed approach, we
can approximate the number of tetratomic complexes formed per second in a
trapped ultracold molecular gas as

$N=\eta_{g}\eta n^{2}VK_{\Delta,E_{b^{\prime}}}(T_{0},I_{0})\;,$ (9)

where $n$ the gas density, $V$ is the total volume of the gas illuminated by
the PA laser of intensity $I_{0}$, $\eta$ is the fraction decaying into the
right geometry as defined in Eq,(6), and $\eta_{g}$ is the fraction of
initially correctly aligned pairs at the mean gas temperature $T_{0}$ (before
the PA process). Based on the existing experiments with KRb and RbCs, we take
the molecular gas density and the PA volume to be $n=10^{12}$ cm-3, and $V=1$
mm3, respectively. We also assume that the polar molecules are aligned before
undergoing the PA, and set $\eta_{g}=1$.

Before describing the numbers of molecules produced, we first discuss the
value of $\eta$. From Eq.(6), we can write

$\eta=\frac{1}{4\pi}\frac{\int d\Omega
A_{\Delta,b^{\prime}}(\Omega)}{A_{\Delta,b^{\prime}}(0)}.$ (10)

As shown in Fig. 7, the probability density of the ground state wave function
$|\psi_{g}|^{2}$ in the vicinity of the outer well (roughly 40-70 $a_{0}$ for
RbCs+RbCs) is 50-100 times larger in the presence of a long-range barrier/well
than without. Hence, the $A$-coefficient will be 50-100 times larger than for
cases without long-range barriers and wells. Furthermore, the $A$-coefficient
does not change significantly as long as the long-range barrier and an outer
well exist, and so in these cases we can assume
$A_{\Delta,b^{\prime}}(\Omega)\approx A_{\Delta,b^{\prime}}(0)$. Consequently,
using the long-range expansion (see Eq.(12) in the next section), we can
estimate $\eta$ by integrating over all geometries for which a long-range
barrier and an outer well exist:

$\eta\approx\frac{1}{4\pi}\int d\phi
d(\cos\theta_{1})d(\cos\theta_{2})f(\theta_{1},\theta_{2},\phi),$ (11)

where $f=1$ if a long-range barrier and well exist, and $f=0$ otherwise. For
the molecules considered, the function $f$ is non-zero in a region roughly
defined by $\phi\in[0,20^{o}]$, and $\theta_{1}+\theta_{2}<50^{o}$. This is
illustrated for KRb+KRb in Fig. 2 of Ref. Byrd  _et al._ (2012b), where the
barrier height is given as a function of the angles $\theta_{1}$,
$\theta_{2}$, and $\phi$.

Figure 7: (Color online) Comparison of the probability density of the ground
state wave function $|\psi_{g}|^{2}$ for different orientations
($\theta_{1}=\theta_{2}=\theta$ and $\phi=0$ corresponding to geometries with
($\theta=0$) and without ($\theta=45^{o}$) long-range barrier/wells for
RbCs+RbCs. The dashed box approximately indicates the outer well and the
region where the transitions take place. The $|\psi_{g}|^{2}$ in that region
will be about 5 to 100 times larger for the orientations for which the long-
range barrier is present.

Using the above parameter ranges, the integral evaluates to $\eta=0.0057$,
indicating that about 0.6% of tetramers will decay to a bound ro-vibrational
state in the outer well. If we assume a conservative estimate of the
production rates, where $\eta=0.005$ and
$K_{\Delta,E_{b^{\prime}}}(T_{0},I_{0})\approx 10^{-14}$ cm-3, a total of
$N=5\times 10^{4}$ molecules per second would be produced. These numbers
compare favorably with estimates of molecular production rates in ongoing
experiments with ultracold molecules.

We note that if the molecules are in a cloud (and not a restrictive trapping
geometry), although perfectly aligned by an external field (i.e.
$\theta_{1}=\theta_{2}=\theta$ with $\phi=0$), they can approach each other
from any direction, leading to all possible values $\theta$. Since only
$\theta<20^{o}$ leads to long-range wells, averaging over all $\theta$ gives
$\eta_{g}\approx 0.04$, for which we get $N=2000$ molecules per second.

A similar order-of-magnitude analysis for (RbCs)2 yields $N\approx 10^{4}$
molecules per second (with $\eta_{g}=1$) for four ro-vibrational states in the
outer well of the complex if the pairs are aligned, and $N\approx 400$ s-1 in
the absence of pair alignment ($\eta_{g}=0.04$).

Finally, in the case of KRb, our calculation for $\theta=0^{\circ}$ suggests
that deeply bound states of a tetratomic in a short-range well could be
produced if the PA laser is red-detuned by more than 60 GHz from the
dissociation limit of the excited ground state. While the computed production
rates, $K_{\Delta,E_{b^{\prime}}}(T_{0},I_{0})<10^{-14}$ cm-3 would be on the
low side, the produced tetratomics would be bound by more than 40 GHz, making
them stable against collisional destruction in the trap.

###  IV.2 Other systems

The results obtained for KRb and RbCs can be generalized to other tetratomics
composed of two dipolar molecules, as well as to larger polyatomic complexes,
if their electronic structure supports a long-range potential well. As we
discuss below, this can be accomplished by aligning polar molecules with
appropriate properties with preferentially collinear orientations
($\theta<\theta_{c}$) . The existence of a barrier can be inferred from the
long-range intermolecular potential expansion

$V(R,\theta_{1},\theta_{2},\phi)\xrightarrow{R\rightarrow\infty}-\sum_{n}\frac{W_{n}(\theta_{1},\theta_{2},\phi)}{R^{n}}\;,$
(12)

where the angles $\theta_{1}$,$\theta_{2}$, and $\phi$ define the relative
orientation of the pair in a general case, and the functions $W_{n}$ contain
electrostatic and/or dispersion contributions Mulder  _et al._ (1979); Byrd
_et al._ (2011).

For two aligned identical molecules in a plane, i.e.,
$\theta_{1}=\theta_{2}\equiv\theta$ and $\phi=0$, Eq. (12) simplifies toByrd
_et al._ (2012b)

$V(R,\theta)\simeq-\frac{W_{3}}{R^{3}}-\frac{W_{5}}{R^{5}}\;.$ (13)

For collinear molecules ($\theta=0^{\circ}$), the terms $W_{3}$ and $W_{5}$
become $W_{3}=2\mathcal{D}^{2}$ and
$W_{5}=-6\mathcal{Q}^{2}+4\mathcal{D}\mathcal{Q}\simeq-6\mathcal{Q}^{2}$,
where $\mathcal{D}$, $\mathcal{Q}$, and $\mathcal{O}$ are the electrostatic
dipole, quadrupole, and octupole moments, respectively.

Consequently, if the ratio $|\mathcal{D}/\mathcal{Q}|\ll 1$, a long-range
barrier and potential well can be present and sufficiently deep (on the order
of several Kelvin for the studied alkali metals dimers) to support bound
states. This occurs when the long-range attractive $R^{-3}$ dipole-dipole
interaction is overtaken by the shorter-range repulsive $R^{-5}$ quadrupole-
quadrupole interaction at shorter intermolecular distances where the chemical
bonding does not yet dominate.

By setting $V=0$ in Eq. (13) and solving for the intermolecular distance, we
obtain $R_{\mathcal{Q}}\simeq\sqrt{3\mathcal{Q}^{2}/\mathcal{D}^{2}}$, the
distance where the $R^{-5}$ quadrupole repulsive interaction overtakes
$R^{-3}$ dipolar attraction. Similarly, we can estimate the intermolecular
distance where the short-range attraction becomes dominant as
$R_{\mathrm{sr}}=-W_{6}/W_{5}$, where $W_{6}$ is the collinear orientation
termByrd  _et al._ (2012b). For the collinear molecules ($\theta=0$), the
$W_{6}$ term simplifies to $W_{6}=C_{6,0}+4\left(C_{6,1}+C_{6,2}\right)\simeq
C_{6,0}$, to the leading term (the $C_{n,i}$ coefficients include dispersion
and induction contributions Byrd  _et al._ (2011, 2012c, 2012b)). If
$R_{\mathcal{Q}}$ is outside the region where the bonds are strongly perturbed
(i.e., $R_{\mathrm{sr}}\ll R_{\mathcal{Q}}$) and higher terms in Eq. (12) can
be neglected, we expect the barrier to be present and the long-range bond
between the two polar molecules to be possible.

The $C_{6,0}$ term can be approximated using Zemke  _et al._ (2010)

$C_{6,0}\simeq\frac{3}{2}\frac{\alpha_{1}\alpha_{2}\Delta E_{1}\Delta
E_{2}}{(\Delta E_{1}+\Delta E_{2})}\;,$ (14)

where $\alpha_{i}$ and $\Delta E_{i}$ are the static dipole polarizability and
the excitation energy to the first electronic excited state of molecule $i$,
respectively. For identical molecules, Eq.(14) reduces to

$C_{6,0}\simeq\frac{3}{4}\alpha^{2}\Delta E\;.$ (15)

Table 2: ab initio values of the electric dipole ($\mathcal{D}$), quadrupole ($\mathcal{Q}$), and octupole ($\mathcal{O}$) electrostatic moments (defined from the center-of-mass), long-range expansion coefficient $C_{6,0}$, and outer turning point $R_{\mathcal{Q}}$ for sample molecules likely to have an outer well in the ground electronic state. All values are in atomic units. The values for diatomics were adopted from Ref.Byrd  _et al._ (2012c). Molecule |  Ab initio calculation |
---|---|---
Diatomic | $\mathcal{D}$ | $\mathcal{Q}$ | $\mathcal{O}$ | $C_{6,0}$ | $R_{\mathcal{Q}}$
LiNa | 0.20 | 10.07 | -47.33 | 3289 | 95
KRb | 0.23 | 16.99 | -3.16 | 13540 | 126
KCs | 0.75 | 13.00 | -105.70 | 17260 | 38
NaK | 1.12 | 10.56 | -26.54 | 7777 | 19
RbCs | 0.55 | 14.19 | -5.39 | 19210 | 60
Triatomic |  |  |  |  |
LiONa | 0.800 | 19.948 | -69.046 | 238 | 43.2
NaOK | 0.597 | 24.088 | -73.389 | 465 | 69.9
KORb | 0.005 | 23.527 | -115.538 | 571 | 8150
RbOCs | 0.612 | 19.063 | 89.463 | 732 | 53.9
Tetraatomic |  |  |  |  |
HCCCl | 0.193 | 5.214 | 28.284 | 315 | 46.8
HCCF | 0.254 | 3.082 | 23.789 | 113 | 21.0
HCCBr | 0.128 | 6.127 | 38.711 | 424 | 82.9
HCCI | 0.004 | 7.071 | 51.056 | 608 | 3062
Figure 8: (Color online) Fits of long-range intermolecular potential expansion
to the ab initio points for (LiONa)2 (a) and (HCCCl)2 (b) molecules. Outer
well region is shown. The fitting form
$V(R)=-(2\mathcal{D}^{2})/(R^{3})+(6\mathcal{Q}^{2}-4\mathcal{D}\mathcal{O})/R^{5}-W/R^{6}$
was used, where the parameter $W$ contains $W_{6}$ and higher-order terms.
Fitting coefficients are: $\mathcal{D}=0.792175$, $\mathcal{Q}=19.8437$,
$\mathcal{O}=-68.143$, and $W=-9743.84$ for (LiONa)2 and $\mathcal{D}=0.184$,
$\mathcal{Q}=4.9879$, $\mathcal{O}=30.6274$, and $W=581.217$ for (HCCCl)2.
$R_{\mathcal{Q}}$ computed using the fitted values
($R_{\mathcal{Q}}^{\mathrm{fit}}=43.38$ (LiONa) and
$R_{\mathcal{Q}}^{\mathrm{fit}}=46.95$ (HCCCl)) are in agreement with the ab
initio results (given in Table 2). All quantities are in atomic units.

In Table 2, we summarize the moments as well as the $C_{6,0}$ terms obtained
using ab initio methods Byrd  _et al._ (2012b) or the above approximate
expression. We list the bi-alkali molecules with appropriate ${\cal D/Q}$
ratios Byrd  _et al._ (2012b), as well as families of triatomic linear polar
XOY molecules (with X and Y being different alkali atoms), and of tetratomic
linear polar molecules HCCZ (with Z being a halogen atom). For each of them,
we also include the value of $R_{\mathcal{Q}}$. These parameters were computed
at the CCSD(T)/def2-QZVPPD level of theory, and as in Byrd  _et al._ (2012b),
we used the center-of-mass to compute the multipole terms. The ab initio
$C_{6}$ coefficients were constructed using the London equationZemke  _et al._
(2010) with the excitation energies and polarizabilities calculated at the
EOM-CCSD/def2-QZVPPD and MP2/def2-QZVPPD level of theory, respectively.

For all the molecules given in Table 2, $R_{\mathcal{Q}}$ is located at large
enough separation to allow the existence of a long-range well separated from
the shorter range surface by a barrier. We illustrate this in two specific
cases, the triatomic LiONa and tetratomic HCCCl molecules. Fig. 8 shows the
long-range well for those two cases in the collinear geometry ($\theta=0$).
The long-range wells are located at large separation, and fairly shallow
(roughly $10^{-6}$ and $10^{-7}$ hartree, respectively). Note that long-range
expansions fit the ab initio points very well for both molecules.

By orienting polar molecules with the appropriate ${\cal D/Q}$ ratio, it
should be possible to form larger polyatomic molecules, as long as one can
keep them in the ultracold regime. It is worth noting that more flexibility
can be achieved if mixed molecular species are interacting. In fact, assuming
two different aligned linear polar molecules, some terms in the long-range
expansion (12), that would otherwise cancel out for identical molecules,
remain. For example, a $R^{-4}$ term will appear, to give (for $\theta=0$)

$V(R)\simeq-\frac{2{\cal D}_{1}{\cal D}_{2}}{R^{3}}+\frac{3({\cal D}_{1}{\cal
Q}_{2}-{\cal D}_{2}{\cal Q}_{1})}{R^{4}}+\frac{6{\cal Q}_{1}{\cal
Q}_{1}}{R^{5}}-\frac{W_{6}}{R^{6}}\;.$

The sign of the $R^{-4}$ term depends on the relative value of ${\cal D}_{i}$
and ${\cal Q}_{i}$ of each molecule; if ${\cal D}_{1}{\cal Q}_{2}>{\cal
D}_{2}{\cal Q}_{1}$ a longer range barrier could exist, while the reverse
might reduce the size of the barrier or even prevent it. Naturally, mixing
ultracold molecules in the same trap is currently experimentally challenging,
but its realization would provide added flexibility in forming more complex
molecules.

##  V Summary and Conclusions

In this study, we theoretically investigated photoassociative production of
cold tetratomic molecules, and potentially larger polyatomic molecules, from
pairs of ultracold polar molecules. Our working hypothesis was that the
standard PA of ultracold atom pairs into diatomic molecules could be extended
to form a polyatomic complex from two ultracold polar molecules. We
demonstrated the concept quantitatively using two ultracold polar molecules
that have been cooled and trapped successfully, namely KRb and RbCs, to form
(KRb)2 and (RbCs)2 tetratomic molecules. We computed realistic electronic
potentials for specific geometries obtained by aligning these molecules using
external fields. Such alignments and restricted geometries, beside being used
for our proof-of-concept, have been realized in laboratory settings Moses  _et
al._ (2015, 2017). These electronic potential energy curves were computed from
first principles for the first several excited states of (KRb)2 and (RbCs)2
(including spin-orbit coupling), and various relative orientations ($\theta$
between $0^{\circ}$ and $45^{\circ}$). The PECs are shown for angles
$0^{\circ}\leqslant\theta\leqslant 20^{\circ}$ in Fig. 2, while the transition
dipole moments between the ground state and the first excited singlet state
are given in Fig. 3. The ab initio points are available on request.

We found that the PA rate coefficients for the formation of excited (KRb)2 and
(RbCs)2 tetratomics are comparable to the rates reported for the PA of atomic
pairs, as expected, based on our knowledge of long-range potential
interactions and general properties of the PA. We also confirmed that the
radiative decay into the electronic ground state of the tetratomic molecule is
highly impacted by the orientation of the diatomic molecules, as dictated by
the Franck-Condon principle. In fact, since the PA takes place at a large
separation $R$ where the excited electronic energy surface is not steep, the
pair of diatomics does not have sufficient time to significantly accelerate
(due also to their large mass) before spontaneously decaying. This allows one
to consider the decay as a vertical transition down to the ground state. Since
the Franck-Condon overlap is enhanced by the existence of a long-range barrier
and well for the appropriate geometry, the polyatomic ground state molecules
will be mostly produced in the collinear geometry (with $\theta\approx 0$),
with the production rate basically proportional to the rate of fully aligned
pairs with $\theta=0$.

Specifically, if the diatomics approach each other in a collinear geometry
($\theta=0$) or near-collinear ($\theta<\theta_{c}\approx 20^{\circ}$), a
long-range outer potential well in the ground electronic state may be present
for molecules whose dipole-to-quadrupole electrostatic moments ratio is small,
i.e., $|\mathcal{D}/\mathcal{Q}|\ll 1$. The spontaneous emission rates, from
the excited state into bound ro-vibrational states of the outer well, were
found to be significant both for (KRb)2 and (RbCs)2. The resulting ground-
state molecular complex (i.e., (RbCs)2 tetratomic molecules) formed in this
way consists of a pair of polar molecules held together at long range by
electrostatic interactions.

We computed the binding energy of the complexes to be in the range of 30 MHz
to 300 MHz (or about 1.4 mK to 14 mK) for (KRb)2, and 20 MHz to 4 GHz (about 1
mK to 190 mK) in the case of (RbCs)2. The long-range tetratomics would be
considered strongly bound at conditions found in ultracold traps and stable
enough against decay through rotation of the pair (estimated to be in
milliseconds) to allow the experimental verification (e.g., through a
resonantly enhanced multi-photon ionization or a similar technique). The
predicted long-range complex is a novel type of tetratomic molecule, or, more
generally, a long-range complex composed of two nearly-collinear polar
molecules.

We generalized the conclusions based on our quantitative study conducted on
RbCs and KRb to other systems. In fact, we identified other diatomic,
triatomic, and tetratomic linear polar molecules with the appropriate
properties, i.e. $|\mathcal{D}/\mathcal{Q}|\ll 1$. Table 2 lists families of
such molecules, namely XY, XOY, and HCCZ, where X and Y are alkali atoms, and
Z a halogen atom. All those molecules exhibit the right behavior, with the
corresponding $R_{\cal Q}$ at large distances. We confirmed this result by
using ab initio computation and obtained long-range potential wells for the
collinear case of LiONa and HCCCl, in good agreement with the estimated value
of $R_{\cal Q}$ (see Fig. 8). These preliminary calculations suggest that
these systems have long-range wells capable of supporting bound states.

Finally, the production of polyatomic molecules using the long-range well will
result in the formation of polyatomic molecules in highly excited ro-
vibrational levels of the molecular complex. The existence of long-lived
states in those long-range complexes is also made possible by the ultracold
temperatures that prevent their break up due to collisions. This type of
states is usually not accessible in more standard settings, and could allow
spectroscopic study of regimes relevant to roaming reactions Bowman (2014).

###### Acknowledgements.

This work was partially supported by the MURI U.S. Army Research Office grant
number W911NF-14-1-0378 (MG, JS), and by the National Science Foundation grant
number PHY-2034284 (RC).

## References

  * Yan  _et al._ (2013) B. Yan, S. A. Moses, B. Gadway, J. P. Covey, K. R. A. Hazzard, A. M. Rey, D. S. Jin, and J. Ye, Nature (London) 501, 521 (2013).
  * Gorshkov  _et al._ (2013) A. V. Gorshkov, K. R. A. Hazzard, and A. M. Rey, Molecular Physics 111, 1908 (2013), arXiv:1301.5636 [cond-mat.quant-gas] .
  * Eisert  _et al._ (2013) J. Eisert, M. van den Worm, S. R. Manmana, and M. Kastner, Phys. Rev. Lett. 111, 260401 (2013), arXiv:1309.2308 [quant-ph] .
  * Richerme  _et al._ (2014) P. Richerme, Z.-X. Gong, A. Lee, C. Senko, J. Smith, M. Foss-Feig, S. Michalakis, A. V. Gorshkov, and C. Monroe, Nature (London) 511, 198 (2014).
  * Doçaj  _et al._ (2016) A. Doçaj, M. L. Wall, R. Mukherjee, and K. R. A. Hazzard, Phys. Rev. Lett. 116, 135301 (2016), arXiv:1512.06177 [cond-mat.quant-gas] .
  * Kruckenhauser  _et al._ (2020) A. Kruckenhauser, L. M. Sieberer, W. G. Tobias, K. Matsuda, L. De Marco, J.-R. Li, G. Valtolina, A. M. Rey, J. Ye, M. A. Baranov, and P. Zoller, arXiv e-prints , arXiv:2001.11792 (2020), arXiv:2001.11792 [cond-mat.quant-gas] .
  * Capogrosso-Sansone (2011) B. Capogrosso-Sansone, Phys. Rev. A 83, 053611 (2011), arXiv:1009.6213 [cond-mat.stat-mech] .
  * Hazzard  _et al._ (2013) K. R. A. Hazzard, S. R. Manmana, M. Foss-Feig, and A. M. Rey, Phys. Rev. Lett. 110, 075301 (2013), arXiv:1209.4076 [cond-mat.quant-gas] .
  * Wall  _et al._ (2015) M. L. Wall, K. R. A. Hazzard, and A. M. Rey, “Quantum Magnetism with Ultracold Molecules,” in  _From Atomic to Mesoscale: The Role of Quantum Coherence in Systems of Various Complexities_, edited by S. Malinovskaya, A. and et al. (World Scientific Publishing Co, 2015) pp. 3–37.
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Micheli  _et al._ (2006) A. Micheli, G. Brennen, and P. Zoller, Nat. Phys. 2, 341 (2006).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008), arXiv:0704.3011 .
  * DeMille (2002) D. DeMille, Phys. Rev. Lett. 88, 067901 (2002).
  * Krems  _et al._ (2009) R. Krems, B. Friedrich, and W. C. Stwalley,  _Cold molecules: theory, experiment, applications_ (CRC press, 2009).
  * Kuznetsova  _et al._ (2011) E. Kuznetsova, S. F. Yelin, and R. Côté, Quantum Information Processing 10, 821 (2011).
  * Côté (2014) R. Côté, “Ultracold molecules: Their formation and application to quantum computing,” in  _Quantum Information and Computation for Chemistry_ (John Wiley & Sons, Ltd, 2014) pp. 403–448.
  * Yelin  _et al._ (2006) S. F. Yelin, K. Kirby, and R. Côté, Phys. Rev. A 74, 050301(R) (2006).
  * Kuznetsova  _et al._ (2010) E. Kuznetsova, M. Gacesa, S. F. Yelin, and R. Côté, Phys. Rev. A 81, 030301 (2010), arXiv:0908.4558 [quant-ph] .
  * Stuhler  _et al._ (2005) J. Stuhler, A. Griesmaier, T. Koch, M. Fattori, T. Pfau, S. Giovanazzi, P. Pedri, and L. Santos, Phys. Rev. Lett. 95, 150406 (2005).
  * Boisseau  _et al._ (2002) C. Boisseau, I. Simbotin, and R. Côté, Phys. Rev. Lett. 88, 133004 (2002), arXiv:physics/0201022 [physics.atom-ph] .
  * Singer  _et al._ (2005) K. Singer, J. Stanojevic, M. Weidemüller, and R. Côté, J. Phys. B 38, S295 (2005).
  * Jaksch  _et al._ (2000) D. Jaksch, J. Cirac, P. Zoller, S. Rolston, R. Côté, and M. Lukin, Phys. Rev. Lett. 85, 2208 (2000).
  * Lukin  _et al._ (2001) M. D. Lukin, M. Fleischhauer, R. Cote, L. Duan, D. Jaksch, J. I. Cirac, and P. Zoller, Phys. Rev. Lett. 87, 037901 (2001).
  * Cote (2005) R. Cote, in  _Active and Passive Optical Components for WDM Communications V_ , Vol. 6014 (International Society for Optics and Photonics, 2005) p. 601415.
  * Stanojevic  _et al._ (2006) J. Stanojevic, R. Côté, D. Tong, S. Farooqi, E. Eyler, and P. Gould, Eur. Phys. J. D 40, 3 (2006).
  * Hazzard  _et al._ (2014) K. R. A. Hazzard, B. Gadway, M. Foss-Feig, B. Yan, S. A. Moses, J. P. Covey, N. Y. Yao, M. D. Lukin, J. Ye, D. S. Jin, and A. M. Rey, Phys. Rev. Lett. 113, 195302 (2014), arXiv:1402.2354 [cond-mat.quant-gas] .
  * 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).
  * Moses  _et al._ (2015) S. A. Moses, J. P. Covey, M. T. Miecnikowski, B. Yan, B. Gadway, J. Ye, and D. S. Jin, Science 350, 659 (2015), arXiv:1507.02377 [cond-mat.quant-gas] .
  * Park  _et al._ (2015) J. W. Park, S. A. Will, and M. W. Zwierlein, Phys. Rev. Lett. 114, 205302 (2015).
  * 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).
  * 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).
  * Reichsöllner  _et al._ (2017) L. Reichsöllner, A. Schindewolf, T. Takekoshi, R. Grimm, and H.-C. Nägerl, Phys. Rev. Lett. 118, 073201 (2017), arXiv:1607.06536 [cond-mat.quant-gas] .
  * 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), arXiv:1602.03947 [cond-mat.quant-gas] .
  * Heo  _et al._ (2012) M.-S. Heo, T. T. Wang, C. A. Christensen, T. M. Rvachov, D. A. Cotta, J.-H. Choi, Y.-R. Lee, and W. Ketterle, Phys. Rev. A 86, 021602 (2012).
  * Son  _et al._ (2020) H. Son, J. J. Park, W. Ketterle, and A. O. Jamison, Nature 580, 197 (2020).
  * Gaubatz  _et al._ (1990) U. Gaubatz, P. Rudecki, S. Schiemann, and K. Bergmann, J. Chem. Phys. 92, 5363 (1990).
  * Thorsheim  _et al._ (1987) H. R. Thorsheim, J. Weiner, and P. S. Julienne, Phys. Rev. Lett. 58, 2420 (1987).
  * Jones  _et al._ (2006) K. M. Jones, E. Tiesinga, P. D. Lett, and P. S. Julienne, Rev. Mod. Phys. 78, 483 (2006).
  * Pellegrini  _et al._ (2008) P. Pellegrini, M. Gacesa, and R. Côté, Phys. Rev. Lett. 101, 053201 (2008), arXiv:0806.1295 [physics.atom-ph] .
  * Pellegrini and Côté (2009) P. Pellegrini and R. Côté, New Journal of Physics 11, 055047 (2009).
  * Krzyzewski  _et al._ (2015) S. P. Krzyzewski, T. G. Akin, J. Dizikes, M. A. Morrison, and E. R. I. Abraham, Phys. Rev. A 92, 062714 (2015).
  * Hai  _et al._ (2020) Y. Hai, L.-H. Li, J.-L. Li, G.-R. Wang, and S.-L. Cong, The Journal of Chemical Physics 152, 174307 (2020), https://doi.org/10.1063/5.0001794 .
  * Dulieu and Gabbanini (2009) O. Dulieu and C. Gabbanini, Rep. Prog. Phys. 72, 086401 (2009).
  * Balakrishnan (2016) N. Balakrishnan, J. Chem. Phys. 145, 150901 (2016).
  * Byrd  _et al._ (2010) J. N. Byrd, J. A. Montgomery, Jr., and R. Côté, Phys. Rev. A 82, 010502 (2010), arXiv:1003.4514 [physics.chem-ph] .
  * Buchachenko  _et al._ (2012) A. A. Buchachenko, A. V. Stolyarov, M. M. Szczȩśniak, and G. Chałasiński, J. Chem. Phys. 137, 114305 (2012).
  * Byrd  _et al._ (2012a) J. N. Byrd, H. Harvey Michels, J. A. Montgomery, R. Côté, and W. C. Stwalley, J. Chem. Phys. 136, 014306 (2012a).
  * Byrd  _et al._ (2012b) J. N. Byrd, J. A. Montgomery, Jr., and R. Côté, Phys. Rev. Lett. 109, 083003 (2012b), arXiv:1207.2797 [physics.chem-ph] .
  * Byrd  _et al._ (2012c) J. N. Byrd, J. A. Montgomery, Jr., and R. Côté, Phys. Rev. A 86, 032711 (2012c), arXiv:1207.3546 [physics.chem-ph] .
  * Lepers  _et al._ (2010) M. Lepers, O. Dulieu, and V. Kokoouline, Phys. Rev. A 82, 042711 (2010).
  * Pérez-Ríos  _et al._ (2015) J. Pérez-Ríos, M. Lepers, and O. Dulieu, Phys. Rev. Lett. 115, 073201 (2015), arXiv:1505.03288 [physics.atom-ph] .
  * Vexiau  _et al._ (2015) R. Vexiau, M. Lepers, M. Aymar, N. Bouloufa-Maafa, and O. Dulieu, J. Chem. Phys. 142, 214303 (2015), arXiv:1502.05636 [physics.atom-ph] .
  * Lepers  _et al._ (2016) M. Lepers, G. Quéméner, E. Luc-Koenig, and O. Dulieu, J. Phys. B 49, 014004 (2016), arXiv:1508.06066 [physics.atom-ph] .
  * Quéméner (2017) G. Quéméner, ArXiv e-prints (2017), arXiv:1703.09174 [physics.atom-ph] .
  * Schnabel  _et al._ (2021) J. Schnabel, T. Kampschulte, S. Rupp, J. Hecker Denschlag, and A. Köhn, Phys. Rev. A 103, 022820 (2021).
  * Lett  _et al._ (1995) P. D. Lett, P. S. Julienne, and W. D. Phillips, Annu. Rev. Phys. Chem. 46, 423 (1995).
  * Stwalley and Wang (1999) W. C. Stwalley and H. Wang, J. Mol. Spectr. 195, 194 (1999).
  * Côté (2010) R. Côté, in  _Proceedings Of The Dalgarno Celebratory Symposium_ (World Scientific, 2010) pp. 262–280.
  * Fioretti  _et al._ (1998) A. Fioretti, D. Comparat, A. Crubellier, O. Dulieu, F. Masnou-Seeuws, and P. Pillet, Phys. Rev. Lett. 80, 4402 (1998).
  * Vanhaecke  _et al._ (2002) N. Vanhaecke, W. de Souza Melo, B. Laburthe Tolra, D. Comparat, and P. Pillet, Phys. Rev. Lett. 89, 063001 (2002).
  * Vanhaecke  _et al._ (2004) N. Vanhaecke, C. Lisdat, B. T’jampens, D. Comparat, A. Crubellier, and P. Pillet, Eur. Phys. J. D 28, 351 (2004).
  * Côté and Dalgarno (1997) R. Côté and A. Dalgarno, Chem. Phys. Lett. 279, 50 (1997).
  * Côté and Dalgarno (1999) R. Côté and A. Dalgarno, J. Mol. Spectrosc. 195, 236 (1999).
  * Hu  _et al._ (2014) X.-J. Hu, T. Xie, Y. Huang, and S.-L. Cong, Phys. Rev. A 89, 052712 (2014).
  * Gacesa and Côté (2014) M. Gacesa and R. Côté, J. Mol. Spectrosc. 300, 124 (2014), arXiv:1402.0494 [physics.atom-ph] .
  * Hai  _et al._ (2020) Y. Hai, L.-H. Li, J.-L. Li, G.-R. Wang, and S.-L. Cong, J. Chem. Phys. 152, 174307 (2020).
  * Sun  _et al._ (2020) Z.-X. Sun, Y. Hai, B.-K. Lyu, G.-R. Wang, and S.-L. Cong, Journal of Physics B: Atomic, Molecular and Optical Physics 53, 205204 (2020).
  * Koch  _et al._ (2006) C. P. Koch, E. Luc-Koenig, and F. Masnou-Seeuws, Phys. Rev. A 73, 033408 (2006), arXiv:physics/0508090 [physics.atom-ph] .
  * Ghosal  _et al._ (2009) S. Ghosal, R. J. Doyle, C. P. Koch, and J. M. Hutson, New J. Phys. 11, 055011 (2009), arXiv:0810.5703 [physics.atom-ph] .
  * Carini  _et al._ (2015) J. L. Carini, S. Kallush, R. Kosloff, and P. L. Gould, Phys. Rev. Lett. 115, 173003 (2015), arXiv:1602.08026 [physics.atom-ph] .
  * Wang  _et al._ (2017) M. Wang, J.-L. Li, X.-J. Hu, M.-D. Chen, and S.-L. Cong, Phys. Rev. A 96, 043417 (2017).
  * Gacesa  _et al._ (2013) M. Gacesa, S. Ghosal, J. N. Byrd, and R. Côté, Phys. Rev. A 88, 063418 (2013), arXiv:1310.7140 [physics.atom-ph] .
  * Hu  _et al._ (2015) X.-J. Hu, J.-L. Li, T. Xie, Y.-C. Han, and S.-L. Cong, Phys. Rev. A 92, 032709 (2015).
  * Mulder  _et al._ (1979) F. Mulder, A. van der Avoird, and P. E. Wormer, Mol. Phys. 37, 159 (1979).
  * van der Avoird  _et al._ (1980) A. van der Avoird, P. E. Wormer, F. Mulder, and R. M. Berns, in  _Van der Waals Systems_ (Springer, 1980) pp. 1–51.
  * Byrd  _et al._ (2011) J. N. Byrd, R. Côté, and J. A. Montgomery, J. Chem. Phys. 135, 244307 (2011), arXiv:1109.4410 [physics.chem-ph] .
  * De Miranda  _et al._ (2011) M. De Miranda, A. Chotia, B. Neyenhuis, D. Wang, G. Quéméner, S. Ospelkaus, J. Bohn, J. Ye, and D. Jin, Nat. Phys. 7, 502 (2011).
  * Bohn  _et al._ (2017) J. L. Bohn, A. M. Rey, and J. Ye, Science 357, 1002 (2017), arXiv:1708.02806 [physics.atom-ph] .
  * De Marco  _et al._ (2019) L. De Marco, G. Valtolina, K. Matsuda, W. G. Tobias, J. P. Covey, and J. Ye, Science 363, 853 (2019).
  * Aymar and Dulieu (2005) M. Aymar and O. Dulieu, J. Chem. Phys. 122, 204302 (2005), quant-ph/0502059 .
  * Purvis III and Bartlett (1982) G. D. Purvis III and R. J. Bartlett, J. Chem. Phys. 76, 1910 (1982).
  * Raghavachari  _et al._ (1989) K. Raghavachari, G. W. Trucks, J. A. Pople, and M. Head-Gordon, Chem. Phys. Lett. 157, 479 (1989).
  * Bartlett and Musiał (2007) R. J. Bartlett and M. Musiał, Rev. Mod. Phys. 79, 291 (2007).
  * Weigend and Ahlrichs (2005) F. Weigend and R. Ahlrichs, Phys. Chem. Chem. Phys. 7, 3297 (2005).
  * Leininger  _et al._ (1996) T. Leininger, A. Nicklass, W. Küchle, H. Stoll, M. Dolg, and A. Bergner, Chem. Phys. Lett. 255, 274 (1996).
  * Stanton and Bartlett (1993) J. F. Stanton and R. J. Bartlett, J. Chem. Phys. 98, 7029 (1993).
  * Werner  _et al._ (2010) H.-J. Werner, P. J. Knowles, F. R. Manby, M. Schütz, _et al._ , “MOLPRO, version 2010.1, a package of ab initio programs,” (2010), see http://www.molpro.net.
  * (89) “The sum in eq. (9) was evaluated up to n=8. the code to evaluate long-range potentials is available on request,” .
  * Werner and Knowles (1985) H.-J. Werner and P. J. Knowles, J. Chem. Phys. 82, 5053 (1985).
  * Knowles and Werner (1985) P. J. Knowles and H. J. Werner, Chem. Phys. Lett. 115, 259 (1985).
  * Napolitano  _et al._ (1994) R. Napolitano, J. Weiner, C. J. Williams, and P. S. Julienne, Phys. Rev. Lett. 73, 1352 (1994).
  * Juarros  _et al._ (2006) E. Juarros, K. Kirby, and R. Côté, J. Phys. B 39, 965 (2006).
  * Gacesa  _et al._ (2016) M. Gacesa, J. A. Montgomery, H. H. Michels, and R. Côté, Phys. Rev. A 94, 013407 (2016), arXiv:1603.08032 [physics.atom-ph] .
  * Juarros  _et al._ (2006) E. Juarros, P. Pellegrini, K. Kirby, and R. Côté, Physical Review A 73, 041403 (2006).
  * Herzberg (1950) G. Herzberg,  _New York: Van Nostrand Reinhold, 1950, 2nd ed._ (1950).
  * Kokoouline  _et al._ (1999) V. Kokoouline, O. Dulieu, R. Kosloff, and F. Masnou-Seeuws, J. Chem. Phys. 110, 9865 (1999).
  * Johnson (1978) B. R. Johnson, J. Chem. Phys. 69, 4678 (1978).
  * Zemke  _et al._ (2010) W. T. Zemke, J. N. Byrd, H. H. Michels, J. A. Montgomery, and W. C. Stwalley, J. Chem. Phys. 132, 244305 (2010).
  * Moses  _et al._ (2017) S. A. Moses, J. P. Covey, M. T. Miecnikowski, D. S. Jin, and J. Ye, Nat. Phys. 13, 13 (2017), arXiv:1610.07711 [cond-mat.quant-gas] .
  * Bowman (2014) J. M. Bowman, Molecular Physics 112, 2516 (2014).