Magnetic tuning of ultracold barrierless chemical reactions

While attaining external field control of bimolecular chemical reactions has long been a coveted goal of physics and chemistry, the role of hyperfine interactions and dc magnetic fields in achieving such control has remained elusive. We develop an extended coupled-channel statistical theory of barrierless atom-diatom chemical reactions, and apply it to elucidate the effects of magnetic fields and hyperfine interactions on the ultracold chemical reaction Li($^2\text{S}_{1/2}$) + CaH($^2\Sigma^+$) $\to$ LiH($^1\Sigma^+$) + Ca($^1\text{S}_{0}$) on a newly developed set of ab initio potential energy surfaces. We observe large field effects on the reaction cross sections, opening up the possibility of controlling ultracold barrierless chemical reactions by tuning selected hyperfine states of the reactants with an external magnetic field.

While attaining external field control of bimolecular chemical reactions has long been a coveted goal of physics and chemistry, the role of hyperfine interactions and dc magnetic fields in achieving such control has remained elusive.
We develop an extended coupledchannel statistical theory of barrierless atom-diatom chemical reactions, and apply it to elucidate the effects of magnetic fields and hyperfine interactions on the ultracold chemical reaction Li( 2 S 1/2 ) + CaH( 2 Σ + ) → LiH( 1 Σ + ) + Ca( 1 S0) on a newly developed set of ab initio potential energy surfaces. We observe large field effects on the reaction cross sections, opening up the possibility of controlling ultracold barrierless chemical reactions by tuning selected hyperfine states of the reactants with an external magnetic field.
Using external electromagnetic fields to control chemical reactivity is a central goal of chemical physics [1,2], which stimulated the development of new research avenues ranging from mode-selective chemistry [1] and coherent control [2] to the study of stereodynamics and vector correlations in molecular collisions [3][4][5] and ultracold controlled chemistry [6,7]. Molecular chemical reactions are most readily controlled at ultralow temperatures, where the reactants can be prepared in single internal and motional quantum states [8], which maximizes the effects of external electromagnetic fields [9] and allows for the manifestation of quantum phenomena, which would otherwise be obscured by thermal averaging, such as threshold and resonance scattering [7,8,10], tunnelling [7,11], and interference [12,13].
Recent examples include the observation of resonance scattering in low-temperature He * + H 2 [10], He + NO [14], and NO + H 2 [15] collisions, stereodynamical control of low-temperature H 2 + HD collisions in merged molecular beams [4,5], and chemical reactions in trapped ensembles of alkali-metal dimers [16,17], and atom-dimer mixtures [18]. The vast majority of the previous control studies have focused on the rovibrational and nuclear spin degrees of freedom of the reactants. In particular, the chemical reaction KRb + KRb → K 2 + Rb 2 can be efficiently suppressed by preparing the reactants in the same rotational and nuclear spin states [19] and stimulated by applying an external electric field, which modifies the pwave centrifugal barrier preventing the reaction of two identical fermionic molecules [16,20].
Recent experimental advances in laser cooling and trapping [21,22] have led to the production of dense, trapped ensembles of molecular radicals (i.e. molecules with nonzero electron spins) such as CaF( 2 Σ + ) [23][24][25][26], SrF( 2 Σ + ) [27], YbF( 2 Σ + ), [28] and SrOH( 2 Σ + ) [29]. Cotrapping of these molecules with ultracold alkali-metal atoms [27,28] would open up the fascinating prospect of studying spin-selective ultracold controlled chemistry [6,30,31]. Specifically, the electron spins of the reactants can be polarized in an external magnetic field to form a nearly spin-pure state in the entrance reaction channel corresponding to, e.g., the maximum possible total spin S of the reaction complex [6,30]. Because such high-S states are typically non-reactive, the chemical reaction of spin-aligned reactants [(A(↑) + B(↑)] will be suppressed compared to that of spin-antialigned However, theoretical studies of the effects of spin polarization, hyperfine interactions, and external magnetic fields on atom-molecule chemical reactions have been limited to reactions of weakly bound Feshbach molecules [32,33], save for a recent model study of ultracold NH-NH reactive scattering in a magnetic field [34], which did not include the hyperfine structure of NH and focused on collisions of fully spin-polarized molecules. As a result, the effects of hyperfine interactions and magnetic fields on ultracold reaction dynamics remain unexplored, limiting our ability to use the fields as a tool to control chemical reactivity at ultralow temperatures.
Here, we develop a theoretical approach to ultracold reaction dynamics in a magnetic field based on a rigorous coupled-channel statistical (CCS) model [1][2][3]38]. The CCS model postulates the existence of a long-lived reaction complex formed temporarily when the reactants get trapped in a long-lived resonance state [1-3, 39, 40], a powerful idea that forms the basis for quantum threshold models [41,42] and quantum defect theories [43][44][45]. The CCS approach rigorously accounts for the multichannel nature of molecular wavefunction in the entrance and exit reaction channels [1][2][3] and it has been successfully applied to calculate low-temperature inelastic [46] and reactive [15,41,45,47,48] collision rates. Building on the previous work, we extend the CCS approach to explicitly include the effects of hyperfine interactions and external magnetic fields in the entrance reaction channel, which allows us to explore the magnetic field dependence of the reaction cross sections. We exemplify the extended CCS approach by applying it to the chemical reaction Li + CaH → LiH + Ca on a newly developed set of ab initio potential energy surfaces (PESs). Our fieldfree results are in good agreement with experiment at T = 1 K [50]. We find that the reaction can be efficiently suppressed by tuning the hyperfine states of the reactants with an external magnetic field, opening up the possibility for controlling ultracold spin-dependent chemical reactions. Minimizing atom-molecule reaction rates is essential for efficient sympathetic cooling, in which molecules are immersed in a gas of ultracold atoms and refrigerated by elastic collisions [9,51,53,54]. Our results thus show that sympathetic cooling of chemically reactive 2 Σ radicals could be facilitated by applying external magnetic fields. Theory.
The original CCS theory [1,2] relates the state-to-state reaction probability P γAγB →γ ′ A γ ′ B to the capture probabilities in the entrance and exit reaction channels p γAγB and p γ ′ subject to capture boundary conditions [1][2][3] as described in the Supplemental Material [56]. Here, R is the atommolecule separation vector, r joins the nuclei in the diatomic molecule, and µ andL are the reduced mass and orbital angular momentum of the collision complex. The asymptotic HamiltoniansĤ A andĤ B account for the rotational, fine, and hyperfine structure of the reactants in the presence of an external magnetic field [56], which is crucial for controlling ultracold reaction dynamics, as shown below.
The atom-molecule interaction operator in Eq. (1) is given byV (R, r) = S,MS V S (R, r, θ)|SM S SM S |, where V S (R, r, θ) are the adiabatic atom-molecule PESs in the entrance reaction channel calculated ab initio as described in the Supplemental Material [56], S is the total spin of the reaction complex and M S is the projection of S on the magnetic field axis. Figure 1 Results. We now apply the extended CCS methodology to explore the effect of tuning the Zeeman states of the reactants on the ultracold chemical reaction Li + CaH → LiH + Ca. Figure 2 shows the collision energy dependence of the reaction cross section calculated for the different spin states of the reactants |Sm S with the hyperfine structure omitted for the moment. The reaction cross section for spin-antialigned reactants CaH decreases with the collision energy E C as expected for the Langevin cross section (σ R ≃ E −1/3 [57]). By averaging the dependence σ R (E) over a Maxwell-Boltzmann distributions of collision energies, we obtain the reaction rate in quantitative agreement with the measured value of 3.6 × 10 −10 cm 3 /s [50].
As shown in Fig. 2, the reaction of spin-aligned reactants | 1 2 , 1 2 Li + | 1 2 , 1 2 CaH is suppressed by four orders of magnitude compared to that of spin-antinaligned initial states. The spin-aligned reaction nevertheless occurs through the intramolecular spin-rotation and intermolecular magnetic dipole-dipole interactions, which flip the total spin of the Li-CaH complex in the entrance reaction channel [30,31]. Because these interactions are weak, the spin-aligned reaction rate is small, and is comparable to that of non-reactive spin relaxation in . We next explore the effects of external magnetic fields and hyperfine interactions on chemical reactivity. Figure 3 shows the magnetic field dependence of reaction cross sections for the different initial hyperfine states of Li and CaH [see Figs. 1(b)-(c)]. We observe that certain combinations of initial hyperfine states are far more reactive than others: In particular, changing the initial state from |11 CaH + |2, 1 Li to |11 CaH + |2, −2 Li enhances the reaction by a factor of 5. This suggests the possibility of controlling ultracold reaction rates by tuning the hyperfine states of the reactants, which could be realized experimentally via radiofrequency and/or optical pumping.
Remarkably, as shown in Fig. 3, the reactivities of selected initial hyperfine states are extremely sensitive to the magnetic field strength, which opens up the prospect of controlling ultracold barrierless chemical reactions with external magnetic fields. To gain insight into the extreme field dependence of the reaction cross sections, we observe that the nuclear spin degrees of freedom (DOF) do not directly participate in the reaction dynamics, which is governed instead by the spin DOF. This implies, in the spirit of the degenerate internal states approximation [58,59], that the matrix elements of the atom-molecule PES are diagonal in the nuclear spin projections m IA and m IB . The reaction cross sections are given in terms of the exact S-matrix elements The initial Zeeman states |γ i m i (i = A, B) are linear combinations of the electron and nuclear spin states We assume that the S-matrix elements on the right are independent of B (which is approximately true as shown in Fig. 3) and they are different from zero only if m SA and m SB correspond to the reactive singlet state (S = 0) as discussed above. The magnetic field dependence of the reaction cross section is thus encapsulated in the hyperfine mixing coefficients C mS i mI i ,γimi (B).
We now illustrate the hyperfine model (3) where σ R 1 2 ,− 1 2 →f is the reaction cross section in the absence of the hyperfine structure (upper trace in Fig. 2).
As shown in Fig. 3, the reaction cross sections predicted by the hyperfine model decrease as a function of the applied magnetic field, in qualitative agreement with the CCS results. The suppression is due to the decoupling of the electron and nuclear spins: As shown in Fig. 1(c), the initial state |2, 1 Li correlates with the nonreactive Zeeman state |m S = 1/2 Li in the high-field limit, leading to a decrease of the contribution of the reactive state |m S = − 1 2 Li . The reaction cross section scales as B −2 due to the mixing coefficient C − 1 2 3 2 ,21 (B) ≃ B −1 . We note that the hyperfine model predicts a less steep decline of the reaction cross sections with the field. This is expected because the hyperfine model neglects the effects of the magnetic field on the S-matrix elements, which are likely to become more pronounced at higher fields [58,59].
In conclusion, we have extended the rigorous CCS model of barrierless chemical reactions [1][2][3] to include the hyperfine structure of open-shell reactants and their interactions with external magnetic fields. We have applied the model to explore the effects of hyperfine interactions and magnetic fields on the dynamics of the prototypical barrierless chemical reaction CaH + Li → LiH + Ca. Our calculated reaction rates agree with experiment [50] and display a dramatic dependence on the external magnetic field, which could be used to facilitate sympathetic cooling of chemically reactive 2 Σ molecules with alkali-metal atoms [9,60].
We expect our approach to be readily applicable to a wide range of ultracold barrierless chemical reactions of current experimental interest, including those involving molecular ions [61][62][63][64] and alkaline-earth halides SrF and CaF [27,54,65].
We are grateful to Gerrit Groenenboom, Roman Krems, and Masato Morita for encouraging discussions. This section provides an overview of the extended CCS approach. We begin by introducing the CC equations in Sec. IA and describing the procedure of applying the boundary conditions in Sec. IB. Sec. IC describes further technical details pertaining to the evaluation of the matrix elements of the atom-molecule interaction and of the orbital angular momentum of the collision complex.

