Non-adiabatic quantum interference in the ultracold Li + LiNa - Li 2 + Na reaction

Electronically non-adiabatic eﬀects play an important role in many chemical reactions. However, how these eﬀects manifest in cold and ultracold chemistry remains largely unexplored. Here for the first time we present from first principles the non-adiabatic quantum dynamics of the reactive scattering between ultracold alkali-metal LiNa molecules and Li atoms. We show that non-adiabatic dynamics induces quantum interference eﬀects that dramatically alter the ultracold rotationally resolved reaction rate coeﬃcients. The interference eﬀect arises from the conical intersection between the ground and an excited electronic state that is energetically accessible even for ultracold collisions. These unique interference eﬀects might be exploited for quantum control applications such as a quantum molecular switch. The non-adiabatic dynamics are based on full-dimensional ab initio potential energy surfaces for the two electronic states that includes the non-adiabatic couplings and an accurate treatment of the long-range interactions. A statistical analysis of rotational populations of the Li 2 product reveals a Poisson distribution implying the underlying classical dynamics are chaotic. The Poisson distribution is robust and amenable to experimental verification and appears to be a universal property of ultracold reactions involving alkali metal dimers.


Introduction
Ultracold molecules and in particular ultracold polar molecules are at the forefront of precision spectroscopy, sensing, controlled studies of chemical reactions, quantum many-body physics, and quantum computing. [1][2][3][4][5][6][7][8][9][10] Polar molecules comprised of heteronuclear alkali metal dimers such as KRb, NaK, NaRb and LiNa have attracted considerable attention in recent years in controlled studies of chemical reactions. 2,5,6,9,10 Electronically non-adiabatic effects are expected to play an important role in atom-dimer reactions involving these molecules. The reactions proceed along a barrierless reaction pathway into a deep attractive potential well. A conical intersection (CI) occurs between the ground electronic state and the first excited doublet electronic state within the attractive well region and this CI is energetically accessible even for collision energies in the ultracold limit for ground state reactants. Thus, a non-adiabatic quantum mechanical treatment is required that includes both electronic states. Explicit quantum calculations for these reactions remain a formidable challenge even for dynamics on a single Born-Oppenheimer adiabatic electronic potential energy surface (PES). [11][12][13] However, calculations on a single adiabatic PES are not accurate for these systems even if a generalized vector potential based approach 14,15 is used to include the geometric phase (GP) 14,16 associated with the CI. Fortunately, we have recently developed a new non-adiabatic quantum reactive scattering methodology that has made it possible to rigorously treat ultracold reactions occurring on two coupled electronic states for the first time. 17 This work presents the first rigorous treatment of non-adiabatic chemistry in the ultracold regime for systems of current experimental interest, advancing the state-of-the-art in the theoretical treatment of ultracold chemistry to new frontiers.
Non-adiabatic dynamics is most often associated with collisions at high (thermal) energies where excited electronic states are typically more accessible. At these high energies quantum effects are much less pronounced and the dynamics of heavy nuclei is often treated classically. In contrast, at ultracold collision energies a fully quantum mechanical treatment is required even for heavy nuclei and fundamentally new reaction mechanisms can also come into play. The unique properties of ultracold collisions therefore provide a novel testing ground for exploring and understanding non-adiabatic quantum dynamics. Quantum phenomena (interference, Fermi-Bose statistics, phase shift quantization, etc.) are dramatically enhanced at ultracold temperatures and the delicate interplay between them can lead to unprecedented sensitivity and control of the collision outcome. Alkali metal systems in particular provide an ideal test bed for ultracold experimental and theoretical studies of non-adiabatic effects and their possible technological applications.
In this work, we present a first principles full-dimensional quantum dynamics study of non-adiabatic effects in the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 , j 0 ) + Na reaction. The methodology also includes a rigorous treatment of the identical particle exchange symmetry between the two 6 Li nuclei. The rotationally resolved rate coefficients are computed as a function of collision energy from 1 nK to 10 K using a coupled two-state diabatic electronic representation. [17][18][19] The non-adiabatic results are compared to a conventional Born-Oppenheimer calculation based on a single adiabatic electronic PES. Both of these calculations are also compared to a universal model which is based on a simple one-dimensional reaction path consisting of a long-range van der Waals (C 6 ) potential. 20 Quantum interference between the two reaction pathways which encircle the CI is shown to significantly enhance or suppress the rate coefficients at ultracold collision energies (E c o 1 mK). The GP which is included in the non-adiabatic calculations reverses the nature of the quantum interference from constructive to destructive and vice versa. [21][22][23] Thus, the non-adiabatic ultracold rate coefficients are significantly enhanced or suppressed relative to the conventional Born-Oppenheimer rate coefficients when quantum interference effects are significant. The quantum interference effects lead to strong fluctuations in the rotationally resolved rate coefficient distributions. A statistical analysis of these fluctuations reveals that they are Poissonian which is consistent with an underlying classically chaotic dynamics. 12,13 The Poisson distributions are shown to be robust with respect to variations in the PES and chemical system and therefore appear to be a universal property of these types of reactions that proceed through a potential well. The quantum dynamics calculations are based on fulldimensional accurate ab initio electronic PESs which are computed for both the ground and first excited states for the first time. A state-of-the-art electronic structure code (MOLPRO) is used to compute the electronic PESs and the non-adiabatic coupling elements. 24 Also included in the PESs is an accurate experimentally based treatment of the long range (van der Waals) interactions and diatomic potentials.
The paper is organized as follows: in Section 2 the nonadiabatic quantum dynamics methodology is reviewed. Section 3.1 presents the Born-Oppenheimer electronic PESs for the ground and first excited electronic states that include the prominent energetically accessible CI. The energetics and reaction path are discussed for the reaction as it proceeds from reactants Li + LiNa to products Li 2 + Na, and notable features of the PESs are also emphasized. The rotationally resolved rate coefficients are then presented in Section 3.2 and the underlying quantum interference mechanisms are discussed in detail. Finally, in Section 3.3 a statistical analysis of the rotational distributions is discussed and the resulting Poisson distributions are presented. The conclusions are presented in Section 4 and additional analysis, figures, and the PES scaling sensitivity studies are included in the Appendix.

First principles non-adiabatic quantum dynamics
First principles non-adiabatic quantum dynamics calculations require solving a generalized Born-Oppenheimer equation for the nuclear motion which includes both the ground and first excited electronic states. That is, the total molecular wavefunction |Ci is expanded as jCi ¼ the diabatic electronic basis functions and the expansion coefficientsc i are the nuclear motion wave functions. The nuclear coordinates relative to the center-of-mass of the molecule are collectively denoted by the six-dimensional vector x which includes the three interatomic degrees of freedom and three rotational (Euler) angles. Different choices for the nuclear coordinates can be used as will be discussed in detail below. For numerical reasons it is advantageous to use a diabatic electronic representation for which the PESs are well-behaved smooth functions of the nuclear coordinates. In addition, the coupling between the two electronic states appears simply as an off-diagonal element in the potential energy matrix (PEM) (i.e., there are no derivative coupling terms in the kinetic energy operator). 17 The relevant diabatic Schrödinger equation for the nuclear motion is given bŷ where the first term in brackets is the diagonal kinetic energy operator for the nuclear motion with matrix elements T = Àh 2 r 2 /(2m),r denotes derivatives with respect to x and m is three-body reduced mass. The second term is the PEM which couples the two diabatic electronic states and is a function of only the three interatomic degrees of freedom (x). That is, each matrix element (i.e., the V d 11 , V d 22 and V d 12 = V d 21 ) of the PEM is a three-dimensional PES. We note that in contrast to the 2 Â 2 diabatic Schrödinger equation given in eqn (1), a conventional Born-Oppenheimer quantum dynamics calculation solves the adiabatic single surface Schrödinger equation where V a 1 denotes the adiabatic ground electronic state PES (which can be expressed in terms of the diabatic PEM elements, see Section 2.1).
A first-principles solution of eqn (1) requires two primary steps: (1) a quantum mechanical ab initio calculation of the diabatic PEM, and (2) a numerically exact quantum mechanical solution of the six-dimensional diabatic Schrödinger equation for the nuclear motion. Each of these steps are highly nontrivial and require extensive calculations as discussed in the following Sections 2.1 and 2.2, respectively.

