A first-principles based quantum dynamics study of the Li + LiNa($v$ = 0, *j* = 0) → Li_{2}($v\u2032$, *j*′) + Na reaction is reported for collision energies spanning the ultracold (1 nK) to cold (1 K) regimes. A full-dimensional *ab initio* potential energy surface for the ground electronic state of Li_{2}Na is utilized that includes an accurate treatment of the long-range interactions. The Li + LiNa reaction is barrierless and exoergic and exhibits a deep attractive potential well that supports complex formation. Thus, significant reactivity occurs even for collision temperatures approaching absolute zero. The reactive scattering calculations are based on a numerically exact time-independent quantum dynamics methodology in hyperspherical coordinates. Total and rotationally resolved rate coefficients are reported at 56 collision energies and include all contributing partial waves. Several shape resonances are observed in many of the rotationally resolved rate coefficients and a small resonance feature is also reported in the total rate coefficient near 50 mK. Of particular interest, the angular distributions or differential cross sections are reported as a function of both the collision energy and scattering angle. Unique quantum fingerprints (bumps, channels, and ripples) are observed in the angular distributions for each product rotational state due to quantum interference and shape resonance contributions. The Li + LiNa reaction is under active experimental investigation so that these intriguing features could be verified experimentally when sufficient product state resolution becomes feasible for collision energies below 1 K.

## I. INTRODUCTION

The experimental techniques for cooling and trapping molecules at cold (*T* < 1 K) and ultracold (*T* < 1 mK) temperatures continue to advance at a remarkable pace.^{1–16} The first experimental demonstration of a controlled ultracold chemical reaction was reported for the benchmark KRb + KRb system several years ago.^{3} Recent experimental techniques now routinely prepare cold molecules in a specific initial rovibrational state and can now also probe and detect specific intermediate complexes and product states.^{14,15} Product angular distributions of cold non-reactive (inelastic) collisions have also been recently measured for the HD + H_{2} and HD + He systems.^{5,12,13} Reactive systems involving heteronuclear alkali metal dimers such as KRb, NaK, NaRb, and LiNa are under active experimental investigation, and we anticipate that the measurement of fully state-resolved initial and product state rate coefficients and angular cross sections will soon become feasible.^{3,6,8,9,14,15} Thus, accurate first-principles based theoretical state-to-state quantum mechanical calculations of these systems are needed. For these systems, the reaction path is barrierless and exoergic so that significant reactivity can occur even for collision energies approaching absolute zero. At these extremely low temperatures, quantum mechanical (e.g., interference) effects dominate and the tiny relative kinetic energy allows for the possibility to steer and control these reactions via the application of external fields,^{17,18} the preparation of a particular initial quantum state,^{3,4,8,9,16} or the choice of the relative orientation of the colliding partners (i.e., stereodynamics).^{13,19–22}

Accurate ultracold quantum dynamics calculations are very challenging for a number of reasons: (1) despite the small relative collision energy, the density of states in the interaction region can be quite large due to the relatively heavy masses of the nuclei and the deep attractive potential wells, (2) the tiny collision energies and long-range van der Waals interactions require a much larger computational grid, (3) an accurate and stable numerical method capable of treating a large dynamic range is required, and (4) an accurate full-dimensional *ab initio* potential energy surface (PES) is required that includes a smooth and accurate treatment of the long-range region as well as accurate asymptotic diatomic potentials. As the collision energy approaches the ultracold regime, by definition, only a single partial wave contributes to the scattering process. This property reduces the size of the calculations significantly. However, at the higher cold collision energies approaching 1 K, several partial waves contribute and these must also be included in the theoretical calculations. For these reasons, the majority of first-principles based quantum mechanical calculations have been limited to a single partial wave, which is sufficient for the ultracold regime: these include O + OH,^{23,24} H + H_{2} (also its isotopic variants),^{25–32} Li + LiYb,^{33} K + KRb,^{34,35} and Li + LiNa.^{36} For the H + H_{2} and O + OH reaction systems, the calculations have been extended into the cold energy regime by including several partial waves. However, until now, the treatment of the heavier alkali metal (non-hydrogen) systems has been limited to the ultracold regime and a single partial wave.

In this work, we report accurate first-principles based quantum reactive scattering calculations for the Li + LiNa($v$ = 0, *j* = 0) → Li_{2}($v\u2032$, *j*′) + Na reaction. The calculations are based on a time-independent coupled-channel (CC) formulation using hyperspherical coordinates. A recently developed state-of-the-art *ab initio* PES for Li_{2}Na is used, which includes an accurate treatment of the long-range interactions and diatomic potentials.^{36} This PES was recently used in a quantum mechanical calculation for the ultracold Li + LiNa reaction but was limited to one partial wave (i.e., *l* = 0).^{36} In the present study, the collision energies span the entire range from *E*_{c} = 1 nK to 1 K and include all the required partial waves (i.e., up to *l* = 10 at *T* = 1 K). Both the total and rotationally resolved rate coefficients are reported as a function of the collision energy. Particular focus of the present work is on the angular distributions or differential cross sections (DCSs), which are also reported as a function of the collision energy and scattering angle. Intriguing quantum interference effects are seen in the DCSs as a function of both energy and scattering angle that give rise to prominent features: channels, sharp dips, and curved ridges. Furthermore, some of these features appear at small collision energies *T* ≈ 1 mK where only two or three partial waves contribute. These features or “quantum fingerprints” are unique to each Li_{2}($v\u2032$, *j*′) product state. As recently shown in a series of theoretical studies,^{23–29} the unique properties of ultracold and cold chemical reactions can lead to a dramatic enhancement or suppression of the reactivity due to quantum interference effects. These quantum interference effects could be exploited by experimentalists to control the reaction outcome. We show here for the first time how these quantum interference effects manifest in the angular distributions of a reactive alkali system under active experimental investigation.

