Non-adiabatic quantum interference and complex formation in ultracold collisions of Rb with KRb

Abstract

Ultracold elastic collisions of ⁸⁷Rb with ⁴⁰K⁸⁷Rb in its ground vibrational and rotational state are investigated using a first principles based theoretical methodology. Full-dimensional ab initio computed potential energy surfaces are reported that include the two lowest-lying electronic states, their conical intersection, non-adiabatic couplings and an accurate long-range behavior. A numerically exact time-independent quantum dynamics method in hyperspherical coordinates is used to compute the elastic scattering cross sections, rate coefficients and collision lifetime spectrum. The quantum scattering calculations include all degrees of freedom and treat both electronic states and their non-adiabatic couplings using a two-state diabatic representation. The theoretically computed elastic rate coefficient is in good agreement with the recently reported experimental value. Significant non-adiabatic quantum interference effects are shown to originate from the unique properties of ultracold collisions and the geometric phase associated with the conical intersection. A high-resolution collision energy grid is used to investigate the origin of the experimentally reported long-lived 3-body collision complexes.

---

1 Introduction

The control and manipulation of ultracold molecules is crucial to the development and application of new technological applications in quantum information, sensing and precision measurement.^1–10 A key requirement in achieving these novel technologies is the formation and trapping of stable ensembles of ultracold molecules with high densities. Unfortunately, the densities are often limited by collisional losses which must be understood and mitigated. Most notable is the recent experimental evidence for the formation of very long-lived complexes during bi-molecular (e.g., NaK + NaK, NaRb + NaRb)^11,12 and atom-molecule collisions (e.g., Rb + KRb).^13–15 These long-lived collision complexes can lead to significant trap losses due to unwanted optical excitation of the complex within the magneto-optical trap.¹⁴ Furthermore, the experimental lifetimes of these complexes are 10³ to 10⁵ times larger than statistical theory estimates.^16–21 Clearly, new theoretical insight is needed to elucidate the origin of these mysterious long-lived collision complexes and to ultimately provide trap loss mitigation strategies.²² A primary motivation of the present study is to investigate the role of an excited electronic state in regards to complex formation. Specifically, could the mysterious long-lived complexes be due to resonances associated with localized vibrational bound states lying within the inverted cone of an excited electronic state? In addition, based on our previous ultracold studies, we expect that quantum interference effects and the geometric phase will play an important role in governing the overall collision outcome.^23–27

Ideally, a first principles based quantum dynamics calculation including all degrees of freedom, electronic and spin interactions as well as external fields is desirable. Unfortunately, such a complete treatment is beyond current computational capabilities even for atom–diatom collisions. For example, recent theoretical treatments that include both spin and external field interactions have been limited to a reduced dimensional (rigid-rotor) approximation.^15,28 On the other hand, full-dimensional theoretical treatments include excited electronic state interactions but ignore spin and external fields.^25,27 Both of these theoretical approaches are providing novel insight into the importance of the different kinds of interactions present in the ultracold collision regime. In the present study, we investigate ultracold collisions of ⁸⁷Rb with ⁴⁰K⁸⁷Rb in its ground vibrational (v = 0) and rotational (j = 0) state using a numerically exact full-dimensional time-independent quantum scattering method in hyperspherical coordinates.^29–34 This methodology enables the treatment of both the ground and first excited electronic states.³⁴ As for Li + LiNa²⁵ and K + KRb,²⁷ the excited electronic state in Rb + KRb is energetically accessible within the interaction region even at ultracold collision temperatures and must be included in the quantum scattering calculations. The excited electronic state exhibits a conical intersection (degeneracy) with the ground electronic state near the bottom of the attractive potential well. A two electronic state diabatic representation is used which includes the geometric phase and other non-adiabatic couplings.³⁴ The diabatic potential energy matrix is computed from first principles using state-of-the-art electronic structure methods and also includes an accurate treatment of the long range interactions.

In summary, the non-adiabatic scattering methodology is briefly outlined followed by a detailed description of the ab initio two electronic state diabatic potential energy matrix calculations (this includes two diagonal and one off-diagonal three-dimensional potential energy surfaces). The quantum scattering methodology is reviewed followed by computational details and the various parameters used in the final production runs are given. Important features of the potential energy surface are then discussed and the scattering results are presented which include the elastic cross sections, rate coefficients and resonance (complex lifetime) spectrum followed by discussion and conclusions.

2 First principles non-adiabatic methodology

An accurate theoretical treatment of ultracold collisions of Rb with KRb requires a first principles based non-adiabatic quantum dynamics calculation including both the ground and first excited electronic states of KRb₂. For computational reasons, it is advantageous to use the diabatic electronic representation in lieu of the adiabatic one.³⁴ Thus, the total molecular wavefunction |Ψ〉 is expanded as , where |φ^d_i〉 denotes the diabatic electronic basis functions and the expansion coefficients [small psi, Greek, tilde]_i are the nuclear motion wave functions. The relevant diabatic Schrödinger equation for the nuclear motion is given bywhere the first term in brackets is the diagonal kinetic energy operator for the nuclear motion with matrix elements T̂ = −ħ²∇²/(2μ), μ is the three-body reduced mass and ∇ denotes derivatives with respect to the six-dimensional set of coordinates x = (x, x̂) relative to the center of mass of the molecule. The x denotes the three internal coordinates (e.g., internuclear distances) and x̂ denotes the three Euler angles (α, β, γ) that specify the orientation of the molecule's body frame (BF) relative to the space fixed frame. The second term is the PEM (potential energy matrix) which couples the two diabatic electronic states. Each matrix element (i.e., the V^d₁₁, V^d₂₂ and V^d₁₂ = V^d₂₁) of the PEM is a three-dimensional PES (potential energy surface) and is a function of the x only. For reference, a conventional Born–Oppenheimer quantum dynamics calculation solves the adiabatic single surface Schrödinger equation

[T̂ + V^a₁(x)]ψ₁(x) = Eψ₁(x). (2)

where V^a₁ denotes the adiabatic ground electronic state PES (which is obtained by diagonalizing the 2 × 2 diabatic PEM, see Section 2.1). A first-principles solution of eqn (1) requires two primary steps: (1) an ab initio based calculation of the PEM and (2) a fully coupled quantum dynamics solution of the nuclear motion using the PEM computed in step (1). These two steps are discussed in Sections 2.1 and 2.2, respectively.

2.1 Ab initio potential energy surfaces

We have performed ab initio electronic structure calculations to determine the relevant adiabatic and diabatic electronic potential energy surfaces of KRb₂. We require both the ground and first-excited adiabatic electronic states with total electron spin 1/2. The potentials of the two electronic states are degenerate when all three atoms are well separated and have a seam of conical intersections (CI) in isosceles triangle geometries. The CI leads to significant non-adiabatic quantum interference effects²⁵ and possibly long-lived collision complexes, which may explain the long lifetimes observed in ref. 14.

Our electronic structure calculations of the KRb₂ complex have been carried out with the MOLPRO software package.³⁵ Here, the chemically inert core electrons of K and Rb are replaced by the Stuttgart–Cologne energy-consistent relativistic pseudopotentials ECP18SDF and ECP36SDF, respectively, thereby reducing the number of “active” electrons to three valence electrons. Core–valence correlation and polarization effects are included by an electron orbital angular momentum or l-independent core-polarization potential (CPP) with Müller–Meyer damping functions³⁶ using literature values for dipole-polarizabilities as well as optimized cutoff radii³⁷ for each atomic core. For the treatment of valence electrons, uncontracted sp basis sets augmented by additional s, p, d, and f polarization functions have been adopted.³⁸ The motion of the three valence electrons is assumed to be non relativistic.

Multi-configurational self-consistent-field (MCSCF) calculations^39,40 generate reference configurations, which are then used in multi-reference configuration–interaction (MRCI) calculations^41,42 to determine the two lowest adiabatic potential energy surfaces V^a₁(x) and V^a₂(x). The potentials only depend on a limited number of coordinates of the six-dimensional vector x. This limited set of coordinates is easiest described in terms of Jacobi coordinate r and R, where three-dimensional vector r is the orientation of and separation between the ⁴⁰K atom and one of the ⁸⁷Rb, while three-dimensional vector R represents the orientation of and separation between the the second ⁸⁷Rb atom and the center of mass of isotopologue ⁴⁰K⁸⁷Rb. The potentials only depend on separations r and R and angle θ between r and R. These Jacobi coordinates are the most convenient for the limited phase-space sampled by the non-reactive ultracold Rb + KRb collision.

MRCI simulations have been performed on a non-equidistant grid of angles θ given by the 13-point Gauss–Legendre quadrature on the interval [0, π]. For the radial degrees of freedom, dense grids are used that span the chemically relevant region and extend into the asymptotic domain. Specifically, r is varied from 3.2a₀ to 10a₀ with a uniform step of 0.2a₀, while R is varied from 10a₀ to 19a₀ using a larger increment of 0.3a₀, respectively. The resulting three-dimensional grid has over thirteen thousand nuclear geometries.

At each geometry, the non-adiabatic coupling terms have been numerically evaluated using the finite-difference DDR procedure⁴³ implemented in MOLPRO. These couplings are then integrated to obtain the mixing angle β(r, R, θ), defining the adiabatic-to-diabatic (ADT) transformation, on our grid of geometries. In fact, we have V^d₁₁ = V^a₁ cos²(β) + V^a₂ sin²(β), V^d₂₂ = V^a₂ cos²(β) + V^a₁ sin²(β) and V^d₁₂ = V^d₂₁ = (V^a₂ − V^a₁)cos(β)sin(β) dropping the arguments of β and the potential matrix elements for clarity.

