Photoassociation of ultracold long-range polyatomic molecules

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 polyatomic molecules.


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 fewand many-body quantum physics. For example, they could be used to engineer and study lattice spin models in strongly interacting many-body hamiltonians [1][2][3][4][5][6], supersolidity [7], unconventional superfluid phases and quantum magnetism [8,9], to name a few. Moreover, the tunability of interactions combined with the extraordinary degree of control available in current experiments [10] make such systems suitable for engineering quantum simulators [11,12], as well as quantum entanglement and information processing [11,[13][14][15][16][17][18]. 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 dipoles [19] or between excited Rydberg atoms [20,21], 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 fields [22][23][24][25][26].
Producing ultracold polar molecules was achieved for only a handful of diatomic species and up to this day remains challenging. Fermionic 40 K 87 Rb [27,28] were the first heteronuclear molecules to be produced in deeply * marko.gacesa@ku.ac.ae † robin.cote@uconn.edu bound molecular states, with NaK [29,30], RbCs [31][32][33], NaRb [34], and LiNa [35,36] 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) [37] was used to transfer the population into the molecular ground state. Alternative approaches to produce ultracold molecules is the photoassociation (PA) of ultracold atomic pairs [38,39] or Feshbach Optimized PA (FOPA) [40][41][42][43].
Thus far, those techniques were limited to diatomic molecules, even though the theory suggests that they could be extended to larger molecules [44,45]. 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 → K 2 +Rb 2 is exothermic) or not well-understood "complex nature of scattering processes" found to occur even for endothermic intermolecular reactions [33]. 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 [46,47]. Byrd et al. [48] 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 [49,50]. 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 arXiv:2007.00261v2 [physics.atom-ph] 10 Jun 2021 critical for experiments' design and selection of suitable atomic and molecular systems. Thus, preliminary theoretical studies have been focused on understanding longrange atom-diatom [51,52] and diatom-diatom interactions [48][49][50][53][54][55] for experimentally accessible species including Cs-Cs 2 , KRb-KRb, and other alkali metal diatomic molecules. Notably, Perez-Rios et al. [52] computed the rate coefficient for the PA of Cs 3 and found that it would be comparable to the formation rates observed for diatomic cases, and, more recently, Schnabel et al. [56] theoretically investigated the prospects for photoassociation (PA) of Rb 3 .
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 longrange 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.

