The rich internal structure of ultracold molecules compared to ultracold atoms lend themselves to many applications in emerging areas of quantum science. Ultracold paramagnetic molecules such as Calcium monofluoride, CaF, whose electronic ground state is characterized by a 2 Σ + term, have long been considered as a promising candidate for a number of applications, in particular, quantum simulation [1][2][3][4] , quantum information [5][6][7][8][9] , and precision spectroscopy 10 . This is mostly due to the presence of an unpaired electron as its resultant non-zero electric and magnetic moment serves as a convenient experimental handle for extra control [11][12][13] , by means of external fields (e.g. Stark and Zeeman effects). Additionally, these systems also provide a unique opportunity to improve upon the fundamental understanding of atom-molecule 14 and molecule-molecule interactions 15 , dipolar interactions 16,17 and collision-induced chemistry at the ultra-low range of kinetic energies [18][19][20][21][22] . In particular, experimental explorations of collision-induced trap-loss rate of molecules with singlet and triplet spin multiplicities, in ultracold conditions, have been available for a while (for systems such as Rb 2 , NaRb, KRb, CsRb, NaK, LiNa, NaRb). However, such studies are less prevalent for doublet molecules 23 .
Slowing the translational motion of CaF molecules down to the capture velocity of a 800 mK deep magneto-optical trap (MOT) has been recently achieved by Doyle and co-workers 24 . This follows similar success with SrF, to our knowledge, the first such molecule to be trapped in a MOT [25][26][27][28] . The original work of Lu et al., since then improved to sub-Doppler temperatures [29][30][31][32][33] , represents an important milestone after the seminal work of Di Rosa 11 , the first to observe that molecules such as CaF, CaH, CaOH, SrF, SrOH, YbF, may possess a rovibrational internal structure with a large one-photon oscillator strength and highly diagonal Franck-Condon factors. This, in turn, unlocks the possibility of light-assisted closed cycling transitions, similar to the laser cooling techniques applied to atoms and atomic ions 34 .
Once CaF molecules in the electronic ground state are properly trapped, as demonstrated by Lu et al., a natural next step is the design and implementation of cold collisions between CaF molecules and, say, co-trapped laser cooled atoms or another CaF molecule. The latter case has been recently realized in a pioneering experiment, in which CaF molecules are loaded from a MOT into optical tweezers and, by varying the relative position of two tweezers, CaF + CaF ultracold collisions have been observed to produce two-body loss, most likely due to yet undetermined chemical reactions, with magnitude comparable to a theoretical universal loss rate 23,35 . The former case of CaF collisions with laser cooled atoms is yet an open prospect and, as we shall see below, one of the underlying motivation of the present work. Among the Alkali metal candidates, whose laser cooling and trapping techniques are nowadays routine procedures, only Li 2 S combined with the electronic ground state of CaF provides an exothermicity of about -4440 cm -1 36 . Other atomic species such as Na, K, Rb and Cs would require a few thousand wavenumber of collision-induced excitation in order to trigger chemical events 36 . However, due to the low (< 1 K) kinetic energies involved in such experiments, in general, these collision-induced excitation are all but forbidden energetically.
Thus, the prospects of Li 2 S + CaF X 2 Σ + ultracold reactive collisions are highly regarded as an opportunity to study cold chemistry as well as collision-induced trap losses due to chemical events. To understand and to establish the limits for sympathetic cooling of CaF X 2 Σ + toward even lower temperatures by means of soft collisions with a Li 2 S coolant buffer 37 as well as collisional shielding 38,39 a detailed investigation of Li + CaF collisions is needed. Prior studies of Li 2 S + CaF X 2 Σ + collisions explored only elastic and (non-reactive) inelastic collisions using model potentials or interaction potentials that do not describe the reactive regions. Foreseeing an upcoming demand for more theoretical support in regard to this matter, in this work, we tackle the challenging task of describing a new LiCaF global potential energy surface (PES) and to perform the first description of the Li 2 S + CaF X 2 Σ + -→ Ca 1 S + LiF X 1 Σ + collisions resorting to state-of-the-art quantum reactive scattering, i.e. a coupled-channel (CC) method. It is worthwhile to note that a novel full six-dimensional PES intended to describe the even more complicated CaF + CaF -→ CaF 2 + Ca chemical reaction has been constructed by Sardar and co-workers 40 .
Until very recently a proper quantum description of the title reaction was not feasible. Today, by employing unprecedented computational resources, it remains a very hard numerical task due to several reasons, namely: (i) the system lacks symmetries that could otherwise be used to ease parts of the computational overload; (ii) it is a somewhat heavy system with small diatomic rotation constants (e.g. the CaF constant is about 177 times smaller than that of H 2 ) and, as we shall see below, possesses a relatively deep potential well at short range, all of which translates into a large amount of spatially delocalized internal states required to properly describe the collision; (iii) it is known to be a very anisotropic system characterized by strong couplings between collisional channels that would be negligible otherwise; and (iv), typical of atom-molecule collisions within the cold domain of kinetic energies, the radial solution of the Schrödinger equation is required to be propagated to unusually large atom-molecule separations, due in part to the extremely long de Broglie wavelength of the colliding partners. Therefore, within the limitations imposed by such aspects, we provide below a first investigation on the optimal parameters required to extract accurate scattering characteristics for these collisions, in a time-independent quan-tum reactive scattering formalism, and discuss the predicted features of the collisional cross sections as functions of the incident energy. To this end, the adiabatically adjusting principal axis hyperspherical (APH) quantum reactive scattering suite of programs (hereafter referred to as APH3D), that has been used to describe a diverse array of reactive collisional problems in our group [41][42][43][44][45][46][47][48][49][50][51] , is also utilized below. While formalisms based on the solution of the time-dependent Schrödinger equation are computationally more efficient they are slow to converge at low collision energies 52,53 . Methods based on statistical quantum approaches [54][55][56] have also been applied to complex-forming ultracold chemical reactions but their accuracy for state-to-state transitions is not fully established.
The paper is organized as follows: Section 2 provides a brief description of the theoretical approach with details of the electronic structure calculations presented in subsection 2.1. A brief outline of the quantum scattering formalism using the APH3D code is presented in subsection 2.2. Section 3 presents the results and section 4 provides a summary of our findings.
The Li 2 S + CaF X 2 Σ + reactants asymptotically correlate with the triatomic electronic states 1 A and 3 A . In what follows we describe the computation of the ground electronic state, X 1 A , of the LiCaF complex using the internally contracted multireference configuration-interaction method with the Davidson correction (MRCI+Q) [57][58][59] , as implemented in the MOLPRO package 60 . The augmented correlation-consistent polarized valence quadruple-zeta basis set (aug-cc-pVQZ) of Dunning [61][62][63] was used for the Li and F atoms, whereas the cc-pwCVQZ-PP basis, in which the core electrons are described with a pseudopotential, was used for the Ca atom 64 . Calculations with full valence active space utilizing a state-averaged (1 1 A , 1 3 A and 1 1 A ) complete active space (10 active electrons in 9 active orbitals) selfconsistent field wavefunction (SA-CASSCF) 65,66 were performed. The active space included the 2s, 2s2p, 4s4p orbitals from Li, F, and Ca atoms, whereas the 1s orbitals for Li and F, along with the 3s3p orbitals of Ca were closed in the CASSCF calculations and further cored in the MRCI calculations.
A total of about 11000 geometries below 4.5 eV relative to the global minimum were selected and fitted using a many-body expansion method 67 V abc (r ab , r ac ,
ab (r ab ) +V (3) abc (r ab , r ac , r bc ) , (1) in which r xy is the internuclear distance between x and y (= a, b, or c); V
ab and V
(3) abc are the one-, two-, and three-body terms, respectively. The one-body terms in Eq. ( 1) are set to zero. The two-body terms correspond to the diatomic potential energy curves (PECs). The three-body energy becomes zero at all the dissociation limits.
The two-body terms, V
CaF and V
where D e is the dissociation energy, r e is the equilibrium distance of the diatoms, and α is a parameter. In the CaF case, D e = 5.45 eV, α = 0.51 a.u.
where ρ xy = r xy e -β xy r xy . The linear parameters, d jkl , can be obtained by the linear least squared method and the nonlinear parameters, β xy , are set to 0.5 a.u. -1 . Moreover, the constraints j + k + l = j = k = l and j + k + l ≤ M are employed to ensure the three-body term V
(3)
abc is going to zero at all dissociation limits. In this work, the value of M = 8 is used, which leads to a total of 140 d jkl linear coefficients. The root mean squared deviation (RMSE) of the three-body short-range fit is 22.7 meV.
The ab initio calculation yielded an exothermicity of -0.37 eV (-2984.2 cm -1 ) for the Li 2 S + CaF X 2 Σ + reaction, which is 0.13 eV (1048.5 cm -1 ) higher than the experimental value of -0.5 eV (-4033 cm -1 ). This error is corrected in the two-body terms which are adjusted to reproduce the experimental exothermicity.
The long-range interaction potential, V LR , in each arrangement is fitted with the following expression:
where
bc and R is the distance between the Li (Ca) atom and the center of mass of the CaF (LiF) molecule. The parameters l and n range from -3 to 3 and 4 to 7, respectively. For n = 4 and m = 1, B 1 4 (θ ) = cos θ ; for n = 5 and m = 1,
The errors in the long-range potential fitting for the Li + CaF and Ca + LiF arrangements are 2.55 and 1.86 cm -1 , respectively. In addition, the long-range and shortrange potentials are connected smoothly with a switch function. Specifically, we have
where the arrangement-dependent switching function, s abc , is defined as
When bc = CaF, γ s = 1 a.u. LiF internuclear distance (a.u.)
CaF internuclear distance (a.u.) Potential energy (cm -1 ) at θ = 180°F
ig. 1 Global PES, in cm -1 , as a function of the internuclear distances of LiF and CaF, in a.u., at the fixed angle θ = 180 • centered at the F atom. Isolines varying every 2000 cm -1 from -9000 cm -1 (innermost) to 11000 cm -1 (outermost).
LiF internuclear distance (a.u.) CaF internuclear distance (a.u.) Figure (1) depicts the contour plot of the PES, produced by the fitting procedure described above, as a function of the respective internuclear distances of LiF and CaF at the fixed bond-bond angle of θ = 180 • centered at the F atom. The red-shaded regions are associated with the lowest (attractive) values of the potential whereas the blue areas correspond to higher energies. For the purpose of the scattering calculations presented below, in what follows, the zero-energy of the PES is shifted to correspond to the energy of LiF at the equilibrium position of r LiF = 3.0 a.u., yellow area in Fig. (1). The light-green region is for the CaF potential well with equilibrium position of r CaF = 3.8 a.u., about 4040 cm -1 above zero, excluding zero point energy (ZPE). The LiCaF potential well, where the three atoms are in close proximity, reaches a minimum value (about -10931.7 cm -1 ) at slightly displaced diatomic distances, namely r LiF = 3.18 and r CaF = 4.06 a.u., in a near T-shape geometry at θ = 104.5 • -not shown but very similar to those contours of Fig. (1). In addition, Fig. ( 2) illustrate the strong anisotropic character of the system in which the PES is mostly repulsive for a collinear approach of the Li atom towards Ca (θ = 0 • ) in contrast with the attractive character whenever approaching on the F side as shown in Fig. (1).
As mentioned above, the APH3D code is utilized to model the title reaction using the PES described in the previous section. A somewhat detailed description of the numerical aspects and convergence criteria is provided in the next section whereas, for the sake of completeness, a brief overview on the implementation of APH3D is given here. The description provided below is along the lines of that given in our recent work on the H + D 2 chemical reaction 69 , however, an in-depth discussion of the hyperspherical coupled-channel equations, as implemented on APH3D, has been given by Kendrick and co-workers in many occasions [41][42][43][44][45][46][47][48][49][50][51] .
The LiCaF Hamiltonian is written in APH coordinates of Pack and Parker 41 . The hyperradius, ρ, describing the radial atomdiatom relative motion is partitioned into an inner region, using Smith-Johnson hyperspherical coordinates 42,70 , where collisioninduced re-arrangement is more likely to occur. In the outer region, where the different atom-diatom arrangement channels are largely decoupled, Delves hyperspherical coordinates [71][72][73] are employed. The six-dimensional three-body problem is reduced to a set of coupled equations along the scattering coordinate, ρ, with ρ discretized in a grid of N sectors. The eigenvalues associated to the remaining five internal degrees of freedom are used as the effective set of coupled potentials driving the relative motion along ρ.
The 5D eigenvalue problem is solved in the APH region by means of an implicitly restarted Lanczos method 74,75 , whereas the corresponding eigenvalues within the Delves region are evaluated using a 1D Numerov propagator 76 . Once a sufficiently large set of coupled potentials are evaluated in both regions for all sectors, as well as all sector-to-sector overlap matrices, the resulting set of radial coupled equations is solved using Johnson's log-derivative method 77 , first from ρ min to ρ match . At ρ match the numerical solutions from the outermost sector of the APH region are projected onto solutions at the innermost sector of the Delves region. The propagation is continued from ρ match to ρ max , a sufficiently large value of ρ where the interaction potential is negligible. At ρ max all channels (from all arrangements) are numerically decoupled, scattering boundary conditions are applied, the log-derivative solutions are projected onto solutions associated with each asymptotic diatomic state, written in ordinary Jacobi coordinates, yielding a scattering matrix 41 . The procedure is repeated independently for each value of the total angular momentum quantum number J and its parities, good quantum numbers in the absence of external fields. However, as explained below, for the present work, only the (even) J = 0 case is addressed. Moreover, the basis sets for both APH and Delves regions are independent of collision energy and, therefore, evaluated only once.
Due to the fact that low-lying collisional channels, with relatively higher kinetic energies, are associated with highly oscillatory components of the scattering wavefunction, particular attention is given below to the number of sectors, N, and the grid step size used, ∆ρ aph and ∆ρ delves . In addition, combined with the usual outward sector-to-sector integration of the Schrödinger equation, an intra-sector subdivision of the grid is employed in both APH and Delves regions. Thus, for the n th sector, of length ∆ρ, defined within the ρ n left and ρ n right boundaries, with
where k max is the maximum value of the wave vector considered, i.e.
In Eq. ( 7) µ is the atom-diatom reduced mass and E max is a fixed parameter whose value is set as high as the asymptotic energy of the highest closed channel included in the set of coupled equations, such that all channels are well described. Therefore, in what follows, we shall also determine the optimal value of N steps such that the grid-dependent physical description of the problem remains unaltered, i.e. a proper description of the smallest periods of oscillation of the wavefunction is included.
Despite a proper time-independent quantum formalism that takes into consideration the doublet (or higher) spin multiplicity is available in the domain of inelastic collisions [78][79][80][81][82][83][84] , an implementation of the reactive counterpart of the problem is not. Therefore, we shall use the formalism for collisions between a 1 Σ + molecule and a structureless atom. Such assumptions have been proven valid in certain context for inelastic collisions 83 . In what follows we make a few considerations for the case-study at hand.
Due to the null projection of the electronic orbital angular momentum of CaF on its internuclear axis, Λ(Σ + ) = 0, and the absence of a nearby electronic 2 Π state, only electrostatic interactions are expected to play a significant role on the internal structure of the molecule. Thus, the CaF effective (angular) Hamiltonian (neglecting vibrational and Stark terms) could be approxi-mated as 85,86
where N is the diatomic rotational angular momentum, S is the electronic spin angular momentum with S z = 1 /2 being the spin component along a given z-axis parallel to the internuclear axis, I is the nuclear spin with an I z = 1 /2 component (due to the 19 F isotope), and B e is the diatomic rotation constant. The parameters γ, b, c and f are the strength coefficients for the electronicspin-rotation, isotropic and anisotropic electronic-spin-nuclearspin, and nuclear-spin-rotation couplings, respectively. For convenience, the strength coefficients for the (υ = 0, N = 0) manifold, as measured by Childs and co-workers 87 , are given in Table (1).
Table 1 Experimental values, in cm -1 , of the molecular parameters of CaF for the υ = 0, N = 0 rovibrational manifold 87 .
As expected, the diatomic rotation constant is much larger than the remaining coupling parameters (about 257 times larger than γ and c, 94 times larger than b) and, therefore, the dominant term. As a consequence, the N = 1 rotational structure, within the υ = 0 manifold, is predicted to lie about 2B e ≈ 0.69 cm -1 (or 900 mK, neglecting higher-order centrifugal distortion contributions) above the N = 0 structure. Using either a Hund's case (a) 78 or (b) 80 notation, both N and S are generally well accepted to be weakly coupled. As a consequence, collision-induced changes in either the magnitude or the direction of the electronic spin, S, are unlikely to happen. However, a collision may induce sudden changes in N and, due to the subsequent recoupling between N and S, changes between the resultant parallel (e parity) and antiparallel ( f parity) coupling schemes may occur 78 . Since we are not properly describing the diatomic rotational structure, we shall address collision energies well below the 900 mK threshold, such that the N = 1 rotational state will remain as a closed channel.
In regard to the N = 0 fine/hyperfine structure, for Λ = 0, I z = S z = 1 /2, the predicted sublevels of H CaF are associated to j = 1 /2 and F = 0 or 1 quantum numbers, where j = N + S (fine structure) and F = j + I (hyperfine structure) 87 . However, as the collisions treated below are explored in the absence of external fields and, given the equally small electronic-spin-nuclear-spin interactions, b and c, alongside the negligible nuclear-spin-rotation interaction, f , the multiplet structure of the entrance channel is not considered henceforth. Even if external fields were taken into account, a collision-induced Zeeman relaxation of the N = 0 rotational struc-ture of CaF is expected to vanish at first-order, being mostly a second-order process 88 and, therefore, it appears reasonable to neglect. However, it is worthwhile to note that, as neither the doublet spin multiplicity of CaF nor that of Li is included, a resultant magnetic dipole-dipole interaction is also disregarded. As such interactions are more prevalent in the ultracold regime of kinetic energies, our assumption of a pseudo 1 S + 1 Σ colliding system implies a lower limit on the range of collision energies that can be studied here without compromises. As we shall see below, 1 mK is the minimal energy treated in this work.
Thus, from now on, we drop the use of the typical Hund's case (b) labeling of N j quantum numbers for the diatomic rotational structure in favor of the j rotational level (an integer, j = 0, 1, . . .), as used in the literature of singlet molecules.
Another aspect that we shall not describe in this work is higher values of the total angular momentum, i.e. J > 0. The amount of computational resources required to perform a J = 0 calculation is already substantial and the inclusion of higher J values would increase it drastically, as we would be now required to handle both even and odd parities of each non-zero J case. Despite the relatively low collision energies intended, yet well above the swave regime, it is likely that a few J values are still required to secure convergence. As a consequence, the J = 0 calculation presented below may not be suited for a direct quantitative comparison with experimental results. However, it is worthwhile to stress that J = 0 calculations have been proven to provide an insightful and accurate qualitative description of collisional problems in the past besides providing also the foundation for the optimization of certain key numerical parameters. )
Fig. 3 Set of the lowest 2500 APH coupled potential curves (cm -1 ), in which every 100 curves are highlighted in blue, as functions of the hyperradius, ρ (a.u.). The bold red line tags the CaF(υ = 0, j = 0) diatomic rovibrational level as a pseudo-1 Σ + molecule.
Thus, in the context of a pseudo-1 Σ + CaF molecule colliding with a structureless Li atom for J = 0, we then proceed to determine the optimal parameters in order to describe the scattering wavefunction in both APH and Delves regions. Figure (3) shows the final set of 2500 coupled potentials used in the APH region, evaluated between ρ min = 4.5 a.u. and ρ match ≈ 19.9 a.u. with 76 sectors varying logarithmically with a step size of ∆ρ aph = 0.02 atomic units. For easy visualization every 100th level is shown in blue. The bold red line represents the (υ = 0, j = 0) entrance channel of CaF. Despite many high-lying blue curves being asymptotically correlated with internal states a few thousand wavenumber above the entrance channel, they are strongly interacting at short range, ρ ≈ 6-7 a.u., and their inclusion is needed to achieve converged results. However, the inclusion of more channels increases the time-complexity of solving the Schrödinger equation by a few factors of O N 3 channels , where N channels is the number of channels to be included in the basis set of the scattering wavefunction. Thus, considering the number of calculations required to probe other convergence aspects (as shown below) and energy-dependent calculations, increasing the number of channels would quickly become impractical. If higher collision energies than those addressed here are of interest, it is desirable to include more channels, mainly in the range between ρ min and ρ ≈ 13.5 a.u., then projecting the solutions into a smaller basis set, as the one used here, and resuming the propagation toward large distances. The choice of ∆ρ aph = 0.02 a.u. will remain fixed in the remainder of this work. This is mostly due the inherent higher computational overhead of optimizing it as evaluation of different sets of coupled potentials would be required. However, it is worth noting that similar values have been used for converged calculations of systems as heavy as the one treated here and for somewhat similar grid parameters, e.g. ∆ρ aph = 0.01 a.u. used for Rb + K 2 by Croft et al. 89 and ∆ρ aph = 0.012 a.u. used for Li + LiNa by Kendrick et al. 90 (logarithmic scales used in all cases).
Beyond ρ match , now within the Delves region, the additional concern of how the scattering characteristics may vary with respect to ρ max , due to the lower collision energies, should be addressed. Fortunately, by inspecting Fig. (3) again, the set of coupled potential curves is seen to present a somewhat parallel behavior with respect to one another, at ρ match , mostly due to the smaller couplings, as ρ → ρ max . This aspect in particular suggests that a much smaller number of channels may be included in the basis set of the Delves region. In the present work, 600 channels are utilized to solve the Schrödinger equation along the Delves part, of which 396 are closed channels. Possibly, the number of basis functions used in the Delves region exceeds the requirements to obtain well converged numerical results at J = 0 and can be utilized also to describe collisions at higher energies.
Asymptotically, at ρ max , the diatomic eigenstates used as basis set comprises up to υ max = 7 and j max = 68 (υ = 0) for LiF, from which (υ = 0, j = 55), (υ = 3, j = 33) and (υ = 4, j = 21) are the highest open rovibrational manifolds. Whereas for CaF, υ max = 3 and j max = 77 (υ = 0) are utilized with all but the entrance channel energetically closed. This smaller basis set, when compared with the one used in the APH region, allows us to explore the convergence criteria with respect to N steps (APH and Delves regions) and ρ max (Delves region only) in great detail as described below. a.u.)
Fig. 4 Upper panel: Elastic component of the cross section (a.u.) for the Li + CaF collisions at 1 mK as a function of ρ max (a.u.) for N steps = 10 (black curve), 50 (red curve), 100 (green curve), 200 (blue curve), 400 (orange curve), 600 (magenta curve) and 800 (brown curve). Lower panel: Reactive component of the cross section (a.u.) for the Li + CaF(υ = 0, j = 0) -→ Ca + LiF(υ = 0, j = 0) chemical reaction at 1 mK as a function of ρ max (a.u.) using the same color code as in the upper panel. Circles are raw CC calculations, curves are Akima splines to enhance visualization.
We have used a reference collision energy of 1 mK relative to the (υ = 0, j = 0) entrance channel of CaF to compute both elastic and reactive cross sections for sets of N steps = 10, 50, 100, 200, 400, 600 and 800 combined with ρ max ≈ 85, 125, 155, 185, 305, 405, 805 and 1250 atomic units. We note that, all parameters except N steps are held fixed in the APH region, whereas in the Delves region, each N steps choice is combined with an increasing number of sectors varying linearly at fixed steps of ∆ρ delves = 0.2 atomic units. The result is presented in Fig. (4), where each curve corresponds to a given value of N steps and the cross section is plotted as a function of ρ max . An inspection of the upper panel (elastic component) suggests a somewhat strong dependence on both parameters, as expected, and the cross section converges from below to its optimal value between ρ max ≈ 250 and 350 a.u. with values of N steps > 600 yielding comparable results. The reactive cross sections for the LiF(υ = 0, j = 0) exit channel are presented in the lower panel of Fig. ( 4) where a similar convergence pattern (now from above) is evident, except that the set of calculations with N steps = 10 and 50 are completely unable to describe the reaction. Thus, N steps > 200 and, ideally, 600 is recommended.
For the sake of simplicity, the convergence behavior of cross sections for other choices of υ and j of LiF are not shown but they possess virtually identical patterns as those observed in the lower panel of Fig. (4). Instead, in Fig. (5), the reactive cross section for all open υ (panels) and j (abscissa) exit channels at the fixed values of ρ max = 145, 305, 405 a.u. and N steps = 600 is presented. A qualitative description of Fig. ( 5) is given in the next section whereas, for now, suffice to observe that each independent calculation (blue, red and brown) captures virtually identical branching ratios over both υ and j implying that the calculations are numerically stable in order to infer the actual physical aspects of these collisions.
First, we address how well the PES can reproduce the asymptotic PECs for the CaF(X 2 Σ + ) and LiF(X 1 Σ + ) subsystems at large atomdiatom separations. The bottom of each PEC is presented in Fig. (6), in which case the global dissociation limit corresponding to Li( 2 S) + Ca( 1 S) + F( 2 P), ≈ 47986 cm -1 , is not shown. Due to our earlier choice of using the LiF energy at the equilibrium position as the zero-energy in the scattering calculations, the 47986 cm -1 limit also corresponds to the relative dissociation energy of LiF, D e . For comparison purposes a list of a few selected values of equilibrium positions and dissociation energies, for both LiF and CaF electronic ground states, are collected in Table (2).
In the particular case of CaF, electronic structure data is somewhat scarce and/or dated. Yet, by inspection of Table (2), we do observe a reasonably good agreement between our calculations and those from literature, in particular, the recent results of Sardar and co-workers 40 . As expected, the dissociation energies appear to vary more broadly, within ≈ 2000 cm -1 among the various studies, with our result well within that range. Despite that the error appears to be relatively small if the actual total depth of the potential is taken into consideration (≈ 47986 cm -1 for LiF), the data collected in Table (2) seem to suggest that it is relatively harder to properly reproduce the LiF well depth than that of CaF. As investigated in great detail by Varandas 91 , the LiF electronic ground state is not trivial, manifesting a predominant ionic character at the equilibrium position and avoid crossing the 2 1 Σ + excited state (whose nature is essentially covalent) at relatively short ranges (r LiF ≈ 14 atomic units). Moreover, it is asymptotically correlated with an additional 1 Π state.
Using the PECs shown in Fig. (6), an energy splitting of about 0.65 cm -1 between the first two rotational levels of CaF is predicted, in the υ = 0 vibrational manifold, which suggests an effective diatomic rotational constant of B e = 0.65/2 = 0.325 cm -1 and, thus, is within 0.02 cm -1 from the value measured by Childs et al. -see Table (1). Similarly, an effective diatomic rotational constant of B e = 1.29 cm -1 is predicted for LiF, which agrees reasonably well with the measured value of 1.3452576 cm -1 92 , but overall these evidence seems to suggest that the shape of the PECs (and their bound states) is equally satisfactory. The energy levels for the vibrational states utilized in the scattering calculations are tagged with horizontal lines in Fig. (6) and, for the sake of clarity, only j = 0 cases are displayed, except for j = 37 (υ = 1) and j = 38 (υ = 5, the highest basis function taken into account). By including the respective ZPEs of each molecule a total exothermicity of about 3880.3 cm -1 is expected and, therefore, it is a few hundred wavenumbers below the earlier prediction of 4440 cm -1 by Kosicki and co-workers 36 .
Fig. 6 Diatomic potential energy curves for CaF(X 2 Σ + ) (black dashed curve) and LiF(X 1 Σ + ) (black solid curve), both in cm -1 , as functions of r with the respective third atoms (Li and Ca) at a distance of 100 a.u. (θ = 0 • ). Horizontal bold (solid and dashed) lines tag the diatomic vibrational levels ( j = 0) for the respective quantum number displayed. Horizontal bold dotted lines tag the vibrational levels corresponding to rotational states j = 37 (υ = 1) and j = 38 (υ = 5).
In the remainder of the paper we describe the Li + CaF -→ Ca + LiF chemical reaction with those parameters described in the previous section. We will perform a scan on collision energy from 1 mK to 200 mK for the (υ = 0, j = 0) entrance channel of CaF, in a grid of 128 points varying linearly, using N steps = 600 and ρ max = 305 a.u.; the result of which is presented in Fig. (7). The choice of the ground rovibrational state of CaF as the entrance channel for these collisions, at sufficiently small collision energies, rule out the occurrence of inelastic processes, such that the only non-elastic pathway is the chemical reaction. In addition, as the PES described above does not take into account the nearby 3 A electronic state (degenerate asymptotically), the influence of singlet-triplet nonadiabatic transitions and/or spin-exchange effects on the reaction presented below, if any, is disregarded. As no actual comparison with a measurement and/or other calculations is possible for now, those results presented below are not scaled by the typical 1 /4 statistical weight factor of singlet entrance channels with respect to their triplet counterparts. That implies a hypothetical scenario in which 100% of the colliding partners are prepared in the electronic ground state of the complex. In an actual experimental scenario, with no control of the initial spin, it is expected that up to 75% of the collisions would undergo elastic and inelastic processess along the triplet PES whereas 25% would Table 2 Equilibrium positions (r LiF and r CaF , in a.u.) and dissociation limits (D e , in cm -1 ) along with the method utilized (see the respective references for details). The original values, in Å and eV, are given between parentheses with the factors 1 a.u. = 1 a 0 = 0.52917721092 Å and 1 a.u. = 27.211385 eV = 219474.63137054 cm -1 applied. a The original D e value is converted from kcal/mol (135.8) and the geometry optimization is made at the MP2(full)/6311+G* level. b Within ± 0.3 eV (± 2419.7 cm -1 ).
Ref undergo the reactive process described here. In the discussion presented below, the absolute value of the cross section is less relevant and we shall address aspects of the relative quantities such as (branching) ratios.
In the upper panel of Fig. (7), the energy-dependence of the cross sections, summed over j , is presented for each open manifold associated to υ = 0-4 of the LiF product (cyan, red, green, blue and orange curves), whereas the total, summed over j and υ , is denoted by the solid black curve. For comparison purposes the elastic component of the cross section is also shown as the black dashed curve. As seen in Fig. (7) the reactive cross sections present a somewhat flat behavior whereas the elastic component is suppressed in the vicinity of 100 mK to 200 mK mostly due to the presence of a resonant feature. However, it is worth noting that it may be premature to consider this resonant structure as an actual observable feature due in parts to the fact that our calculation only represents the J = 0 case. The incoherent summation of contributions associated to higher J values may (and are likely to) wash out these features observed in Fig. (7). Thus, the question whether the resonant structure predicted here shall survive the addition of higher J values will remain open for further theoretical explorations. Likewise, the characterization of the resonance, in terms of angular momentum partial waves, width and lifetime is outside the scope of this work. However, as the entrance channel is associated to j = 0 and J = 0, and therefore only incoming = 0 partial wave contributes (J = j + ), there is no centrifugal term associated to the entrance channel potential curve, whose behavior is of an ordinary attractive potential. As a consequence, it is possible that the resonance-like feature shown in Fig. (7) is associated to a triatomic bound state belonging to another channel, i.e. a Feshbach resonance. This hypothesis is reinforced by the somewhat high density of states that may exist in the vicinity of the entrance channel -see the red bold line in Fig. (3).
The hierarchy of the j -summed cross sections for a given υ level may be understood from the following considerations. The υ = 4 manifold of LiF possess only 22 rotational levels that are open with respect to the (υ = 0, j = 0) entrance channel of CaF. This fact is illustrated, at 1 mK, in panel (a) of Fig. (5). As a result, the summation over j yields the smallest cross sections overall -see the magnitude of the orange curve in Fig. (7). Similarly, the υ = 3 and 2 cases possess the second and third smallest amount of open rotational states (34 and 42), and thus, provides the second and third smallest total cross section, i.e. blue and green curves of Fig. (7). Overall the rotational levels of LiF belonging to the υ = 4-2 cases are predicted to be poorly populated by the collision in the energy range described here. In contrast, the production of LiF in the υ = 0 and 1 manifolds, the largest in terms of open rotational states, are the chemical events with higher likelihood to occur, mostly populating j = 0-20, with a smaller but substantial probability of populating also highly rotationally excited states. Elastic-to-total ratio Fig. 7 Upper panel: Reactive cross section, in Å2 , for the Li + CaF(υ = 0, j = 0) -→ Ca + LiF(υ ) chemical reaction, summed over j , as a function of the collision energy, in mK; where, J = 0, N steps = 600, ρ max = 305 a.u., υ = 0 (cyan curve), υ = 1 (red curve), υ = 2 (green curve), υ = 3 (blue curve), υ = 4 (orange curve), the total summed over υ (black curve) and elastic component (dashed black curve). Lower panel: The elastic-to-total reactive cross section ratio as a function of the collision energy.
Another noteworthy aspect, as the collision energy increases, is the somewhat strong suppression of the elastic cross section that reaches a minimum value in the resonant region, at about 182 mK, as shown in the upper panel of Fig. (7). As a consequence, in the range of energies studied here, the Li + CaF collision may become predominantly reactive at collision energies in the vicinity of 200 mK. This fact is illustrated by the elastic-to-reactive ratio of the cross section presented in the lower panel of Fig. (7), where there are up to 1700 collisions for every chemical event at about 1 mK, remaining somewhat constant for about 100 mK, and quickly dropping to a minimum of 1:1 (or smaller) in the vicinity of 182 mK. For the sake of comparison, the elastic-to-inelastic ratio for a collision-induced Zeeman relaxation of CaF by collisions with He atoms, at much higher temperatures (2 K), has been measured by Maussang and co-workers to be about 10 4 13 . Likewise, in the cases of spin-polarized Li + CaH and Mg + CaH inelastic collisions, investigated by Tscherbul et al. 99 , for which chemical reaction is energetically forbidden, the elastic-to-inelastic ratio is predicted to be about 10 5 at 1 mK, i.e. nearly 60 times larger than the case considered here. In addition, CaH is known to be more amenable for magnetic traps, mostly due to its higher rotational constant compared to CaF 24 . This raises concerns on the prospects of sympathetic cooling of CaF by means of cold colli-sions with Li atoms above 100 mK due to potential trap losses induced by the formation of LiF.
From an experimental point of view, our choice of using CaF in its lowest internal state as the entrance channel could be realized by producing the molecule with either a MOT or a Stark decelerator or a combination of these with a microwave trap. A modern MOT implementation is capable of producing molecules for collisions at energies as low as the Doppler limit whereas a Stark deceleration method is likely to produce molecules with a temperature of a few dozens of mK, being therefore more problematic for the case at hand. In either case, however, typical procedures, such as compressing the molecular cloud in order to improve its overlap with the buffer gas coolant may eventually raise the temperature by a few extra dozens of mK and, therefore, also trigger losses due to LiF formation. Those experimental implementation already reaching sub-Doppler temperatures should not be concerned by losses due to chemical reaction but sub-Doppler heating effects, as those demonstrated by Devlin and Tarbutt 100 , may occur for certain kind of MOTs.
A detailed single-arrangement Lennard-Jones-based, and thus disregarding reactivity, simulation of the thermalization (no inelasticity either) of CaF in the presence of Li and Rb cold buffer gases has been carried out by Lim et al. 37 . Their numerical experiment assumed a practical experimental scenario similar to that given above and found a somewhat strong dependency between the cooling rate and the s-wave scattering amplitude for those scenarios with ultracold Li atoms, when compared to Rb, mostly due to the relatively small reduced masses for the Li + CaF combination. In addition, they have also predicted a slowdown of the collision process, and thus the cooling rate, for the Li + CaF case, due to a minimum in the cross sections in the range of 1-10 mK, similar to that found in the present work at higher energies, about 100-200 mK. In contrast, a similar minimum was predicted by Lim et al. within the µK range of collision energies when ultracold Rb atoms were used. As a consequence, the cooling rate when using Li was found to be an order of magnitude slower than that for the Rb case 37 .
Alternatively, as also pointed out by Lim et al. 37 , the use of a light colliding partner such as Li for sympathetic cooling when CaF is produced in an excited rotational state may be favorable due to potentially higher centrifugal barriers that could, as a consequence, suppress losses due to collision-induced inelastic processes. This scenario is yet to be investigated but it is now possible using the PES we have presented here. Moreover, it may be worthwhile to explore collisions driven by the 3 A electronic state of the LiCaF complex and explore the possibility to control the reaction by means of external magnetic fields as in the case of the Li (Mg) + CaH systems 20,99,101 . Chemical reactions are likely to be suppressed in spin-polarized collisions of 2 S + 2 Σ systems on the 3 A PES, due to the less attractive character of the 3 A electronic state and overall endothermicity. Moreover, there is evidence suggesting that the spin-orbit-induced triplet-to-singlet transition, that could trigger the formation of LiF in the 1 A PES, as shown here, may be either small or negligible 101 . An overview of the 3 A electronic state of the LiCaF complex has been given by Frye and co-workers 102 .
Overall, our results show that cold collisions of Li and CaF favor elastic scattering in the 1-100 mK regime but a sharp decrease in the elastic cross section in the vicinity of 200 mK, possibly due to a Feshbach resonance, makes the elastic/reactive cross section ratio < 1, limiting the efficacy of sympathetic cooling of CaF by collisions with cold Li atoms. Despite the high density of asymptotic diatomic states and bound triatomic states that are involved in the collisions -see Fig. (3) -, our calculations predict a somewhat low density of resonances. This is probably due to the downhill nature of the reaction and presumably the short lifetimes of the LiCaF complexes formed. This aspect, the effect of rotational and vibrational excitation of the CaF molecule, and a proper characterization of the resonance will be addressed in future work. Indeed, a recent quantum close-coupling study of Ca + BaCl + system has shown strong vibrational quenching rates for BaCl + that exceeds rotational quenching rates for low-lying rotational levels 103 .
In this work we have applied state-of-the-art quantum chemistry and quantum reactive scattering to study both the interaction and dynamics of Li 2 S + CaF X 2 Σ + , in the context of cold collisions. To this end we have produced a global potential energy surface for the ground electronic state of the LiCaF system, X 1 A , capable of describing both atom-diatom arrangements, Li + CaF and Ca + LiF. The electronic structure calculations were carried out using an internally contracted multi-reference configuration-interaction method with a state-averaged (1 1 A , 1 3 A and 1 1 A ) complete active space (10 active electrons in 9 active orbitals) self-consistent field electronic wavefunction. A total of about 11000 geometries were evaluated and used to produce the final potential energy surface fit with a many-body expansion method augmented with ab initio parameterized long-range potentials.
Scattering calculations for the Li + CaF(υ = 0, j = 0) entrance channel were performed between 1 and 200 mK of collision energy. At 1 mK the collision-induced formation of rovibrationally excited LiF(υ = 0-1, j ≈ 0-20) molecules is predicted to be the most likely collisional outcome, with a total energy release that could reach up to 3880 K. In the vicinity of 100-200 mK a quantum resonance, likely to be a Feshbach resonance, appears to strongly suppress the elastic component. The reactive cross sections, however, remain largely unaffected in this regime, presumably due to its small magnitude compared to its elastic counterpart. The overall effect is that the elastic-to-reactive cross section ratio falls well below the lower limit of one hundred for collision energies above 100 mK suggesting a somewhat poor cooling rate for sympathetic cooling of CaF by Li and a strong trap loss due to the formation of LiF. At the resonance energy of 182 mK nearly every collision is predicted to be reactive (1:1 ratio or smaller).
It is worthwhile to emphasize, however, that the calculations presented here are not yet accurate for direct comparisons with future experimental observations as most likely a single PES, and the single partial wave (J = 0) used in the scattering calculations are insufficient. However, we believe it will serve as a benchmark for further theoretical works as we provided a detailed description of the potential energy surface and of those numerical aspects re-quired to obtain reasonably well converged scattering characteristics, a substantial improvement upon previous studies that were limited to model potentials and elastic/inelastic collisions, oftentimes considering only the equilibrium geometry of the triatomic complex, and equally limited dynamical models.
1-14 Journal Name, [year], [vol.],
Journal Name, [year], [vol.],