In order to obtain continuous and differentiable potential energy surfaces suitable for quantum scattering calculations, the discrete sets of diabatic potential energies and coupling elements are interpolated by using a pairwise-potential model to extract an intrinsic diabatic three-body potential matrix U^d_ij(x) from V^d_ij(x) and a reproducing-kernel Hilbert-space procedure^44,45 to interpolate U^d_ij(x). For these two steps, we find it convenient to change coordinates from (r, R, θ) to the equivalent coordinates (r_KRba, r_KRbb, θ_RbKRb), where r_KRbk are bond lengths with k = a and b, representing the two rubidium atoms in the tri-atomic molecule, and angle θ_RbKRb is given by ∠Rb_a–K–Rb_b.

As in our previous research,^25,27 the 2 × 2 pairwise potential matrix is constructed from ab initio electron-spin singlet (X¹Σ⁺) and triplet (a³Σ⁺) diatomic potentials of KRb and Rb₂ obtained with MCSCF and MRCI calculations using the same bases as the simulations for the trimer. The singlet and triplet potentials approach the same dispersion potential −C₆/x⁶ − C₈/x⁸ − C₁₀/x¹⁰ for large separations x between the two atoms in the diatom.^46,47 The coefficients C_2n for n = 3, 4, and 5 are positive.

The intrinsic three-body potential matrix elements U^d_ij(r_KRba, r_KRbb, θ_RbKRb) are interpolated by first expanding the angular dependence using Legendre polynomials P_n(cos θ_RbKRb) with n = 0, 1, 2,… for each pair of bond lengths. The resulting intrinsic coefficients A_ij,n(r_KRba, r_KRbb) are then interpolated using the reproducing kernel Hilbert space (RKHS) method with two radial kernels. The matrix elements U^d_ij(r_KRba, r_KRbb, θ_RbKRb) are forced to smoothly and rapidly go to zero for large bond lengths r_KRba, r_KRbb, and r_RbaRbb (computed with the help of angle θ_RbKRb) by a switching function constructed from a product of three corresponding sigmoidal functions, each centered at separation r_f = 28a₀.

Fig. 1(a) shows the two interpolated adiabatic potential surfaces of KRb₂ as functions of bond separations r_KRba and r_KRbb with a bond angle Rb_a–K–Rb_b fixed at 63° degrees. The surfaces are symmetric with respect to the line r_KRba = r_KRbb. In the graph the ultracold Rb + KRb collision has an energy of approximately hc × −4250 cm⁻¹ and “enters” the energetically lowest surface as a superposition of states with either large r_KRba or r_KRbb and the other bond length equal to the equilibrium X¹Σ⁺ KRb separation of approximately 7.7a₀. This superposition exhibits either even or odd symmetry with respect to the exchange of the two identical Rb atoms.

Fig. 1 (a) Ground (red surface) and first-excited (green surface) spin-doublet adiabatic electronic potentials of KRb₂ as functions of the separations or bond lengths of K with either Rb atom, where the bond angle Rb–K–Rb is fixed at θ_RbKRb = 63° degrees. The red and green curves are contour lines with the same potential energy. The zero energy corresponds to the asymptotic limit where all three atoms are infinitely far apart. (b) Parametric graph of the potential energy along the seam of the conical intersections as a function of the bond angle θ_RbKRb and the bond length between the K and Rb atom for isosceles geometries.

Fig. 1(a) also shows a conical intersection where the two adiabatic potential surfaces touch. This occurs at r_KRba = r_KRbb = 8.4a₀ and is energetically accessible for our ultracold collision. In fact, as shown in Fig. 1(b) this conical intersection is part of a one-dimensional curve or seam of conical intersections that is contained within the two-dimensional plane of isosceles geometries with r_KRba = r_KRbb. The conical intersection in the figure has the lowest potential energy of all points on the seam at θ_RbKRb = 63.7° (see Table 1). Not shown in the figure is that the equilibrium geometry of the ground-state adiabatic potential V^a₁(r_KRba, r_KRbb, θ_RbKRb) occurs at an isosceles geometry with the two KRb bond lengths equal to 8.04a₀ and a Rb–K–Rb angle of 80.9° (see Table 1). A subroutine which evaluates the two adiabatic and diabatic electronic potential energy surfaces as a function of the three internuclear distances is available for download in the SI.

Table 1 Energies and geometries for various points of interest on the KRb₂ adiabatic electronic potential energy surfaces. All energies are relative to the minimum energy of the asymptotic ground electronic state of Rb + KRb. All distances are in units of a₀ and angles are in units of degrees. The hyperspherical coordinates are (ρ, θ, ϕ), the three internuclear distances are (r_Rb2, r_KRb, r_KRb) and the Jacobi coordinates are (R, r_Rb2,Θ)

Location	Energy (eV)	(K)	ρ	θ	ϕ	rRb2	rKRb	rKRb	R	Θ
V1 (C2v min)	−0.2358	−2736.3	12.580	36.18	180.0	10.43	8.043	8.043	6.123	90.0
V2 (C2v min)	−0.1682	−1951.9	11.506	19.46	180.0	8.736	8.374	8.374	7.114	90.0
V1 (bent min)	−0.2318	−2689.9	11.662	18.11	127.2	8.358	9.907	8.030	7.991	104.6
—	—	—	—	—	−127.2	8.358	8.030	9.907	7.991	75.4
KRb (min)	0.0	0.0					7.7034
Rb2 (min)	0.024594	285.40				7.9835
Rb + KRb(0,0)	0.004590	53.26
K + Rb2 (0,0)	0.028084	325.90
ΔE	0.023494	272.64

---

2.2 Non-adiabatic quantum dynamics

The numerically exact quantum mechanical solution of the 2 × 2 diabatic Schrödinger eqn (1) is based on a time-independent coupled-channel (CC) approach using hyperspherical coordinates.^29,31,34 This methodology has been previously documented in detail so only a brief summary will be given here. The kinetic energy operator in hyperspherical coordinates can be expressed aswhere the first term is the radial part and the second term is the angular part. The operator [capital Lambda, Greek, circumflex] is the grand angular momentum operator (see ref. 31 and 34 for the explicit expression). This operator contains derivatives with respect to ϕ and the polar angle [small theta, Greek, tilde] = π − 2θ. It also contains the operators J_x, J_y, J_z which 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 KRb molecule is in its ground rotational state j = 0, this corresponds to total angular momentum J = l + j = 0.

We expand the total diabatic nuclear motion wave function in eqn (1) in terms of the angular functions (Φ) for a given sector ρ = ρ_ξ as

where i and t denote the CC indices, p denotes inversion parity, q denotes the identical particle exchange symmetry which is relevant for triatomic molecules with identical atoms, and w ≡ ([small theta, Greek, tilde], ϕ, α, β, γ) denotes the angular coordinates of the 5D hypersphere. The radial coefficients ζ^Jpq_it are computed numerically from the propagation of the CC radial equationswhere the potential coupling matrix is given byThe Φ^JMpq_t are the 5D hyperspherical functions which are solutions to the angular equation at sector ρ = ρ_ξThe V^d in eqn (6) and (7) is the 2 × 2 diabatic PEM given in eqn (1). Once the angular solutions to eqn (7) are computed, they can be used to compute the overlap matrices and potential coupling matrices defined in eqn (6) at each sector. These are then used to solve the CC radial eqn (5) using Johnson's log-derivative propagator.^48,49

As ρ becomes large, the different arrangement channels (i.e., Rb + KRb and K + Rb₂) become isolated and well separated by repulsive energy barriers. Thus, at some optimal large value of ρ = ρ_m (determined by convergence studies), the radial solution is transformed into a Delves coordinate basis using the appropriate overlap matrices between the APH and Delves functions both computed at the matching radius ρ_m. The propagation is then continued to a very large final value of ρ = ρ_f using an equivalent set of Delves CC radial equations.⁵⁰ The Delves hyperspherical coordinates are much more efficient and accurate for large ρ since the primitive basis functions are localized within each arrangement channel. At the final value of ρ_f, the asymptotic scattering boundary conditions are applied and the full state-to-state scattering (S) matrix is computed.^29,31,34 The relevant cross sections σ_fi and rate coefficients K_fi = vσ_fi (where v is the relative collision velocity) can then be computed from the S matrix using standard expressions.²⁹ The indices f, i denote the collective final and initial quantum numbers of the diatomic molecules (i.e., f = (τ′, v′, j′) and i = (τ, v, j)) where τ′ and τ denote the diatomic arrangement channel (Rb₂ or KRb), respectively. In the present study, we consider ultracold Rb + KRb elastic collisions for which the K + Rb₂ channel is closed (i.e., non-reactive so that τ = τ′ = KRb and also v′ = v = 0 and j′ = j = 0).