Ab initio potential energy surface calculations
The PESs of the LiNaLi collisional complex were not available from the literature and their computation required substantial effort due to the complexity of the multi-electron open-shell systems. We have computed the two energetically-lowest nonrelativistic doublet adiabatic potentials, 2 A 0 and 2 B 0 , with the MOLPRO electronic structure package. 24 The core electron shells of Li and Na are described by the Stuttgart/Cologne energy-consistent, single-valence electron, relativistic pseudopotentials ECP2SDF and ECP10SDF 25 leaving only three valence electrons in the active space for explicit treatment. The polarization of the effective cores and residual core-valence correlations are modeled via the l-independent core polarization potential (CPP) with Müller-Meyer damping functions. 26 The CPP parameters, i.e. static dipole polarizabilities of the atomic cores are taken from ref. 27. Cutoff functions with exponents 0.95 a.u. and 0.82 a.u. are employed for Li and Na, respectively (here a.u. stands for atomic unit). Basis sets from ref. 28 describe the three valence electrons, specifically, uncontracted sp basis sets augmented by additional s, p, d and f polarization functions are used for both Li and Na. The multi-configurational self-consistent field (MCSCF) method 29,30 is first used to obtain configuration state functions (CSFs). A multi-reference configuration interaction (MRCI) calculation 31,32 is then performed using a large active space constructed from the CSFs, giving the three-dimensional adiabatic surfaces of the two lowest energy states for LiNaLi, V a 1 and V a 2 . The adiabatic PESs are computed as functions of the separation between Na and first Li atom, R 1 , the separation between Na and second Li atom, R 2 , as well as the angle y between R 1 and R 2 . We have computed the potentials on the ten angular grid points y = 101, 201, 401, 601, 801, 1001, 1201, 1401, 1601, 1781, and radial R 1 and R 2 grids from R min = 3.2a 0 to 10a 0 in steps of 0.2a 0 and from 10a 0 to R max = 19a 0 in steps of 0.3a 0 with the constraint that the separation between the two Li atoms R 3 ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi R 1 2 þ R 2 2 À 2R 1 R 2 cos y p 4 R min (a 0 = 0.0529177 nm is the Bohr radius). In total, more than 40 000 points were calculated.
The accuracy of the LiNaLi trimer potentials has been determined from a comparison of the singlet X 1 S + and triplet a 3 S + potentials of LiNa and Li 2 determined with MOLPRO using the same basis sets as for the trimer calculations and the spectroscopically-accurate potentials from ref. [33][34][35]. For example, our computed well depths were smaller by no more than hc Â 25 cm À1 for both the singlet and triplet dimer potentials corresponding to 0.3% and 9% of their well depths, respectively. Here, h is the Planck constant and c is the speed of light in vacuum.
In order to study the reaction dynamics we need to diabatize the two adiabatic potentials V a 1,2 (R 1 ,R 2 ,y) of LiNaLi, as they have a conical intersection for a single geometry with C 2v symmetry. In the adiabatic representation, the electronic eigenstates |j a 1,2 i are coupled by non-adiabatic first-derivative couplings hj a 1 jrjj a 2 i with respect to the nuclear coordinates. The nonadiabatic coupling matrix element was computed for all of the grid points using MRCI theory and a finite difference method for the derivatives (the DDR procedure). These non-adiabatic derivative couplings are problematic at the CI where they become singular. Fortunately, the singular part of these couplings are removable. That is, the non-adiabatic coupling matrix elements are vector quantities and can be decomposed into longitudinal (removable) and transverse (non-removable or residual coupling) components. 36,37 By an appropriate unitary transformation from the adiabatic to diabatic electronic representation, the problematic singular (longitudinal) part of the derivative coupling at the CI can be removed. For two electronic states this can be achieved with the transformation where the real-valued mixing angle b is a function of R 1 , R 2 , and y and the two diabatic electronic wavefunctions |j d 1,2 i are independent of these three coordinates. Hence, the non-adiabatic coupling between the two adiabatic states can be expressed as is then found assuming that b = 0 for all geometries with C 2v symmetry, i.e. R 1 = R 2 , and spatial integration of eqn (4). 38 The signs of hj a 1 jrjj a 2 i as computed within MOLPRO are arbitrary from geometry to geometry. We chose the convention that leads to a smooth b. The electronic potential matrix in the diabatic basis is the 2 Â 2 matrix where we have suppressed all R 1 , R 2 , and y dependences for clarity. Finally, we fit and interpolate the diabatic PESs V d ij (R 1 ,R 2 ,y) for i, j = 1, 2 with a two step procedure. Using the dimer-inmolecule theory, 39 we separate the diabatic PESs into a pairwise and a three-body component with The pairwise components V pw ij (R 1 ,R 2 ,y) are weighted sums of the spectroscopically-accurate X 1 S + and a 3 S + potentials of NaLi and Li 2 . [33][34][35] The weights follow from the realization that the electron angular momenta of the individual Li and Na atoms in their electronic ground state are h /2 (or spin-1/2) and that the two diabatic trimer electron basis functions have a total electron angular momentum equal to h /2 for doublet states, after coupling the three spin-1/2 atomic electron angular momenta (h is the reduced Planck constant). The three-body component V 3B ij (R 1 ,R 2 ,y) is fit using the reproducing kernel Hilbert space (RKHS) technique [40][41][42] and relies on an expansion in terms of Legendre polynomials P l (cos y) as functions of cos y. That is, where we expand up to l = 8 and the coefficients A l,ij (R 1 ,R 2 ) have analytical representations with as many linear parameters as the combined number of R 1 and R 2 grid points.
In order to extrapolate the three-body component V 3B ij (R 1 ,R 2 ,y) for geometries where one or both separations R 1 or R 2 are larger than R max , we combine two strategies. First, we use the fact that the three-body component of the adiabatic potentials V a i fall off as 43 when all three pair separations are large. Here, R 3 is the separation between the two Li atoms and the three-body dispersion coefficient C 9 still depends on the shape of triangle with sides R 1 , R 2 , and R 3 . Therefore, we extended the (R 1 ,R 2 ,y) grid of the diabatic potentials V 3B ij (R 1 ,R 2 ,y) to all geometries with the following procedure: for R 1 4 R max and fixed grid points R 2 and y the C 9 coefficients are determined from V 3B ij (R max ,R 2 ,y); for R 2 4 R max the C 9 coefficients are determined from V 3B ij (R 1 ,R max ,y); and when both R 1 and R 2 are larger than R max the C 9 coefficients are determined from V 3B ij (R max ,R max ,y). We used a step size of 0.3a 0 from R max up to 50a 0 for both R 1 and R 2 . In fact, the RKHS procedure described in the previous paragraph has been applied to this extended data set. Finally, the second part of our strategy to extrapolate to large separations is to multiply the three-body component with a switching function with R c = 30a 0 . This function is unity when the three atoms are close to each other but switches smoothly to zero when the separation between any two atoms is large.