A. Numerical solution of CC equations: Reaction cross sections and capture probabilities
The CCS capture probability in the entrance reaction channel is given by [1][2][3] where γ A , γ B and γ ′ A , γ ′ B stand for the initial and final Zeeman states of the reactants, l and l ′ are the corresponding orbital angular momenta, and M is the space-fixed (SF) projection of the total angular momentum J of the collision complex on the magnetic field axis, which is conserved for reactions in magnetic fields. The total reaction cross section is obtained by summing the entrance channel capture probabilities (5) over a range of orbital angular momenta l and total angular momentum projections M where k γAγB = 2µE C is the wavevector in the incident collision channel and E C is the collision energy. We note that the reaction cross section can be obtained from the fully state-to-state cross section by summing over the final LiH + Ca product states γ ′ A and γ ′ B . The S-matrix elements in Eq. (5) are obtained from the radial solutions F M αAαB JΩ (R) of the coupled-channel (CC) equations at total energy E [1-3] subject to the capture boundary conditions as described below. The CC equations (7) describe atom-molecule scattering in the entrance reaction channel in the presence of an external magnetic field. In Eq. (7) are body-fixed (BF) basis functions for the overall rotational motion (|JM Ω ) and the internal degrees of freedom of molecule A (α A ) and atom B (α B ), including the rotational angular momentum N A , the electron spinŝ S A andŜ B and the nuclear spinsÎ A andÎ B , with K NA , Σ A , Σ B , Σ IA , and Σ IB being the projections of N A , S A , S B , I A , and I B on the atom-diatom separation vector R chosen as the z-axis of the BF coordinate frame [4]. The matrix elements are evaluated as described in Sec. IC below.