The 5D surface functions in eqn (7) are expanded in terms of Jacobi polynomials in [small theta, Greek, tilde] labeled by l, complex exponential functions in ϕ labeled by m and symmetrized Wigner D functions in the Euler angles.^31,34 The number of basis functions in [small theta, Greek, tilde] and ϕ are specified by the parameters l_max and m_max, respectively. In addition, an energy cut-off is specified at each sector ρ_ξ which is used to significantly truncate the size of the total direct product basis in [small theta, Greek, tilde] and ϕ (via the Sequential Diagonalization Truncation technique⁵¹). Convergence studies are required to determine the optimal values of l_max, m_max and E_cut for each value of ρ_ξ. The sparse Hamiltonian matrix in eqn (7) is diagonalized numerically using the parallelized version of the ARPACK library.⁵² For the 2 × 2 diabatic case (eqn (1)) this matrix is complex Hermitian whereas for the adiabatic single surface case (eqn (2)) this matrix is real symmetric. Thus for the 2 × 2 diabatic case, the radial solutions of the CC eqn (5) and the associated potential coupling and overlap matrices are complex valued.

For the alkali systems (e.g., Rb + KRb, K + KRb, Li + LiNa and Li + Li₂) the excited electronic state is energetically accessible only in the interaction region (small ρ). Asymptotically at large ρ this excited electronic state lies very high in energy and is closed. Thus, the calculations can be greatly simplified by including the excited electronic state and solving the 2 × 2 diabatic eqn (1) only in the interaction region. At some intermediate value of ρ = ρ_DAT where the electronic state becomes sufficiently closed, the 2 × 2 diabatic radial solutions are projected onto the ground electronic state adiabatic solutions obtained from solving eqn (2). This projection is performed using the diabatic-to-adiabatic transformation (DAT) matrix which is the inverse of the 2 × 2 ADT matrix.³⁴ In practice, to preserve unitarity this transformation must be performed using many sub-steps.³⁴ Once the 2 × 2 radial solutions have been projected onto the adiabatic ones at ρ = ρ_DAT, the radial propagation is continued using the CC equations for the adiabatic ground electronic state. As discussed above, this radial propagation continues until ρ = ρ_m where the radial solutions are transformed to the Delves coordinates and the propagation continued to the final asymptotic value ρ = ρ_f. At each projection step additional computational savings can be obtained by reducing the number of CC after each projection. This is possible since many of the open channels in the interaction region become highly closed as ρ increases and can be dropped from the CC equations at large ρ (see Appendix for more details). In summary, we perform two sets of scattering calculations: (1) a standard adiabatic calculation on the ground electronic state PES by solving eqn (2), and (2) a non-adiabatic calculation including two electronic states by solving eqn (1). The non-adiabatic 2 × 2 calculation is only required in the interaction region where the excited electronic state is important ρ ≤ ρ_DAT, for ρ > ρ_DAT both sets of scattering calculations are identical and are performed on the ground electronic state (except that the radial solutions are complex valued for the non-adiabatic case). The radial matching points (ρ_DAT, ρ_m and ρ_f), the number of CC, l_max, m_max, E_cut and other parameters are determined from extensive convergence studies and are discussed in detail in the Appendix.

3 Results

3.1 Conical intersection between ground and first excited states

The ab initio adiabatic Born–Oppenheimer electronic PES for the ground electronic state (V^a₁) is plotted in Fig. 2 for several values of ρ. At large ρ = 30a₀ (panel d), the three diatomic potentials (blue regions) for the different arrangement channels K + Rb₂ and the two identical Rb + KRb are clearly visible and well separated from each other by energy barriers (red regions). Diatomic vibrational motion corresponds to motion perpendicular to the blue contour lines while rotational corresponds to motion along the contour lines. As ρ decreases the Rb + KRb minimum energy collision pathway approaches collinear Rb–KRb geometries (the equator) indicated by the black arrows pointing to location A. At ρ = 20a₀ (panel c), the KRb arrangement channels merge at collinear geometries (dark blue regions and black arrows labeled A–B). We note that the conical intersection between the ground and first excited adiabatic electronic states lies in the highly repulsive (red) regions in panels (c) and (d) but is energetically accessible in the interaction region and its location is indicated by the red dot in panels (a) and (b). As ρ approaches ρ = 16a₀ (panel b), the minimum energy collision pathway moves away from collinear along the C_2v (T-shaped) configuration with the K nuclei located symmetrically between the two Rb (black arrow pointing to location B). Also the low energy green (0.0 eV) region that encircles the conical intersection (red dot) is energetically accessible (from Table 1 the energy of KRb (v = 0, j = 0) is 0.00459 eV (53.3 K)). At ρ = 12a₀ (panel a) the C_2v potential well (dark blue labeled by C) corresponds the global minimum of the 3D KRb₂ adiabatic ground state PES while two additional symmetric wells (dark blue labeled by E) occur on either side of the conical intersection. These two symmetric wells correspond to bent geometries and the saddle point (transition state) between them and the C_2v well is labeled by D. From Table 1 we see that the energy of Rb₂ (v = 0, j = 0) is 0.028084 eV (325.9 K) which lies 272.64 K above the Rb + KRb collision threshold and is closed for the ultracold collisions considered in this work. Thus the minimum energy pathway discussed above corresponds to purely elastic Rb + KRb collisions. Several of the important geometries discussed above and their corresponding energies are also listed in Table 1. Specifically, the global energy minimum on V^a₁ (−0.2358 eV) occurs for the C_2v geometry labeled by C in panel (a). The corresponding geometry is listed in Table 1 in three different sets of internal coordinate systems: hyperspherical, internuclear and Jacobi. The minimum energy (−0.2318 eV) for the bent geometries labeled by E in panel (a) are also tabulated as well as the minimum energy (−0.1682 eV) of the C_2v conical intersection (red dot).

Fig. 2 Contour plots of the stereographic projection of the ground adiabatic electronic state ab initio Born–Oppenheimer potential energy surface (PES) for KRb₂ are shown for several values of the hyperradius ρ = 12, 16, 20, and 30a₀ in panels (a)–(d), respectively. All energies are relative to the minimum energy of the asymptotic ground electronic state of Rb + KRb. The minimum energy pathway for collisions of Rb with KRb is indicated by the black arrows. Points of interest are indicated by the labels A through E. Dark blue corresponds to regions of low potential energy and red corresponds to high potential energy regions. The red dots in panels (a) and (b) denote the location of the conical intersection between the adiabatic ground and first excited electronic state (not plotted). In all panels, the contours start at 0.5 eV and decrease in steps of 0.1 eV down to 0 eV. Additional low-lying energy contours include −0.1, −0.2 and −0.22 eV in panel (a), −0.1 and −0.15 eV in panel (b), and −0.05 and −0.075 eV in panel (c). The radial blue dashed lines label the azimuthal hyperangle ϕ and the blue dashed outer and inner circles indicate the hyperangle values of θ = 53.1° and θ = 10°, respectively. The outer black circle denotes the equator of the hypersphere θ = 90° (collinear geometries).

Fig. 3 plots both the adiabatic (blue) and non-adiabatic (red) APH surface function energies as a function of ρ (i.e., the eigenvalues ε^Jpq_t(ρ_ξ) of eqn (7)). For better visualization, only every 100th adiabat is plotted. Asymptotically for large ρ, the surface function energies approach the diatomic KRb (v, j) and Rb₂ (v, j) rovibrational energies. The solid and dashed horizontal brown lines denote the energy of Rb + KRb (v = 0, j = 0) and K + Rb₂ (v = 0, j = 0), respectively. As expected the non-adiabatic (red) curves approach and eventually become identical with the adiabatic (blue) curves for large ρ > 20.0a₀ where the excited adiabatic electronic state becomes highly closed. For small ρ the excited adiabatic electronic state becomes energetically accessible (i.e., where the short-long dashed black curve lies below the horizontal brown curve for 10 < ρ < 13.5a₀). In addition, the conical intersection is readily encircled even when it lies above the Rb + KRb (v = 0, j = 0) threshold energy for 13.5 < ρ < 19a₀ (e.g., the low energy green region in Fig. 2 panel (b)). In this region we see significant differences between the adiabatic (blue) and non-adiabatic (red) curves in Fig. 3. These differences are due to the geometric phase (sign change) that occurs across the C_2v symmetry line (negative x-axis) between the conical intersection (red dot) and equator (outer black circle) in Fig. 2 panels (a) and (b). The points labeled A through E in Fig. 3 correspond to those in Fig. 2. Specifically, the label A denotes the large ρ region (ρ = 30a₀), B labels the intermediate ρ region where the conical intersection is readily encircled but is energetically inaccessible, C labels the region of the global C_2v minimum energy well, D labels the transition state between the global minimum and the two symmetric wells that occur for bent geometries denoted by E (see Fig. 2 panel (a)). The vertical series of short horizontal lines on the left edge of the plot denote the vibrational bound state energies of the cone states which are localized inside the cone (V^a₂) of the excited adiabatic electronic state. The E and O label the bound states of even and odd exchange symmetry, respectively. These vibrational bound states were computed numerically in full-dimensionality using the Numerov propagator on just the uncoupled adiabatic excited electronic state PES (V^a₂). The inset plot shows a zoomed-in region near the Rb + KRb threshold. The vertical energy scale on the inset plot has been converted to collision energy in units of Kelvin relative to the Rb + KRb (v = 0, j = 0) threshold energy. We see that the density of cone states near threshold (E_c < 1.16 K) is relatively low. A total of six states (four even plus two odd) lie in this energy range. A primary goal of the current study is to investigate the possibility of long-lived Feshbach⁵³ type scattering resonances associated with these cone states (see Section 3.3).