A. Background
The photoassociation (PA) of a pair of ultracold atoms was initially proposed by Thorsheim et al. [38] 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 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 rovibrational 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 articles [14,39,44,45,[57][58][59].
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 emission [57,[60][61][62]. The main difficulty with the PA formation can be traced to the Franck-Condon principle, according to (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 θ. The excited tetramer in the ro-vib. state |b (wave function ψ b ), detuned by ∆ from its dissociation energy, is created by absorbing a photon of energy ω from the laser field of intensity I. The excited state spontaneously decays into the ground state, ro-vib. level |b with binding energy E b (wave function ψ b ). Inset: The longrange region of the ground state. Long-range potential well (shown for θ=0°) capable of supporting several bound states exists for θ < θc, where θc ≈ 20°for alkali metal pairs. 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 [63,64]. Another approach to enhancing photoassociative production efficiency of stable ultracold molecules can be achieved via Feshbach resonances [40,[65][66][67][68], short laser pulses [65,[69][70][71][72], or both [73,74].

B. Photoassociation in specific geometries
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 θ 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 θ.
Byrd et al. [49] 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. [49] 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 included [75][76][77]. 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 principle [14]. 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 θ. Here, the spontaneous emission will populate deeply bound rovibrational levels in the long-range well for nearly-aligned dimers (θ < 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.

C. Approximations and Simplifications
In principle, treating the problem exactly would require the calculation of multi-dimensional potential en-ergy 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 θ 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 [49] and complicate the analysis. However, we assume a strong enough field so that the alignment is possible, and longrange 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 ∼ 100 u, E ∼ k B T ∼ 1 2 mv 2 for temperatures T ∼ 1 nK -1 µK lead to velocities v ∼ 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 × 10 −12 − 1.29 × 10 −10 m, i.e. a fraction to a few Bohr radii 0.08 -2.4 a 0 (1 a 0 = 5.29 × 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 ∆. 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 ω is proportional to ω 3 | ψ g |D|ψ e | 2 , i.e. the overlap between the ground (ψ g ) and excited (ψ e ) wave functions will dictate the most probable decay (assuming the transition dipole moment D ∼ 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 θ 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 θ 1 = 0, and any θ 2 and φ 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.

A. Electronic structure calculation
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 temperatures [27] and remains one of the most studied system in the context of novel physical regimes and phenomena. Although the exothermal reaction KRb+KRb→ K 2 + Rb 2 occurs [46], restricted geometries allowed to trap and study those molecules [28,78,79]. Recently, a degenerate Fermi gas of KRb molecules was produced [80]. 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 [81]: 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, θ) for (KRb) 2 and (RbCs) 2 complexes in the planar geometry, where R and θ are the distance between the center-of-mass of two diatomics and their relative orientation angle, respectively (see Fig. 1). For (KRb) 2 , the ground (X 1 A ) and two lowest excited singlet PECs (b 1 A , and c 1 A ) were constructed. For (RbCs) 2 , we constructed the ground (X 1 A ) and lowest six excited PECs ( 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, θ={0°,5°,10°,15°,18°,20°,45°,90°}, at CCSD(T) level of theory [82][83][84] with the Karlsruhe def2-TZVPP basis set [85]. 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 [86]. The excited state energies were computed using EE-EOM-CCSD (no frozen core) theory [84,87] with the same basis and ECPs as for the ground state. All calculations were performed in MOLPRO 2010.1 [88]. In Fig. 2 (top panels) we show the ab initio PECs for ground and excited states for both molecules. The longrange van der Waals tails of the interaction potentials were constructed from the theory presented in Ref. [50] (Eqs. (9-13), pp. 4). The first-and second-order interaction energies (electrostatic, dispersion, and induction) were evaluated up to 8th order [89]. Final PECs were constructed separately for each angle θ 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 longrange form. The most important characteristic values of the potential energy curves are given in Table I.
Similarly, transition dipole moments for different orientations were constructed from ab initio values and smoothly interpolated to the asymptotic values (Fig. 3).

B. 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 , which is followed by the second step, the spontaneous relaxation that transfers the population to the ground electronic state ( Angle θ 1). For both molecules, the PA step consists of a singlephoton excitation from a pair of diatomics in their rovibrational ground state of the ground electronic X 1 Σ + state, red-detuned by ∆ from a pair in which one of the diatomic is in its ground ro-vibrational level of its excited 2 1 Σ + electronic state. In case of RbCs molecules, as the two diatomics approach each other, the b 1 A undergoes a crossing with the a 3 A and b 3 A states, which for RbCs corresponds to the single excitation to the 1 3 Π state. Consequently, for RbCs, it is necessary to perform a 3 × 3 spin-orbit calculation mixing the b 1 A , a 3 A and b 3 A states. This was performed by computing the 3 × 3 spin-orbit matrix using a three-state multi-configuration self consistent field (MCSCF) [90,91] calculation at every ab initio point. As depicted in Fig. 2, the spin-orbit coupling impacts the b 1 A 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 as [38] K where v rel is the relative velocity of the colliding dimers, and σ PA b is the PA cross section for forming the complex in the ro-vibrational state |b = |v, of the selected excited electronic state. Here, v and are the vibrational and total orbital quantum number, respectively, and the brackets · · · indicate averaging over a Maxwellian distribution of initial velocities at the mean gas temperature T . At ultracold temperatures (s-wave regime, = 0) and assuming a low PA laser intensity I, the optimal PA rate coefficient K PA b (T, I) is given by [92][93][94][95]: where Q T = (2πµk B T /h 2 ) 3/2 is the translational partition function, µ is the reduced mass of the pair, ε is the asymptotic relative collision energy, and D(R) is the transition dipole moment between the initial continuum state |ε, = 0 of the colliding pair and the target bound state |b = |v, = 1 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.
C. Formation rate of ground state molecules As described in Section II C, 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 η discussed at the end of this section.
The rate coefficient K b,b (T, I), for the formation of tetratomic molecules in the ro-vibrational state |b = |v , of their electronic ground state by spontaneous emission from the state |b of an excited electronic state, is given by [93,94] where r (α) b,b is the branching ratio for the spontaneous radiative emission for the branch α = {R, Q, P }. The branching ratio can be expressed in terms of Einstein A coefficients weighted by Hönl-London factors [96]: The Einstein A coefficients for bound-bound and boundfree transitions, A where ω are the Hönl-London factors for the branch α [96]. The largest contribution to the PA rate coefficients comes from small values of ε, 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 ε/k B = 1 mK, as in Refs. [93,94]. We assumed that all optical transitions are Q-branch transitions for = = 0. This approximation has minimal impact on the production rate coefficients because |v ≈ |v, for both v and v states. Thus, we simplify the notation by dropping the index α, i.e., 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 of the ro-vibrational levels b and b . In places where the energy notation is used, we replaced the indices b and b in the rate coefficients with E b and E b , respectively. Moreover, the binding energy E b can expressed in terms of the PA laser detuning ∆, selected such that E b = h∆ (see Fig. 1). In this notation, the PA rate coefficient and formation rate coefficient are expressed as K PA ∆ (T, I) and K ∆,E b (T, I), and the Einstein A coefficients as A ∆,E b and A ∆ (ε), respectively.
The matrix elements in Eqs. (2) and (5) were evaluated numerically by diagonalizing the radial Schrödinger equation for electronic motion in θ-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) [97], to simultaneously obtain bound and quasi-discretized continuum spectrum. The MFGR calculation was performed independently for different relative orientations θ of the pairs (no couplings between different potential curves), assuming a variable grid step size determined by the total box size (R 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 method [98] 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 α, we can rewrite the branching ratio r 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 ∆ to label the bound level b, and write r ∆,b = A ∆,b τ ∆ . We now account for the rotation of the excited molecule using an averaged A-coefficientĀ ∆,b , so that we can introduce the averaged branching ratior ∆,b =Ā ∆,b τ ∆ .
We define a coefficient η as follows where the solid angle Ω stands for all possible orientations, and A ∆,b (0) is the Einstein A-coefficient for the fully aligned case with θ 1 = θ 2 = φ = 0. As we will discuss in the next section on results, A ∆,b (Ω) 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 Acoefficients. We can therefore rewriter ∆,b = ηr ∆,b (0), where r ∆,b (0) is the ratio corresponding to fully aligned case with A ∆,b (0), so that K b,b (T, I) in Eq.(3) becomes where K ∆,E b (T, I, 0) is the rate coefficient for the fully aligned case. ThusK ∆,E b (T, I) is the averaged rate coefficient for the formation of ultracold bound tetratomic (or polyatomic) molecules.

Photoassociation rates and transition probablities
Using Eqs. (2) and (5), we calculated the PA rates and Einstein A coefficients for angles θ = 0 . . . 45°for both molecules. The most interesting relative alignments that characterize two different physical regimes based on the Middle & bottom panels: Einstein A coefficients, A∆,E b , for the spontaneous emission (bound-bound) from the excited b 1 A state to the ground state, shown in energy-notation. Both quantities are expressed as base ten logarithms. Note that binding energies E b = h∆ and E b = hν b given in terms of frequency units.
barrier height are obtained for θ = 0 • (highest barrier) and θ = θ c ≈ 20 • , where θ 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 PA ∆ (T 0 , I 0 ), calculated assuming the PA laser intensity I 0 = 1 kW/cm 2 and average gas temperature T 0 = 1 µK. For (KRb) 2 (Fig. 4, top left), the PA rates vary between about 10 −10 cm 3 /s and 10 −14 cm 3 /s for the detuning ∆ 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 θ as the binding energies of the target states are shifted. Notably, for θ = 20 • , the PA rates are as high as 10 −14 cm 3 /s for the largest detunings considered, ∆ > 2.7 × 10 3 GHz, suggesting that deeper ro-vibrational levels in the excited state could be populated for θ ≈ θ c , in contrast to θ = 0 • . 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 cm 3 /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, θ = 0 • , than for θ = θ 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 ori-entation θ. As expected, the PA rate coefficients are the largest for the formation of weakly bound (∆ < 10 GHz) excited complexes. The cutoff value of the detuning ∆, beyond which the PA rates are negligible, increases with the relative orientation θ, 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 , for all pairs of ro-vibrational levels |b and |b in the ground and excited electronic state, for all considered relative orientations θ. Again, the results for two relative alignments θ = 0 • and θ ≈ θ c = 20 • 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 (θ = 0 • ), 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 ∆ = 10 2 GHz, where the Einstein A coefficients remain significant for the transitions to the groundstate ro-vibrational levels with binding energies as large as E b = 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 θ = 20 • , the Einstein A coefficients remain significant for even larger detunings, up to ∆ = 10 3 GHz, corresponding to the ro-vibrational levels as deep as E b = 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 ∆ < 10 GHz, the Einstein A coefficients for the two alignments become very similar.

Ground-state production rate coefficients
The rate coefficients K ∆,E b (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 ∆ and the binding energy E b . 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 θ = 0 • (top panels) and θ = 20 • (bottom panels). Note that the critical angle θ c , for which the potential barrier dips below the kinetic interaction energy of the pair, is θ c = 20.6 • for (KRb) 2 vs θ c = 19 • for (RbCs) 2 .
The common feature in all cases is the diagonal distribution of distinct vertically-elongated lobes, where K ∆,E b > 10 −16 cm −3 s −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 between 0.5 and 10 5 MHz. The binding energies corresponding to the bound rovibrational 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 θ = 0 • (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 θ = 20 • (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 |ψ v (R)| 2 of the ro-vibrational wave functions ψ v (R) in the outer well (v = 34 . . . 38) for collinear (KRb) 2 (θ = 0 • ).
In the case of (RbCs) 2 , the situation is somewhat less favorable for molecule formation: only the four highestenergy bound levels in the outer well are accessible for θ = 0 • (Fig. 5, top right) and none for θ = 20 • (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 (θ = 0 • ), the production rates are about hundred times greater than for θ = 20 • . The dependence on the orientation is the most pronounced for small binding energies, E b < 50 MHz in (KRb) 2 as well as for 10 < E b < 100 MHz in (RbCs) 2 , where the rates are negligible for θ = 20 • . 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 where n the gas density, V is the total volume of the gas illuminated by the PA laser of intensity I 0 , η is the fraction decaying into the right geometry as defined in Eq,(6), and η 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 mm 3 , respectively. We also assume that the polar molecules are aligned before undergoing the PA, and set η g = 1.
Before describing the numbers of molecules produced, we first discuss the value of η. From Eq.(6), we can write As shown in Fig. 7, the probability density of the ground state wave function |ψ 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 The |ψg| 2 in that region will be about 5 to 100 times larger for the orientations for which the long-range barrier is present.
. Consequently, using the long-range expansion (see Eq. (12) in the next section), we can estimate η by integrating over all geometries for which a long-range barrier and an outer well exist: 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 φ ∈ [0, 20 o ], and θ 1 + θ 2 < 50 o . This is illustrated for KRb+KRb in Fig. 2 of Ref. [49], where the barrier height is given as a function of the angles θ 1 , θ 2 , and φ.
Using the above parameter ranges, the integral evaluates to η = 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 η = 0.005 and K ∆,E b (T 0 , I 0 ) ≈ 10 −14 cm −3 , a total of N = 5 × 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. θ 1 = θ 2 = θ with φ = 0), they can approach each other from any direction, leading to all possible values θ. Since only θ < 20 o leads to long-range wells, averaging over all θ gives η g ≈ 0.04, for which we get N = 2000 molecules per second.
A similar order-of-magnitude analysis for (RbCs) 2 yields N ≈ 10 4 molecules per second (with η g = 1) for four ro-vibrational states in the outer well of the complex if the pairs are aligned, and N ≈ 400 s −1 in the absence of pair alignment (η g = 0.04).
Finally, in the case of KRb, our calculation for θ = 0 • 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 ∆,E b (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.

B. 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 (θ < θ c ) . The existence of a barrier can be inferred from the long-range intermolecular potential expansion where the angles θ 1 ,θ 2 , and φ define the relative orientation of the pair in a general case, and the functions W n contain electrostatic and/or dispersion contributions [75,77]. For two aligned identical molecules in a plane, i.e., θ 1 = θ 2 ≡ θ and φ = 0, Eq. (12) simplifies to [49] For collinear molecules (θ = 0 • ), the terms W 3 and W 5 become W 3 = 2D 2 and W 5 = −6Q 2 + 4DQ −6Q 2 , where D, Q, and O are the electrostatic dipole, quadrupole, and octupole moments, respectively. Consequently, if the ratio |D/Q| 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 Q 3Q 2 /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 shortrange attraction becomes dominant as R sr = −W 6 /W 5 , where W 6 is the collinear orientation term [49]. For the collinear molecules (θ = 0), the W 6 term simplifies to W 6 = C 6,0 + 4 (C 6,1 + C 6,2 ) C 6,0 , to the leading term (the C n,i coefficients include dispersion and induction contributions [49,50,77]). If R Q is outside the region where the bonds are strongly perturbed (i.e., R sr R 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 [99] C 6,0 3 2 where α i and ∆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 In Table II, we summarize the moments as well as the C 6,0 terms obtained using ab initio methods [49] or the above approximate expression. We list the bi-alkali molecules with appropriate D/Q ratios [49], 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 Q . These parameters were computed at the CCSD(T)/def2-QZVPPD level of theory, and as in [49], we used the center-of-mass to compute the multipole terms. The ab initio C 6 coefficients were constructed using the London  Table II). All quantities are in atomic units. equation [99] 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 II, R 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 (θ = 0). The long-range wells are located at large separation, and fairly shallow (roughly 10 −6 and 10 −7 hartree, respectively). Note that longrange expansions fit the ab initio points very well for both molecules.
By orienting polar molecules with the appropriate 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 θ = 0) The sign of the R −4 term depends on the relative value of D i and Q i of each molecule; if D 1 Q 2 > D 2 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 [28,100]. 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 (θ between 0 • and 45 • ). The PECs are shown for angles 0 • θ 20 • 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 θ ≈ 0), with the production rate basically proportional to the rate of fully aligned pairs with θ = 0.
Specifically, if the diatomics approach each other in a collinear geometry (θ = 0) or near-collinear (θ < θ c ≈ 20 • ), 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., |D/Q| 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 longrange 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 nearlycollinear 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. |D/Q| 1. Table II 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 Q at large distances. We confirmed this result by using ab initio computation and obtained longrange potential wells for the collinear case of LiONa and HCCCl, in good agreement with the estimated value of R 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 [101].