Non-adiabatic quantum dynamics
The numerically exact quantum mechanical solution of the 2 Â 2 diabatic Schrödinger eqn (1) is based on a timeindependent coupled-channel (CC) approach. 17,44,45 Adiabatically adjusting Principal axis Hyperspherical (APH) coordinates are used in the interaction region and Delves hyperspherical coordinates are used in the long-range asymptotic region. 17,[44][45][46] The hyperradius r is common to both coordinate systems which facilitates the coordinate transformation from the APH to Delves at an intermediate value of r = r m (determined by numerical convergence studies, see below). Furthermore, hyperspherical coordinates have the advantage of treating all reactant and product channels simultaneously and allow for a natural numerical separation of the problem into radial and angular parts. We choose the body frame (BF) z-axis perpendicular to the plane of the triatomic molecule and the BF x and y axis are chosen to lie along the instantaneous principal axes of inertia (i.e., the Q and q vectors of Pack and Parker, respectively). 44 The collective set of six nuclear coordinates are denoted by x = (x,x) where the three internal coordinates are x = (r,y,f) and x = (a,b,g) (the three Euler angles which specify the orientation of the BF relative to the space fixed frame). The hyperradial coordinate r and hyperangles y and f correspond to a symmetric stretch, bending, and pseudorotational motion, respectively. 47 These three coordinates can be expressed explicitly in terms of the three internuclear distances (see ref. 44 for details). The kinetic energy operator in hyperspherical coordinates can be expressed aŝ where the first term is the radial part and the second term is the angular part. The operatorL is the grand angular momentum operator, and the angular part of eqn (13) is given bŷ where the polar angleỹ = p À 2y, and J x , J y , J z denote the BF components of the total angular momentum operator J. For ultracold collision energies we only need to consider the single orbital angular momentum l = 0 (i.e., s-wave scattering). Since the reactant LiNa is in its ground rotational state j = 0, this corresponds to total angular momentum J = l + j = 0. The potential matrix strongly couples the radial (r) and angular (L, f) coordinates so that eqn (1) and (2) are nonseparable. However, we can perform a numerical separation by dividing the hyperradius r into many ''sectors'' and solving the angular part of the problem with r fixed at the center of each sector. The angular solutions within each sector provide an optimal localized fixed basis for solving the CC radial equations within that sector. The global radial solution is obtained by propagating a set of CC equations from small to large r across each sector. At the boundaries of each sector, the radial solutions are transformed from one sector basis to the next using the overlap matrix between the angular functions computed at the two sectors. Specifically, we expand the total diabatic nuclear motion wave function C ¼c 1 c 2 in eqn (1) in terms of the angular functions (F) for a given sector r = r x as where i and t denote the CC indices, p denotes inversion parity, q denotes the particle exchange symmetry which is relevant for triatomic molecules with identical nuclei, and w (ỹ,f,a,b,g) denotes the angular coordinates of the 5D hypersphere. The radial coefficients z Jpq it are computed numerically from the propagation of the CC radial equations where the potential coupling matrix is given by The F JMpq t are the 5D hyperspherical functions which are solutions to the angular equation at sector r = r x The V d in eqn (17) and (18) is the 2 Â 2 diabatic PEM given in eqn (5)- (8). Once the angular solutions to eqn (18) are computed, they can be used to compute the overlap matrices O tt 0 = hF JMpq t (w;r x )|F JMpq t (w;r x+1 )i and potential coupling matrices defined in eqn (17) at each sector. These are then used to solve the CC radial eqn (16) using Johnson's log-derivative propagator. 48,49 At the matching radius r = r m , the radial solution is transformed into the Delves basis using the overlap matrix between the APH and Delves functions computed at r = r m . The radial propagation in r is then continued to the final asymptotic value of r = r f using an equivalent set of Delves CC radial equations (see ref. 46 for the Delves equations). The Delves basis functions at each sector are computed separately for each diatomic arrangement channel (i.e., for r 4 r m the channels are decoupled). They consist of ro-vibrational wave functions computed numerically using a one-dimensional Numerov propagator for the vibrational motion and a set of analytic spherical harmonics for the rotational part. The Delves basis functions are used to compute the potential coupling matrices and overlap matrices between each sector analogous to the APH propagation. Finally, at the last asymptotic value of r = r f (determined from convergence studies, see below), the Delves radial solution is projected onto the analytic asymptotic radial solution in Jacobi coordinates. From this projection, the full state-to-state scattering (S) matrix is computed for all open reactant and product channels for total energy E. Once the S matrix is computed, cross sections s fi and rate coefficients K fi = vs fi (where v is the relative collision velocity) can be computed using standard expressions. 44 We note that the f,i denote the collective final and initial quantum numbers of the diatomic products and reactants (i.e., f = (t 0 ,v 0 , j 0 ) and i = (t,v, j) where t 0 and t denote the diatomic arrangement channel Li 2 or LiNa), respectively.
In the interaction region (r o r m ), the (5D) surface function Hamiltonian matrix in eqn (18) is diagonalized on a discrete grid in r at 144 logarithmically spaced points between (r i = 6.0a 0 and r m = 33.03a 0 ). The 5D basis functions consist of a hybrid FBR (Finite Basis Representation) in f, a DVR (Discrete Variable Representation) in y, and symmetrized Wigner D functions in a,b,g. 17,45,50 The size of the DVR and FBR basis is specified by the values of the integers ł max and m max , respectively. 17,45 The basis size increases with r and is determined from numerical convergence studies. In this study, four regions in r were used: 6.0 r r r 11.02a 0 , 11.02 o r r 19.54a 0 , 19.54 o r r 23.37a 0 , and 23.37 o r r 33.03a 0 . The corresponding sets of l max and m max are: l max = 99, 119, 127, and 127; and m max = 220, 240, 260 and 280, respectively. For total J = 0 considered in this study, the size of the angular Hamiltonian matrix in each region is given by (l max + 1) Â (2m max + 1) = 44 100, 57 720, 66 688 and 71 808, respectively. Fortunately, the size of these matrices can be dramatically reduced by using the SDT (Sequential Diagonalization Truncation) technique. 50 After applying SDT the dimension of the largest angular Hamiltonian matrix in each of the four regions is 13 008, 9412, 8401, and 6995. An efficient parallel numerical eigensolver (PARPACK) is then used to numerically diagonalize the sparse Hamiltonian matrices constructed at each sector. 51 As mentioned above, the set of angular eigensolutions of eqn (18) form a basis for the one-dimensional CC propagation in r (see eqn (16) and (17)). The number of CCs propagated in this work was 820 in the APH region and 500 in the Delves region. In the Delves region, 338 uniformly spaced r values were used between r m = 33. We note that all of the above quantum dynamics calculations rigorously treat the identical particle exchange symmetry between the two identical 6 Li nuclei. 17 Thus, the calculations are repeated separately for each exchange symmetry, even and odd.