Fig. 3 The APH surface function energies are plotted as a function of the hyperradius ρ. Only every 100th adiabat is plotted. The blue and red adiabats correspond to the adiabatic and 2 × 2 diabatic calculations, respectively. The thick solid and dot-dashed black curves plot the minimum energy of the ground V₁ and excited V₂ adiabatic electronic states at each ρ, respectively. The horizontal solid brown line denotes the energy of the Rb + KRb (v = 0, j = 0) elastic collision channel (53.3 K). For reference, the dashed horizontal brown line denotes the energy of the closed K + Rb₂ (v = 0, j = 0) collision channel (325.9 K). The series of energy levels labeled by the E (even exchange symmetry) and O (odd exchange symmetry) are the three-dimensional vibrational bound state energies computed on the excited electronic state V₂ (i.e., cone states). All energies are relative to the minimum energy of the asymptotic adiabatic ground electronic state of Rb + KRb. The labels A through E denote points for discussion (see also Fig. 2). The inset shows the zoomed region near threshold and the vertical scale is plotted in collision energy (K) relative to the Rb + KRb (v = 0, j = 0) energy.

3.2 Scattering cross sections and rate coefficients

Fig. 4 and 5 plot the elastic cross sections for Rb + KRb (v = 0, j = 0) → Rb + KRb (v = 0, j = 0) as a function of collision energy in Kelvin. The red and blue colored curves correspond, respectively, to the even and odd exchange symmetry of the identical Rb nuclei. All of the cross sections plotted in Fig. 4 were computed on the single adiabatic ground electronic state PES only. The thin solid and dashed curves correspond to total angular momentum J = 0 (i.e., s-partial wave since J = j + l and j = 0) and J = 1⁻ (i.e., p-partial wave), respectively. Summing over the two J values give the total cross sections plotted as the thick solid curves. As expected, in the ultracold limit only the s-partial wave contributes and the thin and thick curves become identical. For higher collision energies, the p-partial wave begins to contribute: around 3 × 10⁻⁵ K for even symmetry and 2 × 10⁻⁴ K for odd symmetry. Small bumps occur in the total cross sections due to the p-wave shape resonances (which lie near 10⁻⁴ K, see also Section 3.3). At higher collision energies many more partial waves contribute to the total cross section but these were not computed in this work since our focus here is on the ultracold regime. The cross sections plotted in Fig. 5 are for J = 0 only and the dashed and solid curves correspond to the adiabatic and 2 × 2 diabatic calculations, respectively. We see that the non-adiabatic ultracold cross sections (solid curves) of both symmetries (red and blue) are suppressed relative to the adiabatic (dashed) ones. As we will show below, this suppression is due to destructive non-adiabatic quantum interference effects due to the geometric phase associated with the conical intersection. We note that in both Fig. 4 and 5 all of the cross sections are smooth in the ultracold regime. At higher energies, significant oscillatory structure is observed due to the appearance of scattering resonances (e.g., collision complexes). These resonances will be analyzed in more detail in Section 3.3.

Fig. 4 Elastic cross sections for Rb + KRb (v = 0, j = 0) → Rb + KRb (v = 0, j = 0) are plotted as a function of collision energy. These cross sections are computed on the single adiabatic ground electronic state PES of KRb₂. The red and blue curves correspond to the even and odd exchange symmetry of the identical Rb nuclei, respectively. The thin solid red and blue curves are computed with total angular momentum J = 0 whereas the dashed red and blue curves correspond to J = 1⁻. The thick solid red and blue curves are the total elastic cross sections summed over J = 0 and J = 1⁻.

Fig. 5 Same as in Fig. 4 for J = 0 except that here the dashed curves correspond to the elastic cross sections computed on the single adiabatic ground electronic state PES and the solid curves correspond to the 2 × 2 diabatic calculations (the red and blue colors correspond to even and odd exchange symmetry, respectively).

To understand the intriguing quantum interference effects observed in Fig. 5, we express the total scattering amplitude in terms of the direct (f̃_direct) and looping (f̃_loop) scattering amplitude contributions: .^54–56 The cross sections are proportional to the square modulus of the total scattering amplitude

where we have expressed the complex scattering amplitudes in polar form f̃_direct = f_direct exp(iδ_direct), f̃_loop = f_loop exp(iδ_loop) and Δ is the relative phase Δ = δ_direct − δ_loop. The direct and looping scattering amplitudes can be computed from the total adiabatic and non-adiabatic scattering amplitudes via and .⁵⁶ If the magnitudes of the direct and looping scattering amplitude are similar in magnitude, we can set them equal (i.e., f_direct ≈ f_loop = f) so that eqn (8) simplifies to

|f̃_total|² ≈ f ² (1 + cos Δ). (9)

We see from eqn (9) that if Δ approaches an even or odd multiple of π then |f̃_total|² ≈ 2f ² or 0, respectively. That is, we have maximum constructive or destructive quantum interference that effectively turns the collision on or off.^23,24 The integral π nature of Δ is due to the unique properties of ultracold collisions. Specifically, it follows from Levinson's theorem for which the scattering phase shifts are shown to preferentially approach integral values of π as the wavenumber k goes to zero.⁵⁷ In practice, the direct and looping magnitudes are not always identical and Δ is not an exact multiple of π. These conditions vary with collision energy and exchange symmetry as well as with the initial and final scattering states involved. They are also sensitive to the PES and its associated bound state structure. External fields can also be used to “tune” the degree of quantum interference and control the collision outcome.⁵⁸ We note that eqn (8) and (9) ignore the geometric phase, when the geometric phase associated with the conical intersection is taken into account there is an additional π phase shift in the relative phase. That is, the sign of the interference term containing cos Δ changes sign (i.e., cos Δ→ −cos Δ) and the nature of the quantum interference reverses (i.e., constructive becomes destructive and vice versa).

Fig. 6 plots the ratio of the square modulus of the direct and looping scattering amplitudes as a function of collision energy for both even (solid red) and odd (solid blue) exchange symmetries. The corresponding cos Δ are also plotted with dashed red and blue curves. These quantities are computed from the adiabatic and non-adiabatic scattering amplitudes for the cross sections plotted in Fig. 5. We see that in the ultracold limit, the cos Δ for both the even and odd exchange symmetries approach +1 (i.e., an even integral multiple of π). On the other hand, for this particular case the ratio of the square modulus of the ultracold scattering amplitudes are smaller than unity (they approach 0.054 and 0.035 for even and odd symmetry, respectively). Thus, the quantum interference effects are not maximum but range between approximately 1.52 (1.408) and 0.588 (0.661) for even (odd) symmetry using the general expression of eqn (8) with ±cos Δ. These ranges correspond to cross section ratios of 1.52/0.588 ≈ 2.58 (1.408/0.661 ≈ 2.13) for even (odd) symmetry. Indeed the ratios of the dashed and solid red (blue) curves in the ultracold limit in Fig. 5 are consistent with these ratios (as they should be): 7.322 × 10⁻¹²/2.835 × 10⁻¹² ≈ 2.58 (4.289 × 10⁻¹¹/2.014 × 10⁻¹¹ ≈ 2.13). For higher collision energies, we see that the cos Δ swing to −1 around 0.5 to 1 mK then rapidly back to +1. Thus, the quantum interference reverses and in the energy ranges where cos Δ < 0 we see that the adiabatic (dashed) cross sections lie below the 2 × 2 diabatic ones (solid) in Fig. 5. Fig. 6 shows that the direct scattering amplitude dominates at ultracold energies where the ratio is small. As the collision energy increases, the ratio for even (odd) symmetry rapidly approaches unity near 1 mK (0.5 mK), increases to approximately 10 (70) by 1.6 mK (0.65 mK) and then rapidly decreases again. In the energy ranges where the ratio is greater than unity, the looping scattering amplitude is dominate. The direct scattering amplitude corresponds to the minimum energy pathway denoted by the C–D–E arrows in Fig. 2 panel (a) that then turns around at E and reverses direction (or more generally encircles the conical intersection an even number of times). The looping scattering amplitude corresponds to the same C–D–E pathway but instead of turning around at E it continues to encircle the conical intersection via the transition state across the x-axis between the two symmetric wells labeled by E (or more generally encircles the conical intersection an odd number of times).

Fig. 6 The solid red and blue curves plot the ratios of the modulus of the looping (f̃_loop) and direct (f̃_direct) scattering amplitudes averaged over the scattering angle as a function of collision energy. The dashed red and blue curves plot the averaged cos Δ where Δ is the phase difference between the two scattering amplitudes. The modulus ratios are specified on the log scale on the left vertical axis and the cos Δ values are specified on the linear scale on the right vertical axis. These results are for zero total angular momentum (J = 0) and the red and blue colors correspond to even and odd exchange symmetry, respectively.