B. Boundary conditions
We solve the CC equations numerically by constructing the log-derivative matrix Y = Ψ −1 Ψ, where Ψ is the wavefunction matrix, and propagating it from a small value of R = R c out to the asymptotic region [1][2][3]. We choose an initial value of the capture radius R = R c inside the reaction complex region and initialize the complex symmetric log-derivative matrix as [1][2][3] 10) using the multichannel Wentzel-Kramers-Brillouin (WKB) boundary conditions [3]. The entrance channels of chemical reactions that occur on multiple PESs (such as the Li-CaH reaction considered here) typically include the highly attractive as well as strongly repulsive PESs, leading to two qualitatively different types of adiabatic channels illustrated in Fig 1. The reactive channels decrease in energy with decreasing R < R c , whereas the nonreactive channels show the opposite trend. Both types of channels can be treated on an equal footing using the Airy boundary conditions [3,5]. Following Ref. [5], we initialize the elements of the diagonal matrix Y d as where W ′ (R c ) 1/3 is the real root of x 3 = W ′ (R c ), and W ′ (R c ) is the derivative of the adiabatic eigenvalue W (R c ) (an eigenvalue of matrix W), which is positive (negative) for reactive (inelastic) channels [see Fig. 1]. The scaled wavefunction φ(x) is a solution of the Airy equation for the adiabatic channels linearly extrapolated into the reaction complex region R < R c [5] For W ′ (R c ) > 0 (reactive channels) and R ≪ R c , we have x → −∞ and the wavefunction ratio in Eq. (11) takes the form [5] where Bi(x) and Ai(x) are the Airy functions, which oscillate in the limit of large negative x.
For W ′ (R c ) < 0 (inelastic channels) and R ≪ R c , we have x → +∞. Retaining only the asymptotically decaying Airy function Ai(x), we obtain the wavefunction ratio in Eq. (11) as [5] In practice, the asymptotic expressions (13)- (14) give sufficiently accurate results for |x| ≥ 5. However, in our numerical calculations a small fraction of adiabatic channels has |x| < 5, making it necessary to apply numerically exact expressions for the Airy functions.
At a large atom-molecule distance R = R a , we match the log-derivative matrix Y(R a ) to the standard incoming and outgoing wave boundary conditions to obtain the scattering S-matrix [1,2] where γ is a compound index for γ A , γ B , k γ = [2µE C ] 1/2 is the incident wavevector, E C is the collision energy, and h (±) l (x) are the spherical Hankel functions. The asymptotic solutions for closed channels are given by [2] [ where i l (x) and k l (x) are the modified spherical Bessel functions.

C. Matrix elements
We now turn to the technical details of the evaluation of the matrix elements in the CC equations (7). The asymptotic Hamiltonian may be written as [4] H as =Ĥ A +Ĥ B , whereĤ A andĤ B are the asymptotic Hamiltonians of the reactants [molecule A( 2 Σ) plus atom B( 2 S)]. The molecular Hamiltonian is given by [6,7] whereN A is the rotational angular momentum,Ŝ A and I A are the electron and nuclear spins with space-fixed (SF) projectionsŜ ZA andÎ ZA [I A = S A = 1/2 for CaH(X 2 Σ)],Î A ⊗Ŝ A is a tensor product ofÎ A andŜ A , Y 2−q (r) is a spherical harmonic describing the orientation of the molecular axisr A in the SF frame, and B e , γ, b, and c are the rotational, spin-rotation, and hyperfine constants. We neglect the weak nuclear spin-rotation interaction [6]. The atomic Hamiltonian includes the hyperfine coupling of the electron and nuclear spins parametrized by the atomic hyperfine constant A B , and the interaction of the atomic spin with an external magnetic field B. Writing the asymptotic Hamiltonian (18) as a sum of field-free and Zeeman termŝ and taking advantage of the direct-product structure of the BF basis set (8), we obtain The matrix elements on the right can be evaluated as described in our previous work [4,6]. Diagonalization of H as at B > 0 produces unphysical Zeeman eigenstates, which do not affect low-temperature collision dynamics provided a sufficient number of total angular momentum eigenstates J max is included in the CC basis [4].
To calculate the matrix elements of the orbital angular momentum operator in Eq. (7) in the body-fixed angular momentum basis, we express the latter in the form WhileL 2 can be expressed in terms of the raising and lowering operators for all angular momenta involved as done, e.g., in Ref. [4], the resulting expressions are rather cumbersome due to the presence of two additional nuclear spin operatorsÎ A andÎ B . To simplify the evaluation of the angular momentum matrix elements, it is convenient to define the orbital angular momentum of the atommolecule system in the absence of the nuclear spin The matrix elements of this operator can be evaluated as described in our previous work [4]. To incorporate the nuclear spins, we combine Eqs. (23) and (24) and use the fact that the nuclear spin operators commute withL 4 to obtain Expressing the scalar products of angular momentum operators via the raising and lowering operators, e.g., L 4 ·Î A =L 4z ·Î Az + 1 2 (L 4+ ·Î A− +L 4− ·Î A+ ), the evaluation of the matrix elements of the operatorL 2 (23) in the basis (8) reduces to straightforward angular momentum algebra described, e.g., in Ref. [4].

II. AB INITIO CALCULATIONS AND PES FITTING
To compute the singlet (S = 0) and triplet (S = 1) PES of the Li-CaH reaction complex, we used highlevel multi-reference configuration interaction (MRCI) and coupled-cluster methods with single, double, and noniterative triple excitations as implemented in the MOLPRO code [8].
The triplet PES is calculated as described in our previous work [9] at the CCSD(T) level of theory [10]. To compute the singlet (S = 0) PES, we took into account the multi-reference character of the electronic wavefunction by using the multi-reference configuration interaction (MRCI) method [11] with single and double excitations (MRCISD) and Davidson corrections (+Q) to approximately account for contributions of higher excitations. The MRCISD+Q calculations were started from the reference orbitals obtained at the state-averaged complete active space self-consistent-field (sa-CASSCF) level treating all S states on the same footing. The active space contained 4 A ′ and 1 A ′′ orbitals and 6 orbitals were correlated but kept doubly occupied (5 A ′ and 1 A ′′ ). The frozen core of the Ca atom was composed of 4 A ′ and 1 A ′′ orbitals. We used the augmented, correlation consistent triple-zeta (aug-cc-pvtz) basis for H [12], a quadruple-zeta basis (aug-cc-pvqz) for Li , and a valence quadruple-zeta (cc-pvqz) basis for Ca [13].
The ab initio calculations were performed on a twodimensional grid of R and θ, with θ ∈ 0 • −180 • in steps of 5 • and R ∈ 3.5−30 a 0 [8]. To facilitate the calculation of the matrix elements of the atom-molecule interaction operatorV (R, r, θ) = SMS V S (R, θ, r)|SM S SM S |, where S is the total spin of the reaction complex, we expand the adiabatic potential energy surfaces V S (R, θ, r) in Legendre polynomials P λ (cos θ) Because we neglect the weak nuclear spin-dependent interactions that depend on R (such as the interaction of the nucelar spins with the overall rotation of the reaction complex), the matrix elements of the interaction potential are diagonal in the nuclear spin quantum numbers Σ IA and Σ IB [see Eq. (8)]. As a result, the matrix elements of Eq. (26) are given by the expressions similar to Eqs. (30) and (31) of Ref. [4]. For use in quantum scattering calculations, the ab initio data points are expanded in Legendre polynomials (26) with λ max = 18 (for S = 0) and λ max = 14 (for S = 1). The resulting radial expansion coefficients V λ (R) are fit using the Reproducing Kernel Hilbert Space method of Rabitz and coworkers [14]. To avoid unphysical distortion of the fit, we damped the very high repulsive energies at small R.
Following our previous work [15], we invoke the rigidrotor approximation by freezing the internuclear distance of CaH at its equilibrium value r e = 2.0025 Å. This approximation provides quantitatively accurate capture probabilities for the Li-CaH chemical reaction on a single adiabatic PES [15] at a much reduced computational cost. However, in the context of the CCS model, the rigid-rotor approximation leads to all adiabatic potentials becoming repulsive (i.e. non-reactive) at sufficiently short R. To address this, we introduce the following modification of the isotropic part of the singlet (reactive) PES where R m is a matching point to the right of the potential minimum, where the first and second derivatives of the potential have opposite signs with the second derivative being negative, so as to ensure the decreasing behavior of Eq. (27) with decreasing R < R m [see Fig. 2(a)]. The modification replaces the short-range repulsive wall of the rigid-rotor potential with a function that decreases with R. This results in a one-parameter family of modified potentials parameterized by the values of R m . We have verified (see Sec. II below) that the calculated capture probabilities are insensitive to the choice of R m , thereby validating the procedure.

III. CONVERGENCE TESTS
We carried out a series of convergence tests to determine the optimal values of the asymptotic matching distance R a and the cutoff parameters N max and J max that determine the sizes of the rotational and total angular momentum basis sets. We use the following values of R a to obtain the capture probabilities and reaction cross sections converged to within 10-20%: E C : R a = 360 a 0 (E C = 10 −4 − 10 −3 cm −1 ), R a = 270 a 0 (E C = 10 −3 − 10 −2 cm −1 ), and R a = 200 a 0 (E C > 10 −2 cm −1 ), with a uniform grid step of 0.02 a 0 . The value of the capture radius R c was set to 3.84a 0 in all of the calculations.
At the lowest collision energies studied in this work (10 −4 < E C < 10 −3 cm −1 ) it is sufficient to truncate the total angular momentum basis at J max = 2 to produce results converged to < 5%. At higher collision energies, progressively higher values of J max were used, up to J max = 10 at E C = 7 cm −1 .
We also carried out convergence tests with respect to the maximum number of rotational states N max included in the basis set. The capture cross sections for the spin-antialigned initial states | 1 2 , 1 2 CaH + | 1 2 , − 1 2 Li are large and remarkably insensitive to N max as shown in Fig.  3. This suggests that anisotropic effects in the entrance channel of the Li+CaH reaction play a minor role. For these initial states, a minimal basis set with N max = 2 was used. In contrast, the small capture probabilities of spin-aligned reactants tend to be highly sensitive to the value of N max , making it necessary to employ much larger rotational basis sets with N max = 55. The large rotational basis sets are required to account for the large anisotropy of the Li-CaH interaction, as shown in our previous work on nonreactive spin relaxation in ultracold Li-CaH collisions [9].