Potential energy surfaces for Li 2 Na
The Born-Oppenheimer electronic PESs are plotted in Fig. 1 for both the ground and first excited electronic states of the Li 2 Na molecule. The PESs in Fig. 1 are two-dimensional slices plotted for a fixed Li 2 bond length of 6.25a 0 (close to its equilibrium bond length) and show the topology of the effective interaction potential experienced by the Na nuclei in the vicinity of Li 2 . Notable features include the two deep attractive wells (blue colored regions) on the ground state surface (black contours) and the inverted cone of the excited electronic state (red contours). All energies are reported relative to the bottom of the asymptotic potential well for the Li 2 + Na product channel. The minimum energy of the symmetric potential wells is À5814 K (see the thick solid black curve in Fig. 2). The ground and excited state PESs exhibit a conical intersection for Tshaped (i.e., C 2v ) geometries (see Fig. 1 inset). The minimum energy of the conical intersection is À3140 K (see the thick solid red curve in Fig. 2). The asymptotic energy of the Li + LiNa(v = 0, j = 0) reactant channel is shown by the thick black contour line at 2228 K (see also the thick horizontal dashed line in Fig. 2). From Fig. 1 we see that for ultracold collisions of Li with LiNa in its ground vibrational and rotational state, both the ground and excited electronic states are energetically accessible in the interaction region. Thus, both electronic states and the couplings between them must be included in the quantum dynamics calculations (see Section 2.2 for details). 17 These non-adiabatic couplings include the GP associated with the conical intersection. As discussed in detail in the following section, the GP can lead to a dramatic enhancement or suppression of the ultracold rotationally resolved rate coefficients. [21][22][23] We note that a traditional GP calculation 21-23 based on only the ground adiabatic electronic state is not applicable for this system since the CI is located below the energy of the incident channel. Thus, a non-adiabatic fully coupled two electronic state calculation is required. 17 Fig. 2 plots the APH surface function energies (i.e., the eigenvalues e Jpq t (r x ) of eqn (18)) for the ground (blue) and excited (red) electronic states as a function of r. For visualization purposes, the surface function energies were computed for each electronic state separately (i.e., on an uncoupled single adiabatic electronic state). The black and red thick solid curves plot, respectively, the minimum energies of the ground (V a 1 ) and excited (V a 2 ) adiabatic electronic state PESs as a function of r. Asymptotically for large r, the ground state adiabatic surface function energies (blue) approach the diatomic Li 2 (v 0 , j 0 ) and LiNa(v, j) rovibrational energies. The dashed horizontal black line denotes the energy of the Li + LiNa(v = 0, j = 0) reactant channel (2228 K) (see also the thick black contour line in Fig. 1). The excited electronic state is energetically open in the interaction region (small r) where the red curves drop below the Fig. 1 Ab initio Born-Oppenheimer potential energy surfaces (PESs) for Li 2 Na are plotted for a fixed Li 2 bond length of 6.25a 0 . The xy coordinates denote the location of the sodium nuclei (red sphere) relative to the center of the bond between the two Li nuclei (blue spheres). The ground (V a 1 ) and excited (V a 2 ) electronic state surfaces are contoured with black and red contours, respectively. Two attractive potential well regions (blue) are clearly visible on the ground state surface. The excited state surface exhibits a conical intersection with the ground state surface for T-shaped (C 2v ) geometries (see inset). The thick black contour line denotes the total energy 2228 K of the reactant Li + LiNa(v = 0, j = 0). The other 20 black contours lie between À5000 K and 5500 K inclusive. The 10 red contour lines lie between À2470 K and 5500 K inclusive. Both surfaces are not plotted above 5500 K. All energies are relative to the minimum energy of the asymptotic adiabatic ground electronic state of Li 2 + Na. black dashed horizontal line. The vertical series of short horizontal lines on the left edge of the plot denote the bound state energies of the cone states which are localized inside the cone (V a 2 ) of the excited adiabatic electronic state (see Fig. 1). The E and O label the bound states of even and odd exchange symmetry, respectively. These bound states were computed numerically in full-dimensionality using the adiabatic excited electronic state PES (V a 2 ) and could lead to Feshbach type scattering resonances. 52 The large solid black and red dots denote particular points of interest and are discussed in the Appendix. Fig. 3 and 4 plot the rotationally resolved rate coefficients for the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 = 0, j 0 ) + Na reaction as a function of collision energy from 1 nK to 10 K for even and odd exchange symmetry, respectively. Unless otherwise stated, all rate coefficients include the appropriate nuclear spin statistical factors of 2/3 and 1/3 for even and odd exchange symmetry (associated with the two identical 6 Li nuclei), respectively. At ultracold collision energies (o1 mK), only a single partial wave (i.e., l = 0 where l is the orbital angular momentum of Li about LiNa) contributes to the collision and the rate coefficient becomes finite (often referred to as the Wigner threshold regime). 53,54 The specific values of the ultracold rate coefficients are sensitive to the PES but values reported here are numerically exact for the PESs computed in this study. The curves in Fig. 3(a) and 4(a) are from the coupled two-diabatic electronic states calculation (2 Â 2) and the curves in Fig. 3(b) and 4(b) are from the calculation on a single adiabatic ground electronic state which does not include the GP (denoted as NGP for No GP). We see that many of the 2 Â 2 ultracold rate coefficients are significantly enhanced or suppressed relative to the NGP ones. The rotationally resolved rate coefficients for the other product vibrational levels v 0 = 1-3 show similar trends. A representative example is plotted in Fig. 5 for v 0 = 3 and j 0 = 5. The red curve in Fig. 5 is from the coupled two-diabatic electronic states calculation (2 Â 2) and the black curve is from the NGP calculation on a single adiabatic ground electronic state. We see that in the ultracold limit the 2 Â 2 rate coefficient (red) is significantly enhanced (by around 50 times) relative to the NGP one (black). As discussed in more detail below, the enhancement (suppression) is due to constructive (destructive) quantum interference between the direct and looping contributions to the total scattering amplitude. 21 The GP associated Fig. 2 The APH surface function energies are plotted as a function of the hyperradius r. The blue and red adiabats are computed on the electronic adiabatic ground (V a 1 ) and excited (V a 2 ) states, respectively. The thick black and red curves plot the minimum energy of the ground and excited adiabatic electronic states at each r, respectively. The horizontal black dashed line denotes the total energy of the Li + LiNa(v = 0, j = 0) reactant (2228 K). The blue adiabats on the right edge of the plot (i.e., large r) which lie below the black horizontal line correlate to the asymptotic ro-vibrational energies of the product Li 2 (v 0 , j 0 ) + Na. The series of energy levels labeled by the ''E'' (even exchange symmetry) and ''O'' (odd exchange symmetry) are the three-dimensional vibrational energies computed on the excited electronic state (i.e., ''cone states''). All energies are relative to the minimum energy of the asymptotic adiabatic ground electronic state of Li 2 + Na. The three large black and red dots highlight points for discussion. Fig. 3 Rotationally resolved rate coefficients for the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 = 0, j 0 ) + Na reaction are plotted as a function of collision energy for even exchange symmetry. In panels (a) and (b) the rate coefficients are computed using the coupled two-state diabatic (2 Â 2) and single surface adiabatic (NGP) methods, respectively. The GP which is included in the diabatic 2 Â 2 calculations reverses the quantum interference which leads to a significantly enhanced or suppressed ultracold rate coefficient relative to the NGP calculation which ignores the GP.