The total rate coefficients summed over the appropriately weighted even and odd cross sections are plotted in Fig. 7 (the black solid and dashed curves). The weights are determined by the number of nuclear spin states of the two identical Rb nuclei for each exchange symmetry. The nuclear spin of ⁸⁷Rb is I = 3/2 so that the total number of nuclear spin states for the two Rb is (2I + 1)² = 16. Out of this total, there are (2S + 1)(S + 1) = 10 symmetric and (2S + 1)S = 6 antisymmetric nuclear spin states under exchange of the two identical Rb nuclei. The electronic wave function for Rb₂ is a ¹Σ⁺_g state which is symmetric under exchange of the Rb. Since the total molecular wavefunction is a direct product of the electronic, nuclear motion and nuclear spin wave functions, we know that for a symmetric electronic wave function the product of the nuclear motion and nuclear spin wave function must be antisymmetric since ⁸⁷Rb is a fermion. Thus, the even (odd) nuclear motion wave functions must be multiplied by the antisymmetric (symmetric) nuclear spin wave functions. Upon averaging over the initial and summing over the final symmetrized nuclear spin states, we find that the weights of the even and odd cross sections are 3/8 and 5/8, respectively. The experimentally¹⁴ measured rate coefficient (6.8 × 10⁻¹¹cm³ s⁻¹) at 1 µK is also plotted (black dot with error bars) in Fig. 7. The adiabatic rate coefficient (dashed black curve) at 1 µK is 5.12 × 10⁻¹¹cm³ s⁻¹ and the corresponding 2 × 2 diabatic rate coefficient (solid black curve) is 2.42 × 10⁻¹¹ cm³ s⁻¹. The ratio of the two theoretical rate coefficients is 2.1 which is consistent with the ratios discussed above due to constructive and destructive quantum interference.

Fig. 7 The total elastic rate coefficients for Rb + KRb (v = 0, j = 0) → Rb + KRb (v = 0, j = 0) are plotted as a function of collision energy for zero total angular momentum (J = 0). The dashed curves correspond to the calculations on the single adiabatic ground electronic state PES and the solid curves correspond to the 2 × 2 diabatic calculations. The experimentally measured rate coefficient is the solid black dot with error bars. The black curves include weighted contributions from both the even and odd cross sections that treats all of the nuclear spin states of the two identical Rb nuclei. The blue curves include only the odd cross section and correspond to a single symmetric nuclear spin state (no weights).

Summing over both the even and odd cross sections with the appropriate statistical nuclear spin weights is the standard approach when all of the nuclear spin states are present and no distinction or measurement of the nuclear spins is made in the initial and final states of the collision. However, in the Rb + KRb experiments the initial states of the Rb and KRb are prepared in specific hyperfine states (and therefore specific nuclear spin states): Rb (F = 1, m_F = 1, m^Rb_I = 3/2) and KRb (m^K_I = −4, m^Rb_I = 1/2).¹⁴ Due to energy constraints and angular momentum conservation, the above mixture is expected to be stable against inelastic, spin-changing collisions.¹⁴ Thus, experimentally the nuclear spin states are the same before and after the collision. However, during the collision, the deep attractive potential of the three-body Rb–K–Rb complex allows for the two identical Rb nuclei to freely move and exchange places with each other (see Fig. 2). Therefore, within the interaction region we also have Rb(m^Rb_I = 1/2) and KRb(m^K_I = −4, m^Rb_I = 3/2) and we can construct symmetric and antisymmetric combinations: ψ_± = |Rb(m^Rb_I = 3/2)〉 |KRb(m^Rb_I = 1/2)〉 ± |Rb(m^Rb_I = 1/2)〉 |KRb(m^Rb_I = 3/2)〉. From the above symmetry discussions, we know that the even and odd cross sections correspond to the antisymmetric (ψ₋) and symmetric (ψ₊) nuclear spin states, respectively. The blue solid (2 × 2) and blue dashed (adiabatic) curves in Fig. 7 are the rate coefficients computed from the corresponding odd cross sections in Fig. 5 for a single symmetric (ψ₊) nuclear spin state (i.e., no averaging over initial and summing over final nuclear spin states). We see that this combination gives the best agreement with the experimental rate coefficient at 1 µK: 7.41 × 10⁻¹¹ cm³ s⁻¹ for the adiabatic and 3.56 × 10⁻¹¹ cm³ s⁻¹ for the 2 × 2 diabatic. The universal value is 7.0 × 10⁻¹¹ cm³ s⁻¹ which is in excellent agreement with the experimental one but it ignores the 3-body interaction region entirely and the associated quantum interference effects that are present in both of the adiabatic and 2 × 2 diabatic calculations. These quantum interference effects are sensitive to the details of the PES and as shown recently for the Li + LiNa reaction, scaling the 3-body ab initio contributions to the PES by a few percent can significantly alter the nature and magnitude of the interference.²⁵ While such a study is beyond the scope of the present work, we anticipate that tuning the ab initio 3-body contribution within its estimated numerical uncertainty will likely reverse the nature of the quantum interference (i.e., swap the dashed and solid curves in Fig. 7) bringing the 2 × 2 diabatic rate coefficient into better agreement with the experimental value.

3.3 Collision lifetimes

The collision lifetimes are computed from the scattering matrix and its energy derivative using Smith's Q matrix defined asThe lifetimes are given by the trace of Q: τ = Tr(Q) and physically represent the difference between the time for the collision to occur with the interaction potential and the time for the Rb and KRb to pass each other with no interaction. The matrix Q is Hermitian and in practice we diagonalize Q and sum over the diagonal eigenvalues to compute the trace.⁵⁹ In the diagonal (eigenphase) representation, the S matrix is given by

S_ij = δ_ij exp(i2κ_i(E)). (11)

Substituting the above expression for S into eqn (10) we obtain the diagonal representation for QNear a resonance the most dominant phase shift κ_i takes the Breit–Wigner formwhere κ_i(0) is the local background phase shift, Γ_r the resonance width and E_r the resonance energy. Assuming the background phase shift changes slowly with energy, we use eqn (13) in eqn (12) to obtainThus Q_i has the well-known Breit–Wigner or Lorentzian form and for most resonances the largest (dominant) eigenlifetime dominates so that Tr Q ≈ Q_i. Computing Tr Q on a sufficiently dense energy grid allows for the fitting of eqn (14) to extract the resonance energy and width parameters. In practice, we apply a polynomial interpolation to the computed set of Tr Q which is very good at identifying the locations of narrow resonance peaks.^59,60 We then add additional energy points near the predicted peaks until the interpolated resonance profile is unchanged and fully characterized. In this study, we are interested in resolving potentially very narrow (long-lived) 3-body complexes including those that might be localized in the excited adiabatic electronic state PES (i.e., cone states). Thus, a dense energy grid was used consisting of 4352 collision energies for each exchange symmetry and both the adiabatic and 2 × 2 diabatic calculations for a total of 4352 × 4 = 17 408. A detailed discussion of the energy grid and computational methods for the log-derivative propagation are given in the Appendix.

The collision lifetimes for Rb + KRb (v = 0, j = 0) computed on the single adiabatic ground electronic state are plotted in Fig. 8 as a function of collision energy. The inset in Fig. 8 shows the low energy range below 0.1 mK. Both the even (red) and odd (blue) lifetime curves are smooth functions of energy with no sign of complex formation (resonances). The negative lifetimes are due to the attractive potential well in the interaction region which accelerates the collision leading to a smaller collision time relative the same process but with no interaction potential. At around 6 mK the first 3-body collision complex (resonance) occurs for even symmetry (red) with a peak lifetime of 200 ns. The first complex of odd symmetry (blue) occurs near 8 mK with a peak lifetime of 1.2 µs. Two more notable complexes of odd symmetry occur near 65 and 100 mK with peak lifetimes of 2.9 µs and 650 ns, respectively. The lifetimes for almost all of the other resonances lie below 100 ns. The negative dip in the lifetime at 0.106 K is due to the opening of the KRb (v = 0, j = 1) channel. Similar dips are also seen at 0.319 K and 0.639 K which are the j = 2 and j = 3 threshold energies. The j = 4 threshold lies at 1.065 K which is above the 1 K maximum collision energy considered in this work. The lifetimes for J = 1⁻ are plotted in Fig. 9. In the ultracold regime they are smooth functions of collision energy except that now we see the appearance of shape resonances near 50 µK due to p-wave (l = 1) scattering. As for J = 0 the 3-body resonances start to appear near 6 mK. The lowest energy and most prominent resonance is the 4.3 µs one of even symmetry (red). Two notable resonances of odd symmetry (blue) occur near 65 and 80 mK with lifetimes of 2.5 and 3.9 µs, respectively. The lifetimes for the majority of the other resonances are smaller than 200 ns. The three negative dips in the lifetimes at the j = 1, 2 and 3 threshold energies are clearly visible. The lifetimes for J = 0 computed using the 2 × 2 diabatic approach are plotted in Fig. 10 as a function of collision energy. As in Fig. 8 for the adiabatic calculation, the inset shows smooth lifetime behavior in the ultracold regime. No 3-body resonances appear until about 6 mK where we see the appearance of an odd symmetry (blue) complex with a lifetime of 120 ns. The most prominent complex is of odd symmetry (blue) and occurs near 30 mK with a lifetime of 1.7 µs. As in Fig. 8 and 9, the three negative dips at the three KRb rotational thresholds are clearly visible. While there are qualitative similarities, the 2 × 2 diabatic resonance spectrum in Fig. 10 is very different than the adiabatic resonance spectrum in Fig. 8. Interestingly, no long-lived complexes are seen that might correspond to the bound vibrational states localized in the excited adiabatic electronic state (i.e., cone states). This is most likely due to a lack of good coupling or overlap between the scattering wave function and the localized cone state wave functions.

Fig. 8 Lifetimes for elastic Rb + KRb (v = 0, j = 0) collisions are plotted as a function of collision energy for zero total angular momentum J = 0. The lifetimes are computed on the adiabatic ground electronic state from Tr(Q) where Q is Smith's collision lifetime matrix. The red and blue curves correspond to even and odd exchange symmetry, respectively. The inset shows the ultracold energy range below 0.1 mK. The lifetimes are smooth functions of energy from 1 nK up to 5 mK. Starting at around 6 mK 3-body collision complexes (resonances) begin to appear and increase in density at higher energies. The maximum lifetimes of the three longer lived collisions complexes of odd symmetry (blue) are denoted near the top of the plot.