Section II discusses the computational methodology, numerical parameters, and the *ab initio* PES used to perform the calculations. The quantum reactive scattering results are then presented in Sec. III, which include the total and rotationally resolved rate coefficients and angular distributions. Section IV concludes the paper with summarizing discussion and future outlook.

## II. METHODOLOGY

The quantum reactive scattering calculations are performed using a time-independent coupled-channel approach in hyperspherical coordinates and are reviewed in Sec. II A. The methodology makes no dynamical approximations and is numerically exact for a given PES. A first-principles based *ab initio* Born–Oppenheimer PES for the ground electronic state of Li_{2}Na is used, which includes an accurate treatment of the long-range interactions and diatomic potentials. This PES is discussed in Sec. II B, and the numerical convergence and various parameters used in the calculations are discussed in Sec. II C.

### A. Quantum reactive scattering in hyperspherical coordinates

In the theoretical description of molecules, the conventional Born–Oppenheimer approximation is typically used, where the Schrödinger equation for the nuclear motion is expressed in terms of a ground state electronic PES, which is a function of the internuclear distances. In the present work, we consider triatomic molecules for which there are three nuclei and three internuclear distances. The three-dimensional electronic PES is assumed to be calculated separately for each fixed internuclear distance of interest and is discussed in detail in Sec. II B. The Born–Oppenheimer equation for the nuclear motion (relative to the center of mass of the triatomic system) is given by

where the gradient operator ∇ spans the six nuclear coordinates denoted by $x=(x,x^)$ and x and $x^$ denote the three internuclear distances and three Euler angles, respectively. The three Euler angles specify the orientation of the molecular body frame relative to the space-fixed frame and are typically denoted by $x^=\alpha ,\beta ,\gamma $. The three-body reduced mass for the nuclei is *μ*, and *E* specifies the total energy of the molecule. The corresponding Born–Oppenheimer total molecular wave function is given by

where *ψ*(**x**) is the nuclear motion wave function of Eq. (1), *ϕ*(**r**) is the electronic ground state wave function, and $\psi \u2009N$ is the total nuclear spin wave function for the three nuclei. The electronic wave function and PES (*V*) are assumed to be previously computed and satisfy the electronic Schrödinger equation,

where *h* denotes the electronic Hamiltonian operator, which is a function of all the electronic coordinates (collectively denoted by **r**) and the three internuclear distances x.

In the present work, we seek numerically exact quantum mechanical solutions of Eq. (1) for a given total energy *E* and ground electronic state PES. We choose to work in hyperspherical coordinates that have several advantages: (1) a simultaneous and democratic treatment of all three arrangement channels, (2) a straightforward treatment of the identical nuclei exchange symmetry, and (3) a single radial coordinate that facilitates a numerical separation of the problem. In addition, two variants of hyperspherical coordinates are utilized: (a) the adiabatically adjusting principal axis hyperspherical (APH) coordinates of Pack and Parker^{37} and (b) the Delves channel-centered hyperspherical coordinates.^{38} The APH coordinates are optimal in the interaction region where all three arrangement channels overlap and are mixed, whereas the Delves hyperspherical coordinates are optimal outside the interaction region where the three arrangement channels are well separated and uncoupled. Both sets of coordinates share the same radial coordinate, which makes the transformation from the APH to Delves hyperspherical coordinates straightforward to implement.^{37}

The APH coordinates are denoted by *ρ*, $\theta \u0303$, and *ϕ*, where *ρ* is the radial coordinate and corresponds to a symmetric stretch motion, $\theta \u0303$ is the polar coordinate (related to the *θ* of Pack and Parker via $\theta \u0303=\pi \u22122\u2009\theta $) and corresponds to a bending motion, and *ϕ* is the azimuthal coordinate and corresponds to an internal pseudo-rotational motion. These three APH coordinates can be expressed in terms of the three internuclear distances.^{37} Collectively, these three plus the three Euler angles specify the six internal degrees of freedom of a triatomic molecule. For more detailed descriptions of these coordinates, see Refs. 37 and 39. Expressing Eq. (1) in terms of the APH coordinates in the interaction region gives

where the operator $\Lambda ^$ is the grand angular momentum operator. The explicit expression for this operator is given elsewhere and will not be duplicated here.^{40} It includes the $Jx2$, $Jy2$, $Jz2$, and the Coriolis coupling term $Jz\u2009\u2202\u2202\varphi $. The PES strongly couples the three hyperspherical coordinates. Nevertheless, we can perform a numerical separation of Eq. (4) by first solving the five-dimensional (5D) angular problem on a discretized grid of fixed *ρ* values. That is, we split Eq. (4) into two parts. First, we solve the 5D angular equation given by

where *ρ*_{ξ} denotes the fixed (discretized) value of *ρ* and the $\Phi tJMpq$ are the 5D hyperspherical (angular) functions. The solutions are labeled by the integer *t*, *p* denotes inversion parity, *q* denotes the particle exchange symmetry that is relevant for triatomic molecules with identical nuclei, and $w\u2261(\theta \u0303,\varphi ,\alpha ,\beta ,\gamma )$. In practice, we expand the angular solutions in terms of a primitive basis of Jacobi polynomials in $\theta \u0303$ and complex exponential functions in *ϕ*.^{40} A hybrid approached is used based on a Discrete Variable Representation (DVR) in $\theta \u0303$ and a Finite Basis Representation (FBR) in *ϕ*. The hybrid approach has several advantages, as discussed in detail in Ref. 40: (1) it allows for the truncation of the basis set via the Sequential Diagonalization Truncation (SDT) method, (2) it accurately treats the Eckart singularities associated with body-frame coordinates, and (3) it gives a real-symmetric Hamiltonian matrix even for non-zero total angular momentum *J*. Once the Hamiltonian matrix is constructed in the hybrid basis, it is diagonalized in parallel using the parallel version of ARPACK.^{41} The lowest *N*_{ch} solutions are computed and used as a basis for solving the one-dimensional (1D) radial equation in *ρ*. That is, we expand the wave function *ψ* in Eq. (4) in terms of the angular solutions of Eq. (5) at each fixed value of *ρ* = *ρ*_{ξ},