Rotationally resolved rate coefficients as a function of collision energy
with the conical intersection shown in Fig. 1 changes the sign of the interference term and hence the nature of the quantum interference from destructive to constructive and vice versa. [14][15][16]21,22,[55][56][57][58][59][60][61][62] Furthermore, due to the unique properties of ultracold collisions, according to Levinson's theorem 63 the scattering phase shifts preferentially approach an integral multiple of p. Thus, the quantum interference often approaches its maximal values effectively turning the reaction on or off (i.e., a quantum switch!). 21,22 To understand this intriguing quantum interference mechanism in more detail, Fig. 6 plots a stereographic projection 17,44 of the hypersphere for the adiabatic ground state PES (V a 1 ) for Li 2 Na at a fixed hyperradius of r = 8a 0 . The north pole (y = 0) is at the origin and the equator (y = p/2) corresponds to the outer black circle. The azimuthal hyperspherical angle (f) is indicated by the dashed blue radial lines. The smaller and larger dashed blue circles indicate the values of y = 101 and y = 531, respectively. The conical intersection is indicated by the red dot and the two symmetric potential wells are clearly visible. These are the same features shown in Fig. 1 plotted in product Jacobi coordinates. The hyperspherical coordinates have the advantage   5 Rotationally resolved rate coefficients for the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 = 3, j 0 = 5) + Na reaction are plotted as a function of collision energy. The red and black curves are the rates computed using the coupled two-state diabatic (2 Â 2) and single surface adiabatic (NGP) methods, respectively. The GP which is included in the diabatic 2 Â 2 calculations gives rise to constructive quantum interference and a significantly enhanced ultracold rate coefficient relative to the NGP calculation which ignores the GP. Fig. 6 The relevant reaction pathways (thick black and red curves) are indicated on top of a contour plot of the Li 2 Na PES (V a 1 ) for a fixed hyperradius r = 8a 0 . A stereographic projection of the hypersphere is plotted with the north pole centered at the origin. The value of the hyperspherical azimuthal angle is indicated by the dashed radial blue lines. The location of the energetically accessible conical intersection is indicated by the red dot and the two symmetric attractive potential wells are clearly visible above and below the horizontal axis to the right of the conical intersection. The contours lie between 15 Â 10 3 and À5 Â 10 3 K inclusive and the dashed contour at 2.23 Â 10 3 K indicates the energy of the Li + LiNa ultracold collision. Quantum interference occurs between the thick black (direct) and red (looping) pathways which can dramatically alter the outcome of the ultracold reaction. The solid (dashed) black and red curves indicate the no-exchange (exchange) contributions. The arrows indicate the direction of the reaction as it proceeds from the reactant Li + LiNa to product Li 2 + Na channels. of showing all channels (reactant and product) simultaneously. The dashed contour at 2228 K indicates the energy of the ultracold reactant Li + LiNa(v = 0, j = 0). It is clear that the CI is energetically accessible even at ultracold collision energies and will be encircled by the quantum mechanical wave function. The reaction pathways are indicated by the thick black and red curves and significant quantum interference can occur between the scattering amplitudes associated with these pathways. For clarity only the pathways from one of the symmetric Li + LiNa reactant channels is shown but both are included in the quantum mechanical calculations which rigorously treat the identical particle exchange symmetry (i.e., the two identical 6 Li nuclei). We note that f -Àf across the horizontal axis corresponds to an exchange of the two Li nuclei. Thus, the pathways from the other LiNa channel can be visualized by reflecting all of the red and black curves in Fig. 6 across the horizontal axis (f = 0). The dashed (solid) black and red curves in Fig. 6 indicate a reactive process with (without) exchange of the identical 6 Li nuclei. There are two primary quantum interference mechanisms occurring simultaneously during the reaction. First, both the exchange (dashed curves) and non-exchange (solid curves) pathways can undergo quantum interference with the contribution from the other LiNa channel. Second, each direct pathway (black) undergoes quantum interference with a looping pathway (red). In the following discussion we will focus on the second mechanism. The first mechanism is always present and is included in the rigorous quantum mechanical treatment of the identical particle exchange symmetry which leads to symmetric and antisymmetric functions with respect to f -Àf. 17 These symmetrized functions correlate asymptotically at large r to the even and odd rotational states of Li 2 (see ref. 22 and 64 for a detailed discussion of the exchange symmetry quantum interference mechanism). The quantum interference between two complex scattering amplitudes (e.g., f direct and f loop ) is given by 57 wheref total ¼ 1 ffiffi ffi 2 pf direct þf loop ; D denotes the relative phase between the two complex scattering amplitudes, and f direct and f loop denote their magnitudes. If the direct and looping scattering amplitudes are of comparable magnitude f direct E f loop = f, then eqn (19) reduces to If D is near an integral multiple of p, D E np where n is an integer, then from eqn (20) we see that for even and odd n, | f total | 2 E 2f 2 and 0, respectively. That is, for even n the reaction is ''on'' and for odd n the reaction is ''off''. As noted above, due to the unique properties of ultracold collisions, 63 the relative phase shifts often approach an integral multiple of p and become effectively quantized. [21][22][23] Also, the direct and looping scattering amplitudes often have similar magnitudes. [21][22][23] Thus, the quantum interference between the black and red pathways shown in Fig. 6 often approaches its maximal values which effectively controls the reaction outcome. 21 It is important to emphasize that this unique quantum interference mechanism (eqn (20)) occurs in ultracold scattering calculations based on a single adiabatic ground electronic state PES, as well as those based on a non-adiabatic coupled diabatic (2 Â 2) electronic state representation. The critical difference is that the GP associated with the CI is included in the 2 Â 2 diabatic calculations 17 but not in the NGP single surface adiabatic calculations. The GP gives rise to an additional p phase shift along the looping (red) pathways so that D -D + p which changes the sign of the cos D term in eqn (20). Thus, when the GP is included the opposite interference occurs (i.e., constructive 2 destructive) and the theoretically predicted reactivity is reversed. That is, for even n the reaction is ''off'' and for odd n the reaction is ''on'' (for additional details see the Appendix). The value of the integer n for a particular scattering process depends upon the details of the PES as discussed below. In practice, it is not computed directly but its parity can be inferred from an analysis of the calculated scattering amplitudes.
The total rate coefficient summed over all final vibrational and rotational states of Li 2 are plotted in Fig. 7. The red and black curves correspond to the 2 Â 2 and adiabatic (NGP) calculations, respectively. The horizontal black dashed line is the ultracold rate coefficient computed using a simple onedimensional universal model based on just a long-range C 6 potential ignoring all reflections. 20 The quantum interference effects tend to wash out in the sum over final states so that the 2 Â 2 and adiabatic total ultracold rate coefficients are similar in magnitude (i.e., K 2Â2 /K adiab E 1.05 at E c = 1.0 nK). Interestingly, both the 2 Â 2 and adiabatic ultracold rate coefficients lie below the universal value (i.e., K 2Â2 /K univ E 0.89 and K adiab /K univ E 0.85 at E c = 1.0 nK). Thus, the smaller 2 Â 2 and adiabatic rates are most likely due to non-reactive (elastic) reflections that are included in the exact quantum mechanical calculations. The sensitivity of the rate coefficients to the accuracy of the PES was also investigated and is discussed in detail in the Appendix. Fig. 7 The total rate coefficient is plotted as function of collision energy for the Li + LiNa(v = 0, j = 0) -Li 2 + Na reaction. The red and black curves correspond to the coupled two-state diabatic (2 Â 2) and single surface adiabatic (NGP) calculations, respectively. The horizontal black dashed line denotes the ultracold rate coefficient computed using a universal model.