Fig. 9 Same as in Fig. 8 except for total angular momentum J = 1⁻. Shape resonances occur near 50 µK and the one for even exchange symmetry (red) is most prominent. Starting at around 6 mK 3-body collision complexes (resonances) begin to appear and increase in density at higher energies. The maximum lifetimes of the three longer lived collisions complexes of odd (blue) and two even (red) symmetry are denoted near the top of the plot.

Fig. 10 Same as in Fig. 8 except that the calculations are based on the 2 × 2 diabatic approach. The inset shows the ultracold energy range below 0.1 mK. The lifetimes are smooth functions of energy from 1 nK up to 5 mK. Starting at around 6 mK 3-body collision complexes (resonances) begin to appear and increase in density at higher energies. The maximum lifetimes of the two longer lived collisions complexes of even (red) and one of odd (blue) symmetry are denoted near the top of the plot.

The number of resonances counted within several energy intervals for each of the three resonance spectrums shown in Fig. 8–10 are tabulated in Table 2. From this table we see that the number of resonances for each exchange symmetry is about the same for all of the calculations and energy intervals. Also the number of resonances within the upper two energy intervals (0.1–0.5 K and 0.5–1.0 K) are about the same for all of the calculations and about two to three times larger than the number of resonances in the lowest energy interval (0.005–0.1 K). The total number of resonances for each calculation summed over both exchange symmetries over the full collision energy range 0.005 to 1.0 K are listed in bold in the bottom row of the table. We note that the total number of resonances for the 2 × 2 diabatic calculation is about 20% larger than the adiabatic one. Also, the total number of resonances for the adiabatic J = 1⁻ calculation is about twice that for the J = 0 adiabatic one. This is due to the doubling of states for J = 1⁻ (due to the two Ω = ±1 contributions) versus J = 0 (for which Ω = 0).³¹ The average density can be computed by dividing the total resonance counts by the full energy interval 0.995 K to obtain ρ = 117.6, 141.9 and 230.1 K⁻¹ for the adiabatic, 2 × 2 diabatic and adiabatic J = 1⁻ calculations, respectively. From Table 2 we see that the adiabatic J = 0 even exchange density is ρ = 58/0.995 = 58.3 K⁻¹ which is about twice the one reported for Rb + K₂ elastic collisions ρ = 29 K⁻¹.⁶⁰

Table 2 The number of resonances counted within several energy intervals are tabulated for the adiabatic, 2 × 2 diabatic and adiabatic J = 1⁻ calculations for each exchange symmetry. The density can be calculated from the total number of resonances reported for each calculation listed in bold in the bottom row of the table. The densities are ρ = 117.6, 141.9 and 230.1 K⁻¹, respectively

Energy range (K)	Adiabatic J = 0			2 × 2 diabatic J = 0			Adiabatic J = 1^−
	Even	Odd	Total	Even	Odd	Total	Even	Odd	Total
0.005–0.1	10	10	20	10	11	21	22	19	41
0.1–0.5	26	26	52	28	31	59	45	48	93
0.5–1.0	22	23	45	32	29	61	48	47	95
0.005–1.0	58	59	117	70	71	141	115	114	229

---

From the density of states, we can estimate an overall average lifetime of the collision complexes via the well-known statistical RRKM lifetime expression^16,17,19

where ρ is the density of states and N₀ is the number of open exit channels (for the elastic collisions considered here N₀ = 1). From the average resonance densities quoted above, the corresponding statistical lifetime estimates are τ = 5.64, 6.81 and 11.0 ns for the adiabatic, 2 × 2 diabatic and adiabatic J = 1⁻ calculations, respectively. The above densities and lifetimes do not include the degenerate magnetic quantum numbers (the current calculations are field free so the magnetic quantum number energies are degenerate). Including the magnetic quantum number degeneracies increases the density of states and lifetime estimates by approximately two orders of magnitude.⁶⁰ Thus, our RRKM τ based estimates including the magnetic quantum number degeneracies are approximately τ = 564, 681 and 1100 ns for the adiabatic, 2 × 2 diabatic and adiabatic J = 1⁻ calculations, respectively. These estimates are consistent with previous statistical theory estimates (237 ns) for Rb + KRb²⁰ and are approximately 10³ smaller than the experimentally reported one 390 µs.¹⁴

4 Conclusions

A first principles based treatment of ultracold elastic collisions of Rb + KRb (v = 0, j = 0) in full dimensionality was used to investigate non-adiabatic effects originating from the first excited adiabatic electronic state which exhibits a conical intersection with the ground electronic state. This excited electronic state is energetically accessible within the interaction region even at ultracold collision energies approaching absolute zero and must be included in the theoretical calculations. The unique properties of ultracold collisions often leads to the realization of maximum possible quantum interference between the pathways that encircle the conical intersection an even or odd number of times. This maximal quantum interference can effectively turn the collision outcome on or off. Furthermore, the geometric phase associated with the conical intersection changes the nature of the quantum interference at ultracold energies from constructive to destructive and vice versa.

The non-adiabatic quantum interference effects were quantified by comparing the non-adiabatic scattering results to the adiabatic ones computed on the single adiabatic ground electronic state. It was found that the non-adiabatic ultracold elastic cross sections and rate coefficients were suppressed by a factor of approximately two relative to the adiabatic ones. This suppression is due to destructive quantum interference. The theoretical ultracold rate coefficients are in good agreement with the experimental value. In particular the adiabatic rate coefficient for a single symmetric nuclear spin state gave excellent agreement whereas the corresponding 2 × 2 diabatic one was approximately two times smaller due to the destructive interference. We expect that improvements in the ab initio 3-body interactions will likely reverse the nature of the quantum interference. A high resolution collision energy grid was used to explore the possible formation of 3-body complexes (resonances) within the interaction region. The non-adiabatic calculations include the excited electronic state where localized complexes might form within the inverted cone shaped potential energy surface (i.e., Feshbach type resonances). Such localized resonant states could lead to very long lifetimes. Interestingly, no evidence for long-lived cone state complexes were found in the present study. The resonance spectrum for both the adiabatic and non-adiabatic scattering calculations were smooth functions of collision energy in the ultracold regime. In the cold regime (T > 1 mK) a few longer lived complexes were observed with lifetimes in the 2 to 4 µs range. The lifetimes for the majority of resonances were below 200 ns. The density of resonances was analyzed and a statistical analysis gave complex lifetime estimates that were consistent with previous theoretical estimates.²⁰ These statistical estimates are approximately three orders of magnitude smaller than the experimentally reported lifetime.¹⁴

A full dimensional first-principles based quantum mechanical treatment including non-adiabatic effects from multiple electronic states is required in order to accurately describe ultracold collisions. However, questions remain as to the origin and mechanism behind the mysteriously long-lived collision complexes reported experimentally for the elastic Rb + KRb collisions.¹⁴ As suggested recently, perhaps spin–spin and/or spin-rotation interactions might be responsible.¹⁵ In future work we plan to include these interactions within the present theoretical framework and investigate this possibility. Advances in both theory and experiment are continuing to provide novel insight into the mechanisms of ultracold collisions. We hope that this work will provide a theoretical benchmark to help unravel the intriguing mysteries of the ultracold regime.

Conflicts of interest

There are no conflicts to declare.

Data availability

Data for this article are available within the body of the article and in the form of Fig. 1–10 and Tables 1, 2.

Supplementary information (SI) is available in the form of a subroutine which evaluates the two adiabatic and diabatic electronic potential energy surfaces as a function of the three internuclear distances. See DOI: https://doi.org/10.1039/d6cp00282j.

Appendix

Computational details

In this appendix we discuss the computational details and converged parameters used in the quantum scattering calculations. A series of full scattering calculations are performed for different values of each convergence parameter and the optimal values are determined. The optimal values typically involve a trade-off between the desired accuracy and computational cost. The optimal values for the matching points were determined to be: ρ_DAT = 27.05a₀, ρ_m = 60.07a₀, ρ_f = 325.1a₀. The corresponding radial grid in ρ spanned 432 logarithmetically spaced sectors starting at ρ₀ = 7.0a₀ and extending to ρ_m. In the Delves region (ρ_m < ρ ≤ ρ_f), a logarithmetically spaced radial grid was used with 3379 sectors. The optimal number of CC used in the adiabatic single-surface calculations were 3600 for ρ ≤ ρ_DAT, 800 for ρ_DAT < ρ ≤ ρ_m, and 300 for ρ_m < ρ ≤ ρ_f. For the J = 1⁻ propagation 7125 channels were used for ρ ≤ ρ_DAT and twice the number of the J = 0 CC were used for ρ > ρ_DAT. Within each sector the computational grid in ρ for the CC propagation is subdivided into many smaller sub-steps. The number of sub-steps within each sector varies with ρ depending upon the local de-Broglie wavelength and sector width.³² The resulting propagation grid across each sector is uniform with a grid spacing of Δρ = 2.2 × 10⁻⁴a₀ in the APH region and Δρ = 3.5 × 10⁻⁴a₀ in the Delves region. The optimal number of CC used in the 2 × 2 diabatic calculations were 4600 for ρ ≤ 24.97a₀. The number of CC were then truncated to 2000 and the propagation continued to ρ = ρ_DAT where the DAT is applied. The non-adiabatic propagation beyond ρ = ρ_DAT is identical to the adiabatic one except that a complex valued version of the log-derivative propagator must be used. The DAT transformation at ρ_DAT is performed by applying a sequence of n_DAT = 160 unitary transformations involving overlap projections of the 2 × 2 diabatic solutions with the last projection being onto the adiabatic ground state solutions. To preserve unitarity, the mixing angle β is scaled by a parameter λ which various from 0 to 1 and is uniformly discretized into n_DAT values.³⁴ Thus, an additional set of n_DAT 2 × 2 diabatic surface function calculations are required at ρ_DAT. Each of these calculations uses a slightly different PEM computed using an appropriately scaled β → λ_iβ where i = 1 to n_DAT. We note that λ = 1 corresponds to the diabatic basis and λ = 0 corresponds to the adiabatic basis with a diagonal PEM.