where *i* and *t* denote the coupled-channel indices (*i* and *t* = 1, 2, 3, …, *N*_{ch}). The radial coefficients $\zeta itJpq$ are computed numerically from the propagation of the resulting set of coupled-channel (CC) radial equations,

where the potential coupling matrix is defined as

The coupled-channel equations (7) and (8) are solved using Johnson’s log-derivative method^{42,43} and the optimal number of coupled channels (*N*_{ch}) is determined from convergence studies. This value and other parameters used in the calculations [such as the size of the primitive basis sets used to solve the angular equation (5)] will be discussed in Sec. II C.

The Delves hyperspherical equations are similar to the APH equations (4)(5)(6)(7)(8) and are explicitly given in Ref. 38 and will not be duplicated here. By construction, in the Delves region, the three arrangement channels are well separated and decoupled so that the solutions for each channel can be computed separately. Furthermore, the angular or rotational part of the Delves hyperspherical functions is analytic and can be expressed in terms of the well-known spherical harmonics. Thus, only a 1D numerical solution of the vibrational degree of freedom is required. An efficient Numerov method^{44} is used to compute the vibrational wave functions for each value of *j* (the diatomic rotational quantum number) and *l* (the orbital angular momentum quantum number). The rovibrational solutions are then used as the basis for solving the Delves version of the coupled channel equations (7) and (8).

To summarize, the quantum reactive scattering calculations are performed in a series of steps: (1) the APH surface (angular) functions are computed on a discretized grid from *ρ* = *ρ*_{i} to *ρ* = *ρ*_{match}, (2) the Delves rovibrational functions are computed for each channel on a discretized grid from *ρ* = *ρ*_{match} to *ρ* = *ρ*_{f}, (3) the overlap matrix between the APH and Delves functions is computed at *ρ* = *ρ*_{match}, (4) the APH log-derivative matrix is propagated from *ρ* = *ρ*_{i} to *ρ* = *ρ*_{match}, (5) the overlap matrix from step (3) is used to transform the APH log-derivative matrix into the Delves basis, (6) the Delves log-derivative matrix is then propagated from *ρ* = *ρ*_{match} to *ρ* = *ρ*_{f}, (7) the reactance **K** matrix is computed at *ρ* = *ρ*_{f} using the Delves log-derivative matrix, Delves basis, and asymptotic Jacobi basis, and, finally, (8) the full state-to-state **S** matrix is computed from the **K** matrix. The integral cross sections, rate coefficients, and DCSs can then be computed from the **S** matrix using standard expressions.^{37} The APH and Delves solutions in steps (1) and (2) are independent of the collision energy and only need to be computed once for each value of total angular momentum *J*, inversion parity *p* = ±, and exchange symmetry *q* = even/odd. The log-derivative propagation (steps 4–6) must be repeated for each collision energy and used in steps 7 and 8 to compute the scattering matrix as a function of collision energy. Many of the calculations (i.e., for each *J*, *p*, and *q*) are independent of each other and can be performed simultaneously if sufficient computational resources are available.

### B. *Ab initio* electronic potential energy surface for Li_{2}Na

A state-of-the-art *ab initio* based electronic PES for the Li_{2}Na molecule was recently reported in Ref. 36. These first-principles based theoretical calculations include the treatment of both the ground and first excited electronic doublet states (^{2}A′ and ^{2}B′) of the Li_{2}Na molecule. We give only a brief summary of the *ab initio* calculations here (see Ref. 36 for a more detailed description). The *ab initio* calculations were performed using the MOLPRO^{45} quantum chemistry package at over 40 000 internuclear geometries. Core potentials were used for the filled electron shells and include polarization and core–valence interactions. The valence electrons were treated using a multi-configurational self-consistent field (MCSCF) method to compute an initial set of configuration state functions (CSFs). The CSFs were then used with a large active space to perform multi-reference configuration interaction (MRCI) calculations for the two lowest adiabatic PESs of Li_{2}Na. The *ab initio* points were fit using a reproducing kernel Hilbert space technique. An accurate treatment of the long-range interactions was also included through the introduction of a three-body dispersion term and a smooth switching function that connects the short range three-body potential with the long-range three-body dispersion term. The dynamical effects of using a non-additive three-body long-range potential are discussed in Ref. 33. In addition, the MOLPRO computed *ab initio* potential curves for Li_{2} and LiNa were replaced with spectroscopically accurate diatomic potentials. The non-adiabatic derivative coupling matrix was also computed and used to obtain a two-state diabatic potential matrix. In the present work, we restrict our quantum dynamics treatment to the single adiabatic ground electronic ^{2}A′ state of Li_{2}Na.