Ultracold rate coefficient distributions
All of the rotationally resolved rate coefficients are plotted in Fig. 8 at the ultracold collision energy of 1 nK for each final vibrational product state of Li 2 from v 0 = 0 to 3. The red and black rate coefficients (vertical bars) correspond to the 2 Â 2 and NGP calculations, respectively. Many of the 2 Â 2 rate coefficients are significantly enhanced or suppressed relative to the NGP rate coefficients. As discussed above, this effect is due to the GP which is included in the 2 Â 2 calculations but not in the NGP calculations. The sign change associated with the GP alters the nature of the quantum interference and hence the magnitude of the rate coefficients. For v 0 = 0 (panel a) particularly large GP effects are seen in the product rotational states j 0 = 4, 7, 15, 23, 30, 35-37 and 41 for which the 2 Â 2 rate coefficients are suppressed relative to the NGP ones. In contrast, the product rotational states for j 0 = 24, 34, and 38 show significantly enhanced 2 Â 2 rates coefficients. For v 0 = 1 (panel b) notably suppressed 2 Â 2 rate coefficients are observed for the product rotational states j 0 = 12, 20, 27, and 35 whereas notably enhanced 2 Â 2 rate coefficients occur for j 0 = 14, 21, 24, 26, 30, 32, and 33. For v 0 = 2 (panel c) notably suppressed 2 Â 2 rate coefficients are observed for the product rotational state j 0 = 17 whereas notably enhanced 2 Â 2 rate coefficients occur for j 0 = 1, 11, 24, 27, and 28. Finally, for v 0 = 3 (panel d) notably suppressed 2 Â 2 rate coefficients are observed for the product rotational states j 0 = 4, 8,9,13, and 17 whereas notably enhanced 2 Â 2 rate coefficients occur for j 0 = 3, 5, and 15. In summary, the magnitude of the GP effect on the ultracold rotationally resolved rate coefficients varies significantly across all values of the product ro-vibrational states of Li 2 (v 0 , j 0 ). Fig. 9 plots the normalized distributions s = K/hKi where hKi denotes the average value of the rate coefficients K for a given data set. The probability distributions are computed by binning the K v 0 j 0 into eight equally spaced intervals up to five times the average value. Four normalized data sets are plotted. The red and black data points denote the 2 Â 2 and NGP rate  coefficients, respectively. The circles and squares correspond to the results of even and odd exchange symmetry. The four data sets span all of the vibrational and rotational states shown in Fig. 8. For reference, the Poisson distribution (e Às ) is also plotted (solid black curve). We see that on average all four data sets are consistent with the Poisson distribution. Thus, a statistical analysis of the seemingly random rotational rate coefficient distributions of Fig. 8 provides a unified description of all the results. We note that the Poisson nature of the rotational distributions was also reported previously for the ultracold K + KRb reaction from an NGP calculation on a single adiabatic electronic PES. 12,13 This property appears to be very robust and is independent of the details of the PES and occurs for both the 2 Â 2 and NGP results. For example, in the Appendix Fig. 16

Conclusions
Many ultracold chemical reactions under active experimental investigation, such as Li + LiNa -Li 2 + Na, K + NaK -K 2 + Na, KRb + KRb -K 2 + Rb 2 , and NaRb + NaRb -Na 2 + Rb 2 5,6,9,10 have a barrierless reaction pathway and a deep attractive potential well. In addition, the three atom systems are known to exhibit a C 2v CI between the ground and first excited electronic states in the interaction region. This CI is energetically accessible even for ultracold collisions involving reactant diatomic molecules in their ground ro-vibrational state (e.g., LiNa(v = 0, j = 0)). The four atom systems are also expected to have CIs (possibly multiple CIs) that are energetically accessible. Thus, a standard Born-Oppenheimer calculation based on a single adiabatic electronic PES or even a generalized calculation which includes the GP either by the vector potential approach 14,15 or a double-valued basis approach 56 is not accurate for these systems. A non-adiabatic quantum mechanical dynamics calculation 17 is required that includes both coupled electronic states using full-dimensional accurate ab initio PESs with experimentally based long-range and asymptotic diatomic potentials. In addition, a rigorous quantum mechanical treatment of the exchange symmetry due to the identical nuclei is also required. 17 To the authors' knowledge, the first nonadiabatic calculations of this kind are reported in this work for the ultracold Li + LiNa(v = 0, j = 0) -Li 2 (v 0 , j 0 ) + Na reaction. Indeed, we believe this work significantly advances the state-ofthe-art in the theoretical treatment of experimentally relevant ultracold chemical reactions beyond the Born-Oppenheimer approximation.
During the course of the reaction, two quantum interference mechanisms occur which ultimately govern the reaction outcome. One mechanism is the quantum interference between the two contributions from identical particle exchange symmetry. The second mechanism is the quantum interference between the two reaction pathways (direct and looping) which encircle the CI. Both of these mechanisms are included in all of the quantum mechanical calculations presented in this work (2 Â 2 and NGP). Due to the unique properties of ultracold collisions, the quantum interference often approaches its maximal constructive or destructive values which leads to a significantly enhanced or suppressed rotationally resolved rate coefficient (i.e., the reaction is effectively turned on or off). Furthermore, the GP associated with the CI changes the sign on the interference term between the direct and looping pathways which reverses the nature of the quantum interference. Thus, a nonadiabatic calculation which includes the excited electronic state and its associated GP is crucial for obtaining the correct theoretical prediction of the rate coefficients. A conventional Born-Oppenheimer calculation based on a single adiabatic ground electronic state PES will give the opposite (incorrect) prediction whenever significant quantum interference occurs. Furthermore, a rigorous quantum mechanical treatment of the identical particle exchange symmetry (i.e., the exchange of the two identical 6 Li nuclei) is also required in order to capture the quantum interference associated with exchange symmetry. Several PES scaling studies were undertaken which show that the large quantum interference effects are indeed due to interference between the direct and looping pathways. Furthermore, these PES scaling studies also confirm unambiguously that the GP is responsible for reversing the nature of the interference (i.e., constructive 2 destructive) which leads to the dramatic enhancement or suppression of the 2 Â 2 ultracold rate coefficients relative to the NGP ones. The novel quantum interference mechanism associated with ultracold collisions represents a realization of a molecular quantum switch. The large dynamic range of this quantum switch might be exploited by experimentalists to control the reaction outcome via the application of external fields and/or the selection of a particular initial quantum state. [21][22][23] Interestingly, the interference effects mostly cancel out in the total rate coefficients (both 2 Â 2 and NGP) when summed over all product rotational states. The total rate coefficients also exhibit non-universal behavior.
The rotationally resolved rate coefficient distributions are shown to exhibit Poisson behavior. The S matrix for open chaotic quantum systems obeys the statistics of unitary symmetric random matrices, one of which is the Poisson law behavior of the squares of off-diagonal matrix elements. 65,66 Since state-to-state rates are directly proportional to the square of the corresponding S matrix element, this Poisson law behavior follows directly from the underlying classically chaotic motion of the reaction. 67 Chaotic classical trajectories are extremely complicated and tangled for reactions with longlived intermediate complexes, as such these results show that the ultracold LiNa + Li reaction proceeds via complex formation. Such intermediate complexes can be observed experimentally using a combination of mass spectrometry and velocity map imaging, as was recently demonstrated for the ultracold KRb + KRb -K 2 + Rb 2 reaction. 10 As shown explicitly in this work for the first time, the Poisson nature of these rotational distributions is robust to variations in the PES, occurs for different chemical systems (i.e., both light Li 2 Na and heavy K 2 Rb 12 ) and theoretical methods (i.e., both nonadiabatic (2 Â 2) and adiabatic (NGP)). The robust and universal nature of the Poisson ultracold rotational rate coefficient distributions makes this property an ideal experimental observable. Considering the sensitivity of individual rotationally resolved rate coefficients to fine details of the PES a direct quantitative comparison of these with experimental measurements will be challenging and will likely require improvements to the PES.
In contrast to standard statistical theories which predict uniform product distributions for many complex-forming reactions, 68 we find that at ultracold collision temperatures the product rotational distributions are significantly nonuniform (i.e., non-statistical). This is due in part to the unique properties of ultracold collisions which involve only a singlepartial wave. At higher (thermal) energies many partial waves contribute which tend to wash-out the variations in the product distributions resulting in a uniform distribution. Statistical theories also ignore the detailed dynamics of the interaction region and therefore do not explicitly include the short-range quantum dynamics, non-adiabatic, and interference effects which are significantly amplified at ultracold temperatures and give rise to the large variations in the rotational distributions. Nevertheless, statistical theories can often provide a reasonable estimate of the total ultracold rate coefficients due to the large density of states associated with complex forming reactions with a deep attractive potential well and heavy nuclei. 69 In an exact quantum mechanical calculation, the fluctuations in the rotational distributions average out in the sum over the large number of product states so that the total ultracold rate coefficient is often similar to the statistical one. However, we expect that the underlying distributions are all Poissonian and also occur in ultracold diatom-diatom reactions such as KRb + KRb and NaRb + NaRb. For additional analysis comparing statistical theories to exact quantum mechanical calculations, we refer the reader to recent work on the ultracold Li + LiYb reaction. 11 The unique properties of ultracold collisions are still largely unexplored. We hope that the detailed theoretical results presented in this work will help stimulate new experimental and theoretical studies into the intriguing ultracold energy regime. For example, the extreme non-linear sensitivity of the quantum interference effects on the depth of the PES might be used for very accurate experimental probing and benchmarking of the molecular interactions, beyond anything used in the field at the present time. One could also imagine exploiting this ultra sensitivity in the search of new physics such as tightening the experimental bounds on possible variations in the physical constants (e.g., the proton to electron mass ratio). Clearly, ultracold molecules continue to show exceptional promise for future technological applications in quantum control, sensing and precision measurements.