The optimal size of the primitive basis sets for the 5D angular solutions of eqn (7) increases with ρ. Eight regions in ρ were used: 7.0 ≤ ρ ≤ 12.06a₀, 12.06 < ρ ≤ 20.05a₀, 20.05 < ρ ≤ 23.99a₀, 23.99 < ρ ≤ 28.15a₀, 28.15 < ρ ≤ 35.41a₀. 35.41 < ρ ≤ 41.12a₀, 41.12 < ρ ≤ 50.2a₀. and 50.2 < ρ ≤ 60.07a₀. The optimal l_max and m_max in each region are: l_max = 159, 191, 223, 255, 287, 287, 319, and 351; and m_max = 380, 410, 440, 480, 520, 560, 600, and 640, respectively. For the J = 0 adiabatic calculations, the dimension of the angular 5D Hamiltonian matrix in each region is given by (l_max + 1) × (2m_max + 1) = 121 760, 157 632, 197 344, 246 016, 299 808, 322 848, 384 320 and 450 912, respectively. The size of these matrices is reduced dramatically using an energy cut-off and applying SDT at each value of ρ. The corresponding E_cut values at the center of each region in ρ are respectively: 3.92, 0.73, 0.91, and 1.07 eV for the first four regions and 1.19 eV for the last four regions. After applying SDT the dimension of the largest angular 5D Hamiltonian matrix in each of the eight regions is 48 242, 47 620, 67 344, 93 225, 121 393, 133 552, 164 333, and 197 054, respectively. For the J = 1⁻ adiabatic calculations, the dimension of the primitive basis is twice that for J = 0 due to contributions from the two Ω = ±1 blocks.³¹ Thus, after applying SDT the dimension of the largest angular 5D Hamiltonian matrix for J = 1⁻ in each of the eight regions is twice the values quoted above for J = 0. For the J = 0 2 × 2 diabatic calculations, the size of the primitive basis is also twice that of the adiabatic ones due to the second electronic state. After applying SDT the dimension of the largest angular 5D Hamiltonian matrix in each of the four regions within the interaction region (ρ < 27.05a₀) is 94 077, 79 173, 124 909, and 174 953, respectively. These matrices are smaller than twice the size of the J = 0 adiabatic ones due to the localized (conical) nature of the adiabatic excited electronic state PES.

As mentioned above, the angular 5D Hamiltonian matrices at each sector (ρ_ξ) are diagonalized using the parallel sparse matrix eigensolver library PARPACK.⁵² The J = 0 and J = 1⁻ adiabatic matrices are real symmetric but the 2 × 2 diabatic matrices are complex Hermitian. Fortunately, the explicit construction of these large matrices is not required, only the action of the matrix on a vector is needed. The sparse matrix-vector multiply has been efficiently parallelized and discussed in detailed in previous work.^31,34 The eigensolutions were computed using a 128 node CPU machine where each Intel Sapphire Rapids node contains two 56-core 1.9 Ghz CPUs with 128 Gb of High Bandwidth Memory (HBM). The HBM was found to significantly accelerate the calculations giving a 4× speed-up relative to a conventional non-HBM CPU node. The total number of eigensolutions were computed in blocks using the projection technique.³² For example, the 3600, 4200 and 7125 eigensolutions for the J = 0 adiabatic, J = 0 2 × 2 diabatic and J = 1⁻ adiabatic calculations were computed in blocks of 240, 105, and 125 at a time using 15, 40, and 57 projection steps, respectively. The CPU time for the sparse eigensolver scales as t_cpu = N^d × n_sol where n_sol denotes the number of desired solutions, N is the dimension of the matrix and d is the scaling law³³ which is nominally d = 2. However, due to the banded structure of the sparse matrix and efficient parallelization of the matrix-vector multiply, the actual scaling exponent on the HBM nodes is d = 1.25. For a given exchange symmetry, the computation time for the J = 0 adiabatic surface function calculations for the interaction region (ρ < 27.05a₀) was 1.2 days. For the intermediate range (27.05 ≤ ρ ≤ ρ_m = 60.07) the calculations required an additional 2.8 days for a grand total of 4 days. The J = 1⁻ adiabatic surface function calculations required approximately twice the CPU time as the J = 0 calculations for a total of 8 days. The J = 0 2 × 2 diabatic surface function calculations in the interaction region required approximately five times the CPU time as the J = 0 adiabatic calculations for a total of 6 days. The DAT at ρ = ρ_DAT = 27.05a₀ required an additional 160 surface function calculations totaling 7 days. Thus, a grand total of 13 days was required for the 2 × 2 diabatic surface function calculations. In summary, for a given exchange symmetry the entire set of surface function calculations required 4, 8, and 13 days for the J = 0 adiabatic, J = 1⁻ adiabatic and J = 0 2 × 2 diabatic calculations, respectively. Thus, the total computation time including both even and odd exchange symmetries was approximately 2 × (4 + 8 + 13) = 50 days. The computational cost for the CC radial propagation is discussed below.

The long-range (60.07 < ρ ≤ 325a₀) Delves basis calculations are quite fast since the rotational wave functions are the well known spherical harmonics. Thus, only a one-dimensional Schrödinger equation in the internuclear distance must be solved numerically. A Numerov propagator⁶¹ is used to compute the vibrational wave functions for each diatomic arrangement channel and rotational quantum number. At ρ = ρ_m = 60.07a₀ the 300 Delves ro-vibrational basis functions include the first four vibrational levels of KRb: v = 0, 1, 2, 3 and for each vibrational state the corresponding maximum rotational quantum numbers for KRb are 84, 72, 57 and 36, respectively. For the Rb₂ arrangement channel, the v = 0 and 1 vibrational states are included with corresponding maximum rotational quantum numbers of 58 and 32 for even exchange symmetry, respectively. Asymptotically at ρ = ρ_f the open ro-vibrational Jacobi basis functions for a 1 K collision energy include the single KRb(v = 0, j = 0) elastic collision channel and the higher-lying KRb(v = 0, j = 1, 2, 3) rotational channels with threshold collision energies of 0.106, 0.319, and 0.639 K, respectively. The Delves basis, overlap and potential coupling matrix calculations are parallelized and can be completed within one day using a single CPU node with 128 cores.

The CC radial propagation for J = 0 and a given exchange symmetry was performed using a high-resolution energy grid consisting of 4352 energies partitioned in the following manner: 10, 20, 1440, 1440, and 1442 points spaced logarithmically between 1 to 10 nK, 10 to 100 nK, and linearly between 100 nK to 1 mK, 1 to 10 mK, and 0.1 to 1 K, respectively. The J = 1⁻ calculations were restricted to just the adiabatic case and were done at lower resolution using half the number of collision energies: 2176 for each exchange symmetry for a total of 4352. A GPU version of the log-derivative propagator was implemented using Nvidia's CUDA BLAS library. The GPU acceleration was crucial for enabling the treatment of a large number of collision energies and the relatively large number of coupled channels in the interaction region (3600, 4200 and 7125 for the adiabatic, 2 × 2 diabatic and J = 1⁻ adiabatic, respectively). A large GPU machine (Chicoma) was used which has 118 nodes with 4 Nvidia A100 GPUs each for a total of 472. The A100 tensor cores are ideally suited for the dense BLAS operations (e.g., matrix–matrix multiplication and matrix inversion) required by the log-derivative propagator. A single A100 GPU gives a factor of 55 speed-up relative to a conventional CPU node running a parallel ScaLAPACK version of the log-derivative propagator on 16 cores. The double-precision floating point performance of the matrix-matrix multiply routine (dgemm) on a single A100 was measured to be 17.4TFLOPS. Comparing this value to a single CPU performance of 57GFLOPS or 16 CPUs running ScaLAPACK 450GFLOPS, the GPU advantage is significant. Since each compute node on Chicoma contains four A100 GPUs, we divided the 4352 energies into eight groups of 544 and each set of 544 energies were distributed across 34 compute nodes. Thus, each of the four GPUs per node was assigned four energies (i.e., 16 energies total per node). Since each energy calculation is independent, the energies can be distributed across all available GPU resources. In practice, due to resource availability, we typically submitted the 8 batches of 544 energies sequentially two batches at a time using 68 nodes (4 × 68 = 272 A100 GPUs). The two batches of 544 energies for the adiabatic 3600 CC propagations across the interaction region took approximately 24 hour. Since the number of CCs were truncated at large ρ (see Section 2.2), the long-range propagation was significantly accelerated and took approximately one day as well. Thus, two days were required for two batches of 544 energies which resulted in eight days total for the full set of 8 × 544 = 4352 energies. The 2 × 2 diabatic calculations require complex matrix operations and that together with the slightly larger number of CCs (4200) resulted in about a factor of two longer propagation times. The J = 1⁻ propagations required nearly twice as many CCs (7125) but were done at half as many energies so that they required about three times longer propagation times. In summary, a single exchange symmetry required about 8, 16, and 24 days for the adiabatic, 2 × 2 diabatic and J = 1⁻ adiabatic propagations, respectively. Thus, about three months were required to complete the entire set of propagations for the 4352 × 5 = 21 760 collision energies.