Figure 1 plots the adiabatic ground electronic ^{2}A′ state PES of Li_{2}Na at several fixed values of the hyperradius *ρ*. All energy contours are relative to the bottom of the asymptotic potential of the Li_{2} + Na product and the dashed contour line at 2.23 × 10^{3} K denotes the total energy of the reactant LiNa(*v* = 0, *j* = 0). The PES is symmetric across the horizontal axis (*ϕ* = 0) due to the symmetry associated with the two identical Li nuclei. Thus, wave functions of even and odd exchange symmetry correspond to even or odd symmetry with respect to *ϕ* → −*ϕ*. Asymptotically, these wave functions correlate with the even or odd rotational levels of Li_{2}. Starting with large *ρ* = 25.0 a_{0} [panel (d)], we see that the three arrangement channels are isolated and well separated. As *ρ* decreases to *ρ* = 16.0 a_{0} [panel (c)], the LiNa and Li_{2} arrangement channel regions have spread out, and we see that the barrier between them is smallest near the equator, which corresponds to collinear geometries. Thus, the minimum energy reaction path lies along a collinear approach of Li with LiNa. For smaller *ρ* = 12.0 a_{0} [panel (b)], the arrangement channels have merged and the region surrounded by the dashed contour encompasses all three arrangement channels. In the interaction region *ρ* = 8.5 a_{0} [panel (a)], the three arrangement channels are now entirely mixed. The region lying within the dashed contour spans a large area surrounding the north pole of the hypersphere and includes two symmetric deep potential wells. The red dot in panel (a) denotes the location of the C_{2v} conical intersection (CI), which occurs between the ground ^{2}A′ and first excited ^{2}B′ electronic states. In the present work, we ignore the excited electronic state and its associated non-adiabatic and geometric phase effects (see Ref. 36 for a non-adiabatic quantum dynamics calculation that includes these effects). In summary, the Li + LiNa → Li_{2} + Na reaction proceeds along a barrierless collinear reaction path. In the interaction region, the reactants “fall” into a deep attractive potential well and form an intermediate complex. This complex subsequently decays into the various open Li_{2}(*v*′, *j*′) product rovibrational channels. The quantum mechanical calculation of the associated product state-resolved rate coefficients and angular distributions is presented in Sec. III.

### C. Numerical parameters