Conflicts of interest
There are no conflicts to declare.
A.1 Electronic ground and first excited state adiabatic potential energy surfaces dashed blue circle about the origin indicates y = 101 (531). The barrierless reaction pathway is shown by the thick black curves with arrows indicating the reaction direction. The reaction proceeds along a collinear approach NaLi-Li (panel a) and then enters the attractive potential wells in the interaction region (panel b). The identical particle exchange symmetry (i.e., f -Àf) associated with the two identical Li nuclei is rigorously accounted for in the quantum mechanical calculations. The black dots denote the location of the minimum energy of the PES for the specified fixed value of hyperradius: r = 14.5a 0 (panel a) and r = 8.0a 0 (panel b). The red dot in panel (b) denotes the location of the C 2v conical intersection. The dashed contour line at À2.23 Â 10 3 K corresponds to the energy of the reactant Li + LiNa(v = 0, j = 0) channel. The conical intersection lies at À3.14 Â 10 3 K and is therefore energetically accessible even at ultracold collision energies. All energies are relative to the asymptotic Li 2 + Na adiabatic ground electronic state energy minimum. We note that the sharp feature in the thick black curve (minimum energy curve) near r = 9.5a 0 in Fig. 2 is due to a submerged barrier along the reaction path shown in Fig. 10(b). Fig. 11 is the same as Fig. 10 except that the adiabatic excited electronic state PES (V a 2 ) is plotted at r = 8.0a 0 . The contours exhibit the typical inverted cone shape associated with a conical intersection. The excited state becomes degenerate with ground electronic state at the conical intersection (i.e., the red dots in Fig. 10 and 11). As also shown in Fig. 10, the dashed contour at À2.23 Â 10 3 K lies above the energy of the conical intersection. Thus, the Born-Oppenheimer approximation breaks down necessitating a non-adiabatic coupled two electronic state quantum mechanical treatment.
A.2 Reaction rate coefficients and potential scaling studies Fig. 12 plots the ultracold total rate coefficient for the Li + LiNa(v = 0, j = 0) -Li 2 + Na reaction as a function of the three-body potential scaling factor (l). The ab initio computed three-body contribution to the potential energy surface (PES) was scaled by l but the long-range and pairwise two-body potentials were left unchanged. This ensures that the asymptotic energies and long-range interactions are unchanged and that only the effective depth of the Li 2 Na PES is altered. Specifically, the 3-body components in the diabatic potentials  Fig. 10 except that the contour plot of the adiabatic excited electronic state of Li 2 Na (V a 2 ) is plotted in hyperspherical coordinates for a fixed hyperradius r = 8a 0 . The contours lie between 15 Â 10 3 and À2 Â 10 3 K inclusive and the dashed contour at 2.23 Â 10 3 K indicates the energy of the Li + LiNa ultracold collision. As in Fig. 10, the red dot indicates the location of the energetically accessible conical intersection.  2 -Na (i.e., C 2v ) geometries at À3.14 Â 10 3 K. All energies are relative to the minimum energy of the asymptotic adiabatic ground electronic state of Li 2 + Na. Fig. 12 The total rate coefficients for the Li + LiNa(v = 0, j = 0) -Li 2 + Na reaction are plotted as a function of the 3-body potential scaling parameter at the ultracold collision energy of 1 nK. The red and black data correspond to the coupled two-state diabatic (2 Â 2) and single surface adiabatic (NGP) calculations, respectively. The large circular data points are the total rate coefficients computed by adding the statistically weighted rate coefficients for even exchange symmetry (squares) and odd exchange symmetry (triangles). To guide the eye, the data points for the total, even, and odd rate coefficients are connected by thick lines, thin lines and dashed lines, respectively. V 11 , V 22 , and V 12 (see eqn (9)) were scaled according to the following equation where the scale factor is given by l = 1 + l 0 f switch with a standard switching function f switch = 1/2[1 À tanh(a(r À r 0 ))]. The parameters a = 0.75 and r 0 = 16a 0 were chosen to minimize the effects of the scaling on the long-range component of the PES. The parameter l 0 was chosen to reflect the estimated uncertainty (AE3%) in the ab initio computed 3-body interaction PES. Thus, the scaling studies were chosen to span a range of 25 equally spaced scaling factors between l 0 = AE3% (i.e., l = 0.97 to l = 1.03). A total of 100 = 4 Â 25 separate l calculations (which include l = 1 no scaling) were performed for the 2 Â 2 and NGP methods and each exchange symmetry even and odd. The 2 Â 2 and NGP results are presented in red and black, respectively. The total, even, and odd rate coefficients are plotted with solid circles, squares, and triangles, respectively. The NGP total rate coefficient oscillates between 2.16 Â 10 À10 and 3.74 Â 10 À10 cm 3 s À1 (i.e., by À20% and +38% relative to the unscaled NGP rate coefficient). The 2 Â 2 total rate coefficient oscillates between 2.46 Â 10 À10 and 4.49 Â 10 À10 cm 3 s À1 (i.e., by À13% and +58% relative to the unscaled 2 Â 2 rate coefficient). While computationally demanding, the scaling studies show that the 2 Â 2 and NGP total rate coefficients are robust to variations in the PES. The total rate coefficients of even exchange symmetry are on average two times larger than the odd ones due to the statistical weights of 2/3 (even) versus 1/3 (odd). The computational time for each of the 100 l calculations is primarily dominated by the surface function solutions (i.e., the diagonalization of eqn (18) at each r x ) and the radial logderivative propagation (i.e., the solution of the radial eqn (16) from r i to r f ). All steps of the calculations are parallelized and can be efficiently distributed across a large number of CPUs simultaneously. For comparison, each set of 2 Â 2 surface function calculations was performed using 288 CPUs and required t sf 2Â2 = 6.2 h whereas the NGP surface function calculations on the same number of CPUs required t sf NGP = 0.6 h. Each CPU is an Intel Haswell (Xeon E5-2650 v3) running at 2.3 GHz. The 2 Â 2 and NGP log-derivative propagations used 16 CPUs (one for each collision energy) and required t log 2Â2 = 8.8 and t log NGP = 3.5 h, respectively. Thus, the total wall-clock time for one l calculation was t total 2Â2 = t sf 2Â2 + t log 2Â2 = 15.0 h and t total NGP = t sf NGP + t log NGP = 4.1 h. The non-adiabatic 2 Â 2 calculations require about t total 2Â2 /t total NGP = 3.6 times more CPU time than the NGP ones. The 100 l calculations were separately distributed across approximately 20 compute nodes (each node consists of 32 CPUs and 128 Gb memory) and were completed over the course of approximately three weeks. Fig. 13 plots the ultracold rotationally resolved rate coefficient for the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 = 3, j 0 = 5) + Na reaction as a function of the three-body potential scaling factor (l). The 2 Â 2 and NGP results are presented in red and black, respectively. As the scaling factor increases or decreases away from l = 1 (i.e., no scaling), dramatic and often sudden changes in both the 2 Â 2 and NGP rate coefficients are observed.
These large changes are due to the change in nature of the quantum interference from constructive to destructive and vice versa. The quantum interference occurs between the direct and looping contributions to the total scattering amplitude. The direct pathway proceeds straight from reactants to products through the Li 2 Na potential well whereas the looping pathway ''loops around'' and encircles the conical intersection (see Fig. 6 in the paper).
Let f direct = f direct e iddirect and f loop = f loop e idloop denote the complex scattering amplitudes associated with the direct and looping pathways, respectively. The f direct and f loop are the real valued magnitudes and the d direct and d loop are the real valued phases. The total scattering amplitude is the sum of the two: The cross sections and rate coefficients are computed from the modulus of the total scattering amplitude given by eqn (19) where the relative phase D = d loop À d direct . The third term on the right hand side of eqn (19) is the interference term. If the looping and direct scattering amplitudes are comparable in magnitude: f direct E f loop = f, then we obtain eqn (20). At ultracold collision energies the scattering phase shifts d loop and d direct have a propensity to approach an integral multiple of p (i.e., Levinson's theorem): d direct E pn direct and d loop E pn loop where the n direct and n loop are integers. Thus the relative phase also approaches an integral multiple of p: D E p(n loop À n direct ) E np where the integer n = n loop À n direct . This implies that the interference term cos D in eqn (20) approaches: cos D E 1 or cos D E 0 for even and odd values of n, respectively. From Levinson's theorem we know that the integers n direct and n loop denote the number of bound states supported along the direct and looping reaction pathways, respectively. Thus, the integer n is the relative number of bound states between these two pathways which in general is either even or odd. In practice, we never need to explicitly compute these integers (or D) but can infer the even or odd nature of n from the scattering amplitudes. Eqn (19) and (20) are valid for both the 2 Â 2 and NGP calculations. However, the 2 Â 2 Substituting the above expressions for the looping and direct phase shifts for d 2Â2 and d NGP into the D in eqn (20) (and ignoring e), we obtain where the + sign corresponds to NGP and the À sign corresponds to 2 Â 2 (or GP). If n is an even integer, then eqn (A.2) gives If n is an odd integer, then eqn (A.2) gives In the ultracold limit, eqn (A.3) and (A.4) show that the quantum interference can approach its maximal values of 2f 2 (constructive) or zero (destructive). The constructive or destructive nature of the interference is always opposite for the NGP and 2 Â 2 calculations and can reverse if the integer n changes from being even to odd or vice versa. We note that the maximal values given by eqn (A.3) and (A.4) are not fully realized in practice since the direct and looping scattering amplitudes are not exactly equal in magnitude and their scattering phase shifts are not exactly an integer multiple of p. However, as discussed below, these limiting cases are useful for understanding and interpreting the large differences observed between the 2 Â 2 and NGP results.
Armed with eqn (A.3) and (A.4), we can now understand and interpret the changes in the rate coefficients plotted in Fig. 13. First consider the scaling range between l = 0.975 and 1.0025. In this region, the 2 Â 2 rate coefficient is significantly enhanced relative to the NGP one (by 50 times for l = 1). This enhancement is due to constructive quantum interference as shown in eqn (A.4) which implies n must be odd in this region. Fig. 14 plots the ratio of the magnitudes of the looping and direct scattering amplitudes averaged over the scattering angle as a function of the PES scaling factor. We see that for 0.98 r l r 1.0025 the ratio is near unity which ensures maximal quantum interference. Fig. 15 plots cos D averaged over the scattering angle as a function of the PES scaling factor. In the region 0.975 r l r 1.0025 we see that hcos Di lies near À1.0 (except for a few values between 0.9825 and 0.99). The negative cos D is consistent with an odd value of n. As l is decreased below 0.98, the ratio in Fig. 14 decreases rapidly to below 10 À1 . Thus, the direct scattering amplitude dominates and quantum interference effects become small. In this case, the 2 Â 2 and NGP rate coefficients approach each other. Also in this region, the cos D term reverses sign for l r 0.9725 so that the nature of the quantum interference is reversed and now eqn (A.3) is relevant (i.e., n becomes even). Indeed, in this region the NGP rate coefficient is now larger than the 2 Â 2 one. The change in n from being odd to even (or vice versa) is due to the sudden change in the relative number of bound states between the direct and looping pathways. Presumably this change is due to a bound (continuum) state leaving (entering) the well and becoming a continuum (bound) state as the well depth is decreased (increased) by the l scaling. This explains the sudden change in sign of cos D between l = 0.9725 and 0.975. Now consider the region 1.0025 o l r 1.0125. In this region the ratios in Fig. 14 lie between 1 and 10 and the cos D are positive and approach + 1 for l = 1.01. Again the positive cos D implies an even n (i.e., eqn (A.3) is relevant) and the NGP rate coefficient Fig. 14 The ratios of the modulus of the looping (f loop ) and direct (f direct ) scattering amplitudes averaged over the scattering angle are plotted as a function of the 3-body potential scaling parameter at the ultracold collision energy of 1 nK. The ratios correspond to the rotationally resolved rate coefficients plotted in Fig. 13. is larger than the 2 Â 2 one. The quantum interference and hence the difference between the rate coefficients is largest in this region for l = 1.01 where cos D E 1. Between l = 1.0125 and 1.015, the cos D term changes sign and becomes negative again reversing the nature of the quantum interference (i.e., n becomes odd again). The 2 Â 2 rate coefficients are now larger than the NGP ones in Fig. 13. As l is increased further, the cos D term oscillates back to being positive again with a peak value at l = 1.0275. In the region of large positive l, the ratio in Fig. 14 continues to decrease rapidly which indicates that the direct pathway is dominant. Thus, the quantum interference decreases considerably and the 2 Â 2 and NGP rate coefficients become comparable in magnitude as seen in Fig. 13. We note that a similar analysis was done for all of the other rotationally resolved rate coefficients and they are all consistent with eqn (A.3) and (A.4). The PES scaling studies show unambiguously that the significant enhancement or suppression of the 2 Â 2 rate coefficients relative to the NGP ones is due to the GP. Fig. 16 and 17 plot the ultracold probability distributions of the rotationally resolved rate coefficients for the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 , j 0 ) + Na reaction for the NGP and 2 Â 2 calculations, respectively. The normalized distributions s = K/hKi (where hKi denotes the average value of the rate coefficients K for a given data set) are computed by binning the rate coefficient K v 0 j 0 into equally space bins up to 14 times the average value. In Fig. 16 and 17 the even and odd curves are plotted in dark and light blue (red), respectively. There are 25 curves for each exchange symmetry in each plot which correspond to the 25 values of the scaling parameter l. On average, all of the curves are consistent with the Poisson distribution e Às (black curve). Fig. 16 The probability distributions for all of the rotationally resolved rate coefficients for the Li + LiNa(v = 0, j = 0) -Li 2 (v 0 , j 0 ) + Na reaction are plotted at the ultracold collision energy of 1.0 nK. The distributions correspond to the single surface adiabatic (NGP) calculations. The results for even and odd exchange symmetry are plotted in dark and light blue, respectively. Each curve plots the distribution for a different 3-body scaling parameter (see Fig. 13). The solid black curve is the Poisson distribution.