Acknowledgements

B. K. K. acknowledges that part of this work was done under the auspices of the US Department of Energy under Project no. 20240256ER of the Laboratory Directed Research and Development Program at Los Alamos National Laboratory. This work 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). S. K., H. L., and J. K. acknowledge support from the US Air Force Office of Scientific Research Grant No. FA9550-21-1-0153, US National Science Foundation Grant No. PHY-2409425, and the Gordon and Betty Moore Foundation.

Notes and references

1. N. Y. Yao, A. V. Gorshkov, C. R. Laumann, A. M. Läuchli, J. Ye and M. D. Lukin, Phys. Rev. Lett., 2013, 110, 185302.
2. K.-K. Ni, T. Rosenband and D. D. Grimes, Chem. Sci., 2018, 9, 6830–6838.
3. S. F. Yelin, K. Kirby and R. Côté, Phys. Rev. A: At., Mol., Opt. Phys., 2006, 74, 050301.
4. W. B. Cairncross, D. N. Gresh, M. Grau, K. C. Cossel, T. S. Roussy, Y. Ni, Y. Zhou, J. Ye and E. A. Cornell, Phys. Rev. Lett., 2017, 119, 153001.
5. N. Balakrishnan, J. Chem. Phys., 2016, 145, 150901.
6. M. S. Safronova, D. Budker, D. DeMille, D. F. J. Kimball, A. Derevianko and C. W. Clark, Rev. Mod. Phys., 2018, 90, 025008.
7. W. B. Cairncross and J. Ye, Nat. Rev. Phys., 2019, 1, 510–521.
8. S. L. Cornish, M. R. Tarbutt and K. R. A. Hazzard, Nat. Phys., 2024, 20, 730–740.
9. Y. Liu and K.-K. Ni, Annu. Rev. Phys. Chem., 2022, 73, 73–96.
10. T. Karman, M. Tomza and J. Pérez-Ríos, Nat. Phys., 2024, 20, 722–729.
11. P. Gersema, K. K. Voges, M. Meyer zum Alten Borgloh, L. Koch, T. Hartmann, A. Zenesini, S. Ospelkaus, J. Lin, J. He and D. Wang, Phys. Rev. Lett., 2021, 127, 163401.
12. R. Bause, A. Schindewolf, R. Tao, M. Duda, X.-Y. Chen, G. Quéméner, T. Karman, A. Christianen, I. Bloch and X.-Y. Luo, Phys. Rev. Res., 2021, 3, 033013.
13. J. N. Byrd, J. A. Montgomery and R. Côté, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 82, 010502.
14. M. A. Nichols, Y.-X. Liu, L. Zhu, M.-G. Hu, Y. Liu and K.-K. Ni, Phys. Rev. X, 2022, 12, 011049.
15. Y.-X. Liu, L. Zhu, J. Luke, M. C. Babin, M. Gronowski, H. Ladjimi, M. Tomza, J. L. Bohn, T. V. Tscherbul and K.-K. Ni, Nat. Chem., 2025, 17, 688–694.
16. R. A. Marcus, J. Chem. Phys., 1952, 20, 352–354.
17. R. A. Marcus, J. Chem. Phys., 1952, 20, 355–359.
18. J. C. Light, Discuss. Faraday Soc., 1967, 44, 14–29.
19. R. D. Levine, Molecular reaction dynamics, Cambridge University Press, 2009.
20. M. Mayle, B. P. Ruzic and J. L. Bohn, Phys. Rev. A: At., Mol., Opt. Phys., 2012, 85, 062712.
21. M. Mayle, G. Quéméner, B. P. Ruzic and J. L. Bohn, Phys. Rev. A: At., Mol., Opt. Phys., 2013, 87, 012709.
22. R. Bause, A. Christianen, A. Schindewolf, I. Bloch and X.-Y. Luo, J. Phys. Chem. A, 2023, 127, 729–741.
23. B. K. Kendrick, J. Hazra and N. Balakrishnan, Nat. Commun., 2015, 6, 7918.
24. B. K. Kendrick, J. Hazra and N. Balakrishnan, Phys. Rev. Lett., 2015, 115, 153201.
25. B. K. Kendrick, H. Li, M. Li, S. Kotochigova, J. F. E. Croft and N. Balakrishnan, Phys. Chem. Chem. Phys., 2021, 23, 5096–5112.
26. B. K. Kendrick, J. Chem. Phys., 2021, 154, 124303.
27. H. J. da Silva, B. K. Kendrick, H. Li, S. Kotochigova and N. Balakrishnan, J. Phys. Chem. Lett., 2025, 16, 6171–6177.
28. M. Morita, M. B. Kosicki, P. S. Zuchowski, P. Brumer and T. V. Tscherbul, Phys. Rev. A, 2024, 110, L021301.
29. R. T. Pack and G. A. Parker, J. Chem. Phys., 1987, 87, 3888–3921.
30. B. Kendrick and R. T. Pack, J. Chem. Phys., 1996, 104, 7475–7501.
31. B. K. Kendrick, R. T. Pack, R. B. Walker and E. F. Hayes, J. Chem. Phys., 1999, 110, 6673–6693.
32. B. K. Kendrick, J. Chem. Phys., 2000, 112, 5679–5704.
33. B. K. Kendrick, J. Chem. Phys., 2001, 114, 8796–8819.
34. B. K. Kendrick, J. Chem. Phys., 2018, 148, 044116.
35. H.-J. Werner, P. J. Knowles, G. Knizia, F. R. Manby and M. Schütz, et al., MOLPRO, version 2015, a package of ab initio programs, 2015, see https://www.molpro.net.
36. W. Müller, J. Flesch and W. Meyer, J. Chem. Phys., 1984, 80, 3297–3310.
37. P. Fuentealba, L. Szentpály, H. Stoll, F. Fraschio and H. Preuss, J. Mol. Struct.: THEOCHEM, 1983, 93, 213–219.
38. P. S. Zuchowski and J. M. Hutson, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 81, 060703.
39. H. Werner and P. J. Knowles, J. Chem. Phys., 1985, 82, 5053–5063.
40. P. J. Knowles and H.-J. Werner, Chem. Phys. Lett., 1985, 115, 259–267.
41. H. Werner and P. J. Knowles, J. Chem. Phys., 1988, 89, 5803–5814.
42. P. J. Knowles and H.-J. Werner, Chem. Phys. Lett., 1988, 145, 514–522.
43. D. Simah, B. Hartke and H.-J. Werner, J. Chem. Phys., 1999, 111, 4523–4534.
44. T. Ho and H. Rabitz, J. Chem. Phys., 1996, 104, 2584–2597.
45. O. T. Unke and M. Meuwly, J. Chem. Inf. Model., 2017, 57, 1923–1931.
46. A. Pashov, O. Docenko, M. Tamanis, R. Ferber, H. Knöckel and E. Tiemann, Phys. Rev. A: At., Mol., Opt. Phys., 2007, 76, 022511.
47. C. Strauss, T. Takekoshi, F. Lang, K. Winkler, R. Grimm, J. Hecker Denschlag and E. Tiemann, Phys. Rev. A: At., Mol., Opt. Phys., 2010, 82, 052514.
48. B. Johnson, J. Comput. Phys., 1973, 13, 445–449.
49. Z. Darakjian and E. F. Hayes, J. Chem. Phys., 1990, 93, 8793–8797.
50. G. A. Parker, R. B. Walker, B. K. Kendrick and R. T. Pack, J. Chem. Phys., 2002, 117, 6083–6102.
51. J. C. Light, I. P. Hamilton and J. V. Lill, J. Chem. Phys., 1985, 82, 1400–1409.
52. R. Lehoucq, D. Sorensen and C. Yang, ARPACK Users' Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods, Society for Industrial and Applied Mathematics, 1998.
53. H. Feshbach, Ann. Phys., 1958, 5, 357–390.
54. J. C. Juanes-Marcos, S. C. Althorpe and E. Wrede, Science, 2005, 309, 1227–1230.
55. S. C. Althorpe, J. Chem. Phys., 2006, 124, 084105.
56. S. C. Althorpe, T. Stecher and F. Bouakline, J. Chem. Phys., 2008, 129, 214117.
57. N. Levinson, Phys. Rev., 1953, 89, 755–757.
58. J. Hazra, B. K. Kendrick and N. Balakrishnan, J. Phys. Chem. A, 2015, 119, 12291–12303.
59. B. Kendrick and R. T. Pack, J. Chem. Phys., 1996, 104, 7502–7514.
60. J. F. E. Croft, N. Balakrishnan and B. K. Kendrick, Phys. Rev. A, 2017, 96, 062707.
61. B. R. Johnson, J. Chem. Phys., 1977, 67, 4086–4093.

---