The accuracy of the APH surface function solutions computed by diagonalizing the Hamiltonian operator in Eq. (5) depends upon the size of the primitive basis sets in $\theta \u0303$ and *ϕ* specified by the integers *l*_{max} and *m*_{max}, respectively.^{40} The accuracy also depends on the energy cutoff (*E*_{cut}) used in the SDT procedure. A series of convergence studies are performed at several fixed values of *ρ* to determine a set of optimal values for these parameters. Typically, larger basis sizes are required as *ρ* increases due to the localization of the channels (see Fig. 1). Several full scattering calculations are then performed to also verify the convergence of the cross sections with respect to *l*_{max}, *m*_{max}, and *E*_{cut}. The convergence studies are performed for zero total angular momentum *J* = 0, and the same parameters are then used for all *J* > 0. The optimal values for *l*_{max} used in the calculations are 99, 119, 127, and 127, and for *m*_{max}, they are 220, 240, 260, and 280 for the four ranges of *ρ*: 6.0 ≤ *ρ* ≤ 11.024 a_{0}, 11.024 < *ρ* ≤ 19.544 a_{0}, 19.544 < *ρ* ≤ 23.374 a_{0}, and 23.374 < *ρ* ≤ 33.034 a_{0}, respectively. The basis size in $\theta \u0303$ is given by $n\theta \u0303=lmax+1$. For *ϕ*, the basis size is given by *n*_{ϕ} = 2 *m*_{max} + 1 and *n*_{ϕ} = 2 *m*_{max} + 2 for even and odd *J*, respectively. In the APH region, the hyperradius was discretized at 144 logarithmically spaced points between *ρ*_{i} = 6.0 and *ρ*_{f} = 33.034 a_{0} inclusive. The total dimension of the APH angular Hamiltonian matrix that is diagonalized at each *ρ* before applying the basis reduction via SDT is given by $n=n\theta \u0303\u2009n\varphi \u2009n\Omega $, where *n*_{Ω} denotes the number of Ω blocks for a given value of *J*. The value of Ω denotes the projection of *J* on the body-frame *z* axis and is even (odd) for even (odd) inversion parity *p*.^{40} We note that in this work, the initial reactant rotational state is *j* = 0 so that *J* = *j* + *l* = *l* (where *l* is the orbital angular momentum of Li about LiNa) and *J* + *p* is always even. Thus, only the *J*^{p} values given by *J*^{p} = 0^{+}, 1^{−}, 2^{+}, 3^{−}, … contribute to the cross sections. The values of Ω for both even and odd *J* span the range Ω = +*J*, *J* − 2, *J* − 4, *J* − 6, …, −*J* and *n*_{Ω} = *J* + 1. Thus, the total dimension of the angular matrix increases linearly with *J*. For *J* = 0, the total matrix dimensions in the four regions of *ρ* given above before applying the SDT are 44 100, 57 720, 66 688, and 71 808, respectively. After SDT, the maximum matrix size in each region is significantly reduced to 13 008, 9412, 8401, and 4814, respectively. All these matrix dimensions scale as *n*_{Ω} (i.e., linearly with *J*). Thus, for the maximum value of *J*^{p} = 10^{+} considered in this work, the largest matrix dimension was 13 008 × 11 = 143 088. Clearly, the SDT is critical to making the calculations feasible. The value of the SDT energy cutoff (*E*_{cut}) is also optimized at each *ρ* and is typically about 1.5 times larger than the highest lying eigenvalue of interest. For reference, the values of *E*_{cut} at the center of each of the four regions of *ρ* listed above are 1.09, 0.658, 0.752, and 0.746 eV, respectively (relative to the bottom of the asymptotic Li_{2} + Na potential). Fortunately, the large matrices are very sparse, and only, the diagonal Ω blocks and the off-diagonal coupling between their nearest neighbors (i.e., Ω ± 2) are dense. In addition, only a relatively small percentage of the lowest-lying eigensolutions are needed for the coupled-channel equations. Thus, these matrices can be diagonalized using an efficient sparse matrix diagonalization technique such as the Implicitly Restarted Lanzcos Method (IRLM), as implemented in the ARPACK library.^{41} The IRLM diagonalization scales as $nSDT2\xd7Nch$, where *n*_{SDT} is the dimension of the reduced SDT matrix and *N*_{ch} is the number of required eigensolutions (or coupled channels, see below). The matrix diagonalization at all 144 *ρ* values can be done independently. Furthermore, the calculations for all values of *J*^{p} and each exchange symmetry *q* = *even*/*odd* can also be performed separately. Thus, all these matrix diagonalizations can be performed simultaneously if sufficient computational resources are available. Once the eigensolutions have been computed and stored on disk, the overlap and potential coupling matrices are then computed at each *ρ* for use in the log-derivative coupled-channel propagation [see Eqs. (7) and (8) and Refs. 30, 40, and 46 for additional details].

The primary convergence parameter in the log-derivative propagation is the number of coupled-channels *N*_{ch} (i.e., the number of angular eigensolutions described in the preceding paragraph). The *J* = 0 scattering calculations are performed for several values of *N*_{ch} to determine the optimal value. For *J* > 0, the optimal value for *N*_{ch} scales approximately linearly with respect to *n*_{Ω}, that is, $Nch\u2248Nch0\xd7n\Omega $ (where $Nch0$ denotes the *J* = 0 value). The dense linear algebra required by the propagation involves matrix–matrix multiplication and matrix inversion at each step. Thus, the computational cost of the log-derivative propagation scales as $Nch3$. However, the calculations at each collision energy and for each value of *J*^{p} and exchange symmetry *q* = *even*/*odd* can be done separately. In addition, for large values of *N*_{ch} > 2500, the linear algebra can be done more efficiently using the parallel ScaLAPACK library.^{47} If sufficient resources are available, the parallel ScaLAPACK library enables the treatment of much larger values of *N*_{ch} and therefore larger *J*. At *ρ*_{match} = 33.034 a_{0}, the APH log-derivative matrix is transformed into the Delves basis using the overlap matrix computed between the APH angular functions and the Delves channel-centered rovibrational functions. The log-derivative propagation is then continued into the asymptotic region until *ρ* = *ρ*_{f}. The final value of *ρ* is also determined from convergence studies, and in this work, it was chosen to be *ρ*_{f} = 144.574 a_{0}. Several full scattering calculations are performed for *J* = 0 and various values of *ρ*_{f} (and also *ρ*_{match}) to verify convergence. Thus, the coupled-channel propagation in the Delves region spans 338 uniformly spaced values of *ρ* between *ρ* = 33.034 and *ρ* = 144.574 a_{0} inclusive. Table I lists the number of channels used in the log-derivative propagation in the APH and Delves regions as a function of *J*. From Table I, we see that the number of Delves channels is less than the number of APH channels. In the APH interaction region, the deep potential wells give rise to a larger number of open channels relative to the asymptotic Delves region. Another convergence parameter is the number of substeps used in each *ρ* sector. For a uniform grid, the sector at *ρ*_{χ} is defined as the region *ρ*_{χ} ± Δ*ρ*/2. Each sector is subdivided into smaller propagation steps Δ*ρ*/*n*_{step}, where *n*_{step} is an integer determined from convergence studies. For the APH and Delves regions, the optimal values of *n*_{step} are 100 and 200, respectively. The propagation across all of the APH and Delves sectors was performed at 56 logarithmically spaced collision energies spanning the range from 1 nK to 1 K (see the triangles in Fig. 2). The total scattering energy is given by *E* = *E*_{c} + *E*_{0}, where *E*_{c} denotes the relative collision energy of LiNa + Li and *E*_{0} denotes the rovibrational ground state energy of LiNa(*v* = 0, *j* = 0) relative to the bottom of the asymptotic Li_{2} + Na potential well (in this work, *E*_{0} = 0.192 002 eV or 2.228 × 10^{3} K). The largest collision energy considered here is *E*_{c} = 1 K, which gives *E* = 0.192 088 3 eV or 2.229 × 10^{3} K.

J^{p}
. | APH region . | Delves region . |
---|---|---|

0^{+} | 800 | 500 |

1^{−} | 1600 | 1000 |

2^{+} | 2400 | 1500 |

3^{−} | 3200 | 2000 |

4^{+} | 4000 | 2500 |

5^{−} | 4800 | 3000 |

6^{+} | 5600 | 3500 |

7^{−} | 6400 | 4000 |

8^{+} | 7000 | 4500 |

9^{−} | 7600 | 5000 |

10^{+} | 8000 | 5250 |

J^{p}
. | APH region . | Delves region . |
---|---|---|

0^{+} | 800 | 500 |

1^{−} | 1600 | 1000 |

2^{+} | 2400 | 1500 |

3^{−} | 3200 | 2000 |

4^{+} | 4000 | 2500 |

5^{−} | 4800 | 3000 |

6^{+} | 5600 | 3500 |

7^{−} | 6400 | 4000 |

8^{+} | 7000 | 4500 |

9^{−} | 7600 | 5000 |

10^{+} | 8000 | 5250 |

The Delves hyperspherical functions used in the Delves propagation step discussed above are computed using an efficient 1D Numerov method for the vibrational wave functions at each value of *j* and *l*. As mentioned in Sec. I, the rotational part of the wave functions is analytic. The number of rovibrational quantum numbers is determined from energy considerations. All the locally open channels and a significant number of closed channels must be included to ensure good unitarity. The Delves basis for *J* = 0 includes 500 rovibrational functions, which span nine vibrational levels of Li_{2}(*v* = 0–8) and eight for LiNa(*v* = 0–7) up to a total energy of ∼0.394 213 eV. Table II lists all the rotational quantum numbers used for each vibrational level. The rovibrational states computed for the asymptotic Jacobi basis are also listed. The number of rotational states for each vibrational level of Li_{2} is given by *j*_{max}/2 + 1 and (*j*_{max} + 1)/2 for even and odd symmetry, respectively. For LiNa, all values of *j* are included for each symmetry, so the number of states in each case is given by *j*_{max} + 1. From Table II, we see that by summing up these expressions for each Delves vibrational level, there are 200 even states for Li_{2} and 300 for LiNa giving a total of 500. Similarly, for the odd states, we find 196 for Li_{2} and 304 for LiNa. The Jacobi basis functions are computed in the same way as the Delves functions but satisfy the asymptotic Jacobi diatomic Schrödinger equation for each arrangement channel.^{37} The Jacobi rovibrational quantum numbers are also listed in Table II. The open Jacobi channels label the initial and final states of the scattering matrix $Sf\u2190iJ\u2009p\u2009q(E)$, where *i* = (*τ*_{i}, *v*_{i} *j*_{i}, *l*_{i}) and *f* = (*τ*_{f}, *v*_{f} *j*_{f}, *l*_{f}) label the initial and final quantum numbers: *τ* is the arrangement channel, *v* and *j* are the diatomic vibrational and rotational quantum numbers, and *l* is the orbital angular momentum. In this work, *τ*_{i} and *τ*_{f} label the reactant LiNa and product Li_{2} arrangement channels, respectively. In Sec. III, the product Li_{2} quantum numbers will also be denoted by primes *v*′ and *j*′. Asymptotically, for *J* = *l* = 0, there are 65 and 63 open Jacobi channels for *E*_{c} ≤ 1 K and even and odd symmetry, respectively. These open channels include the even and odd states of the initial reactant channel LiNa(*v* = 0, *j* = 0). The rovibrational quantum numbers for the open Li_{2} channels are listed in the third column of Table II. We note that the number of Delves channels also increases approximately linearly with *J* (see Table I), but the number of Delves channels is much smaller than the number of APH channels (the large number of APH channels is due to the deep potential well in the interaction region).

^{a}Li_{2}
. | Delves v, j_{max}
. | Jacobi v, j_{max}
. | Open v, j_{max}
. | ^{b}LiNa
. | Delves v, j_{max}
. | Jacobi v, j_{max}
. |
---|---|---|---|---|---|---|

0, 62(63) | 0, 68(69) | 0, 40(41) | 0, 62(63) | 0, 72(72) | ||

1, 58(59) | 1, 64(65) | 1, 34(35) | 1, 57(57) | 1, 67(68) | ||

2, 54(55) | 2, 60(61) | 2, 28(27) | 2, 51(52) | 2, 63(63) | ||

3, 50(49) | 3, 56(57) | 3, 18(17) | 3, 45(45) | 3, 57(57) | ||

4, 44(45) | 4, 52(53) | 4, 37(38) | 4, 52(52) | |||

5, 40(39) | 5, 48(47) | 5, 28(28) | 5, 45(45) | |||

6, 34(33) | 6, 42(43) | 6, 13(14) | 6, 38(38) | |||

7, 26(25) | 7, 36(37) | |||||

8, 14(15) | 8, 30(31) |

^{a}Li_{2}
. | Delves v, j_{max}
. | Jacobi v, j_{max}
. | Open v, j_{max}
. | ^{b}LiNa
. | Delves v, j_{max}
. | Jacobi v, j_{max}
. |
---|---|---|---|---|---|---|

0, 62(63) | 0, 68(69) | 0, 40(41) | 0, 62(63) | 0, 72(72) | ||

1, 58(59) | 1, 64(65) | 1, 34(35) | 1, 57(57) | 1, 67(68) | ||

2, 54(55) | 2, 60(61) | 2, 28(27) | 2, 51(52) | 2, 63(63) | ||

3, 50(49) | 3, 56(57) | 3, 18(17) | 3, 45(45) | 3, 57(57) | ||

4, 44(45) | 4, 52(53) | 4, 37(38) | 4, 52(52) | |||

5, 40(39) | 5, 48(47) | 5, 28(28) | 5, 45(45) | |||

6, 34(33) | 6, 42(43) | 6, 13(14) | 6, 38(38) | |||

7, 26(25) | 7, 36(37) | |||||

8, 14(15) | 8, 30(31) |

^{a}

The *j*_{max} value indicates the maximum value of *j* = 0, 2, 4, …, *j*_{max} for even symmetry and the value in parenthesis indicates the maximum value for odd symmetry *j* = 1, 3, 5, …, *j*_{max}.

^{b}

The *j*_{max} value indicates the maximum value of all *j* = 0, 1, 2, …, *j*_{max}.

## III. RESULTS

The full-dimensional numerically exact *ab initio* based quantum reactive scattering results for the Li + LiNa(*v* = 0, *j* = 0) → Li_{2}(*v*′, *j*′) + Na reaction are presented below. The calculations include all the partial waves needed to obtain fully converged state-to-state cross sections at 56 collision energies between 1 nK and 1 K. The rate coefficients are computed by multiplying the cross sections *σ*_{fi} by the relative collision velocity *v*: *K*_{fi} = *v σ*_{fi}. The exchange symmetry associated with the identical ^{6}Li nuclei is rigorously treated by using the appropriately symmetrized quantum mechanical wave functions. All the results include the appropriate nuclear spin statistical factors: 2/3 and 1/3 for even and odd exchange symmetry, respectively. We note that asymptotically, the even (odd) exchange symmetry correlates with the even (odd) rotational levels of Li_{2}. The total rate coefficients are presented in Sec. III A, and convergence with respect to the partial wave sum is demonstrated. A representative set of rotationally resolved rate coefficients are presented in Sec. III B. Notable shape resonances and other features are identified and the primary partial wave(s) that are responsible for the features are indicated. The angular distributions or DCSs are presented in Sec. III C as a function of both the collision energy and scattering angle. Significant oscillations and other striking features (i.e., sharp spikes, deep channels, and ridges) due to quantum interference and shape resonances are observed.

### A. Total rate coefficients

The total rate coefficients for the Li + LiNa(*v* = 0, *j* = 0) → Li_{2} + Na reaction summed over all product *v*′ and *j*′ states of Li_{2} are plotted in Figs. 2(a) and 2(b) for even and odd exchange symmetry, respectively. The individual contributions from each partial wave (*l* = 0–10) are also plotted in color: red (blue) for the even (odd) partial waves. The 56 collision energies are indicated by the black triangles superimposed on the *l* = 0 (red) curve. The total rate coefficients (summed over all 11 partial waves) are the thick black solid curves in panels (a) and (b). The rate coefficients for each partial sum are also indicated by the thin solid black curves. The plots in Fig. 2 show that summing over the 11 partial waves (i.e., *l* = 0–10) are sufficient for obtaining fully converged results up to 1 K. It is clear that the Wigner regime where only the *s*-wave (i.e., *l* = 0) contributes occurs for collision temperatures below 10^{−4} K. At these ultracold temperatures, the rate coefficients approach a constant (non-zero) value. We note that the rate coefficients of odd symmetry are approximately a factor of two smaller than the even ones due to the factor of two difference in the nuclear spin statistical factors. The total rate coefficients increase more or less monotonically with the collision energy. A small but noticeable undulation (bump) is seen near 50 mK due to an *l* = 3 shape resonance (see Sec. III B). The inset in each panel of Fig. 2 plots a zoomed in view of the total rate coefficient near the resonance (note the linear scale in both the x and y axes).

### B. Rotationally resolved rate coefficients

A representative set of rotationally resolved rate coefficients are plotted in Figs. 3 and 4 for even and odd exchange symmetry, respectively. Each of the four panels (a)–(d) in these two figures corresponds to the product Li_{2} vibrational states *v*′ = 0, 1, 2, and 3, respectively. The product Li_{2} rotational quantum number *j*′ is indicated for each curve. As expected, the rate coefficients approach a constant value in the Wigner regime where only *s*-wave scattering occurs (i.e., for collision energies below about 0.1 mK). As shown in several previous ultracold studies, the ultracold rate coefficients can be dramatically enhanced or suppressed due to quantum interference effects.^{23,25,36} Thus, the smallest and largest ultracold rate coefficients in Figs. 3 and 4 are the result of significant destructive and constructive interference, respectively. This quantum interference occurs between two primary reaction pathways: (1) a direct pathway and (2) a looping pathway that encircles the conical intersection [see the red dot in Fig. 1(a)]. Due to the unique properties of ultracold collisions, the relative phase shift of the scattering amplitudes associated with these two pathways approaches an integral multiple of *π*.^{23,25} Thus, the quantum interference can approach its maximum constructive or destructive values if the magnitudes of these two scattering amplitudes are comparable.^{23,25} At higher collision energies, the rate coefficients typically increase in magnitude and notable features such as bumps, steps, and shoulders are observed due to underlying shape resonances associated with the various partial waves. Several of the most notable features are labeled by the partial wave(s) responsible. For example, in panel (a) of Fig. 3, a prominent *l* = 3 shape resonance occurs for *j*′ = 14 (black long-short dashed curve) near 50 mK. A series of two shape resonances due to *l* = 2 and 5 appear for *j*′ = 2 (black dashed curve) near 20 mK and 0.2 K, respectively. A prominent *l* = 2 shape resonance also occurs for *j*′ = 24 (red solid curve) in the 20–40 mK range. Shape resonance features can also be seen in panels (b)–(d) of Fig. 3. In particular, a high energy resonance feature is predicted for *j*′ = 0 in both panels (c) and (d) near 0.6 K due to *l* = 6, 7, and 8. Highlighting a few notable features for odd symmetry, we see in panel (a) of Fig. 4 that a notable shape resonance occurs for *j*′ = 13 (black dashed-dotted curve) due to the *l* = 2 partial wave near 20 mK. In panel (b) of Fig. 4, a series of two *l* = 1 and *l* = 3 shape resonances are seen for *j*′ = 3 (black dashed curve) near 5 and 50 mK, respectively. An *l* = 3 shape resonance is also seen for *j*′ = 15 (black dashed-dotted curve) and a *l* = 2 resonance occurs for *j*′ = 21 (red solid curve). As shown in Fig. 3, a high energy resonance feature is predicted near 0.6 K for *v*′ = 3 and *j*′ = 1 [black solid curve in panel (d)] due to *l* = 7 and 8.

### C. Differential cross sections

The appearance of intriguing quantum phenomena (i.e., interference and shape resonances) reported above in the rotationally resolved rates coefficients become even more noticeable in the angular distributions or DCSs. The integration of the DCSs over the scattering angle in computing the integral cross sections and rate coefficients averages out a significant portion of these quantum effects. Figure 5 plots the DCSs for *v*′ = 0 and even symmetry for the same product rotational states presented in panel (a) of Fig. 3. That is, panels (a)–(f) in Fig. 5 correspond to *j*′ = 0, 2, 14, 24, 34, and 36, respectively. The DCSs are plotted as a function of both collision energy and scattering angle. The scattering angle is the relative angle between the reactant Li–LiNa and product Na–Li_{2} **k** vectors. Thus, forward and backward scattering correspond to *θ* = 0° and 180°, respectively. A quick scan over the six panels immediately reveals fascinating features in the DCSs that are unique to each product state: high-frequency choppy oscillations are seen in panel (a), sharp downward spikes in the forward scattering region are seen in panels (b) and (c), also notable bumps are seen in the backward regions in panels (b) and (c) due to the *l* = 2 and *l* = 3 shape resonances, a deep channel (due to the *l* = 1 partial wave) is seen in panel (d) that starts in the forward scattering region at low energies and curves into larger angles at higher energies, a significant dip is seen in the backward scattering region in panels (e) and (c) due to *l* = 1 and (a) due to *l* = 3, and a very large downward spike is seen in the forward scattering region in panel (f) due to the *l* = 1 partial wave. These features are due to constructive and destructive quantum interference between the two reaction pathways mentioned above, shape resonances, and quantum interference between the partial waves themselves.

The DCSs for *v*′ = 1 and even symmetry are plotted in Fig. 6, and the product rotational states correspond to those in panel (b) of Fig. 3. That is, panels (a)–(f) in Fig. 6 correspond to *j*′ = 0, 2, 10, 14, 16, and 34, respectively. Again, we see unique structure in the DCSs for each product state: in panel (a), a notable bump in the backward region associated with the *l* = 2 shape resonance, a depression near *θ* = 130° and 20 mK is also seen in panel (a) due to primarily *l* = 1 and 4, a broad symmetric ridge (bump) is seen in panel (b) in both the forward and backward regions associated with an *l* = 1 shape resonance, a notable dip due to *l* = 1 and a ridge due to *l* = 2 and 3 are seen in the backward region in panel (c), several sharp downward spikes are seen in the forward region in panels (d) and (e) (the lowest energy spikes are due to *l* = 1), and a broad channel due to *l* = 1 and a ridge due to the *l* = 2 shape resonance are seen in the forward region in panel (f).

The DCSs for even symmetry with *v*′ = 2 and *j*′ = 0, 8, and 16 are plotted in the panels (a)–(c) of Fig. 7. These three product states correspond to those plotted in panel (c) of Fig. 3. Two bumps are seen in the backward region of panel (a). The broad bump at a low energy (near 5 mK) is due to the *l* = 1 shape resonance and the narrow bump at a high energy (≈0.2 K) is due to the overlapping *l* = 6 and 7 shape resonances. The low energy channel feature in the forward region in panel (b) near 0.2 mK is due to *l* = 1 and the sharp downward spike at 50 mK is due to *l* = 3. The broad dip at a low energy in the forward region in panel (c) is due to *l* = 1 and the sharp dip in the backward region is due to *l* = 2. The bump in the backward region of panel (c) is due to *l* = 3. The DCSs for even symmetry with *v*′ = 3 and *j*′ = 0, 8 and 18 are plotted in panels (d)–(f) of Fig. 7. These three product states correspond to those plotted in panel (d) of Fig. 3. The first shallow dip at a low energy in the forward region in panel (d) is from *l* = 1 and the sequence of ripples extending to higher energies corresponds to several higher partial waves. The first dip and narrow downward spike in the forward region of panel (e) correspond to *l* = 1 and 3, respectively. In panel (f), the first and second downward spikes correspond to *l* = 2 and 3, respectively.

## IV. CONCLUSIONS

In summary, a recently developed first-principles based *ab initio* PES for the adiabatic ground electronic state of Li_{2}Na was used to perform numerically exact quantum reactive scattering calculations for the Li + LiNa(*v* = 0, *j* = 0) → Li_{2}(*v*′, *j*′) + Na reaction. This PES includes an accurate treatment of the long-range three-body interactions and spectroscopically accurate asymptotic diatomic potential curves.^{36} A time-independent coupled-channel methodology based on hyperspherical coordinates was used to perform the full-dimensional scattering calculations at 56 collision energies spanning nine orders of magnitude from ultracold (1 nK) to cold (1 K) temperatures. The scattering calculations include all the partial waves necessary to obtain fully converged cross sections (with respect to the partial wave sum) over the full energy range for the first time. The calculations also include a rigorous quantum mechanical treatment of the identical particle exchange symmetry due to the two identical ^{6}Li nuclei.

The total and rotationally resolved rate coefficients are computed as a function of collision energy for both even and odd exchange symmetries. The total rate coefficients approach a constant (non-zero) value in the ultracold limit and increase monotonically with the increase in collision energy. A small bump is observed in the total rate coefficients for collision energies near 50 mK due to an *l* = 3 shape resonance. The rotationally resolved rate coefficients exhibit more structure and several bumps and other features are observed as a function of energy due to the contributions from various partial wave shape resonances. The unique properties associated with ultracold collisions leads to significant quantum interference effects that ultimately determine the asymptotic value of the rotationally resolved rate coefficients as the collision energy approaches absolute zero. Particular focus of the present work is on the investigation of the angular distributions or DCSs as a function of both the collision energy and scattering angle. Significant quantum interference and resonance structure are observed in the DCSs, which are unique to each product Li_{2} rovibrational state. These unique “quantum fingerprints” encode the underlying quantum dynamics and molecular interactions (e.g., potential wells, multiple reaction pathways, rotational anisotropy, and long-range interactions) and are experimentally measurable. Although we expect that improvements to the PES will be required to obtain quantitative agreement with future state-resolved experiments, the current theoretical study clearly demonstrates that fascinating quantum phenomena and structure will be present in product state-resolved experimental cross sections. In future theoretical studies, we plan to investigate the effects of the excited electronic state and its associated geometric phase and non-adiabatic coupling on the angular distributions. This will require the use of a two-state diabatic approach.^{30,36}

The quantum effects reported here are unique to the ultracold/cold energy regime and provide an unprecedented opportunity for probing and understanding fundamental molecular interactions, reaction mechanisms, and quantum phenomena. One can envision exploiting these quantum effects to produce extremely sensitive sensors and measuring devices. They also provide a potential avenue for the realization of quantum control technologies via the selection of a particular initial and/or final state, the preparation of a specific orientation between the colliding partners (i.e., stereodynamics), or the application of external fields. We hope that this work will help stimulate future theoretical and experimental studies on the fascinating energy regime below 1 K.

## SUPPLEMENTARY MATERIAL

See the supplementary material for plots of the differential cross sections of the odd exchange symmetry.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

This work was performed under the auspices of the U.S. Department of Energy under Project No. 20170221ER of the Laboratory Directed Research and Development Program at Los Alamos National Laboratory. This research used resources provided by the Los Alamos National Laboratory Institutional Computing Program. Los Alamos National Laboratory is operated by Triad National Security, LLC, for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001).

## REFERENCES

*ab initio*programs,