Quantum dynamics of tunneling dominated reactions at low temperatures

Quantum state-resolved chemistry with cold and ultracold molecules continues to be a topic of active experimental and theoretical investigations [1–3]. This allows complete control on the initial state preparation, including, electronic, vibrational, rotational, spin and hyperfine levels as well as initial collision energy and angular momentum partial wave. Such controlled studies of chemical reactions between two KRb molecules or a KRb and a K atom have recently been reported by the JILA group [4]. Several other research groups are pursuing similar experiments on a variety of alkali metal, alkaline earth metal and other non-alkali metal atom systems [5–8]. Also, an experimental study of the benchmark F + H₂ reaction has recently been reported at temperatures as low as 12 K [9]. However, widespread studies of chemical reactions in the sub-Kelvin regime continue to be a challenge, in particular, for systems composed of non-alkali metal atoms.

Over the last 10–15 years numerous theoretical calculations of inelastic and reactive collisions in atom–molecule and molecule–molecule systems have been reported in the cold and ultracold regimes [1–3, 10, 11]. The vast majority of these studies on chemical reactions have resorted to a single or a few partial waves in the reactant channel. While only s-waves contribute in the limit of zero collision energy (p-wave for fermionic systems), for many systems, a large number of partial waves contribute to the cross sections in the mK and K regime. For systems with deep potential wells or moderately strong van der Waals interaction potential, the low temperature rate coefficients may be influenced by resonance structures supported by the interaction potential [10]. Whether these resonances survive the partial wave summation is not known for many systems. Such calculations are extremely demanding at low temperatures as large basis sets and propagation to large radial separations are required to yield converged results. Though the time-dependent wave packet methods [12–15] offer computational advantages for reactive scattering calculations they are problematic at very low energies and have not been proven to be accurate for cold and ultracold collisions.

The goal of this paper is two-fold. First, we present state-to-state reaction dynamics of the Li + HF $({{v}_{i}},\;{{j}_{i}})\;\to \;$ LiF(v_f , j_f ) + H reaction using a quantum close-coupling (CC) method focusing on moderate to low collision energies in light of recent experimental data [16] that showed qualitatively different behavior at low collision energies compared to previous calculations based on a wave packet method [17]. Second, we demonstrate usefulness of a recently developed approach [18] for ultracold reactions that invoke the quantum CC method at short-range and the multichannel quantum defect theory (MQDT) at long-range. The method yields full state-resolved cross sections at comparable level of accuracy as full CC calculations but at much reduced computational expense. Results are presented for HD formation in both D + H₂ and H + D₂ reactions.

The Li + HF reaction has been the topic of numerous experimental and theoretical studies over the last four decades. Experimental studies [16, 19–28] focused on measurements of its reaction cross section as well as determination of product angular distributions and alignment effects. Theoretical work included extensive investigation of its potential energy surfaces (PESs) [29–39] and quantum dynamics calculations based on both time-dependent and time-independent methods [17, 40–50]. It has long served as a prototype for reactions involving an alkali metal atom and a hydrogen halide. Due to the presence of three light atoms, high-level electronic structure calculations of its PESs have become possible in recent years [35–39]. The reaction is characterized by deep van der Waals potentials both in the incident and outgoing channels. We have previously investigated [48, 49] the dynamics of this reaction in the cold and ultracold regimes by including only the s-wave contribution in the incident Li + HF channel. These studies employed the PES of Aguado et al [37]. A more detailed dynamics study, including contributions from higher partial waves, is reported here for the collision energy range of 10⁻⁷ − 0.4 eV. Further, we investigate the sensitivity of the reaction dynamics to the three-body forces in the interaction potential by comparing results from two different ab initio PESs, one employed in our previous work and a newer potential reported in [38]. Results show very strong sensitivity to the choice of the potential, especially for low energies where quantum tunneling may play an important role in the reaction dynamics.

The choice of Li + HF system is motivated by a recent experimental study of the reaction by Bobbenkamp et al [16] who reported cross sections in the energy range 25 meV ≤ E_trans ≤ 131 meV (290−1520 K). Though the collision energy is much higher than in cold collisions, the experimental results indicated a rapid rise in cross section at low energies suggesting a vanishing threshold or an energy threshold less than 22 meV. Bobbenkamp et al [16] indicated the possibility of carrying out lower energy measurements by employing a magneto optical trap for preparing ultracold lithium atoms combined with a source of cold molecules based on rotating nozzle and electrostatic guiding.

Despite the large number of electronic structure calculations and theoretical investigations of the reaction dynamics, there seems to be considerable discrepancy between the predicted reaction cross sections at energies close to the threshold. This is, due in part, to the extreme sensitivity of the dynamics to fine details of the interaction potential, especially in the threshold regime where tunneling may play an important role. The computations of Laganá and coworkers [47] seem to predict a vanishing threshold for the reaction consistent with the experimental results of Bobbenkamp et al [16]. However, their results differed from the experimental data for collision energies above 0.2 eV. A more recent work from their group [51] that included limited CC calculations at two energies close to the threshold and the J-shifting approximation yielded results that are consistent with the experimental findings. In contrast, the time-dependent quantum calculations of Zanchet et al [17] predict a sharp threshold at energies below 10 meV though their cross sections agree with experiments for energies above 0.2 eV. These two calculations employed two different PESs for the reaction which differed in the height of the energy barrier and depth of the van der Waals potential well. The computations of Zanchet et al adopted a PES computed by Aguado et al [38] (referred to as the APW PES) which is also adopted in other recent theoretical studies [50] of this reaction that explored the effects of reagent rotational excitation on the reaction dynamics. The above results suggest that the dynamics of the reaction is very sensitive to the energy barrier. As the reaction occurs by tunneling at low energies, very small changes in the height and width of the barrier can strongly influence the dynamics in the threshold regime.

In an attempt to accurately characterize the reaction profile of the Li + HF reaction, recently, Fan et al [39] carried out CCSD(T) and CCSDT calculations of key reaction intermediates using basis sets as large as core correlated quintuple zeta. Unfortunately, their calculations were restricted to the reaction energetics and not the full PES surface suitable for dynamics studies. It was found that stationary points along the reaction path are very sensitive to the level of theory. They reported a value of −6.1 kcal/mol (−0.265 eV) for an entrance channel van der Waals complex (relative to Li + HF reactants) and a bent transition state at 4.5 kcal/mol (0.195 eV) above the reactants with an Li⋯H⋯F angle of ${{72.2}^{\circ }}$. The barrier height for the APW PES of Aguado et al [38] is 0.221 eV and the van der Waals minimum is at −0.243 eV. All energies are relative to energy zero for the bottom of the HF potential and excluding zero-point energy correction. The reaction is endoergic without the inclusion of the zero-point energy but becomes exoergic when it is included. The exorgicity for the APW PES is −0.08 eV which is the same as that obtained by Fan et al. These values are somewhat different from a previous version of the PES reported by Aguado et al [37] for the LiHF system. The barrier height of this PES is 0.251 eV and it has a smaller exoergicity of −0.01122 eV. This PES was used in our previous studies of Li + HF and its reverse reaction [48, 49] which indicated very small reactivity for Li + HF(v_i = 0, j_i = 0) initial state at energies below 1 K. It was also employed by Tscherbul and Krems [52] in their study of the reaction in an external electric field but the basis set adopted was not adequate to yield converged results. On the otherhand, our calculations present fully converged scattering results on APW PES of Aguado et al [38]. Table 1 compares the key features of the three PESs. The barrier height and depth of the van der Waals complex are given without the zero-point energy correction while the exoergicity includes zero-point energy correction.

Here we provide a more detailed investigation of the Li + HF system on the PES of Aguado et al [37] and the presumably more accurate APW PES to explore sensitivity of the reaction dynamics to details of the interaction potential. Further, we map out the energy dependence of the cross section in the threshold regime to compare against the experimental results of Bobbenkamp et al [16].

The paper is organized as follows. In section 2 we give a brief overview of the quantum scattering approach adopted for the reactive scattering calculations of the Li + HF reaction. Convergence studies of key reaction attributes as well as comparisons with previous calculations are discussed in section 3. Section 4 presents a detailed discussion of the results, including comparisons with experimental data. A brief discussion of the CC–MQDT method and its application to D + H₂ and H + D₂ reactions are given in section 5. A summary and conclusions are presented in section 6.

The quantum scattering calculations are carried out using the ABC code of Skouteris et al [53]. The ABC code solves the time-independent Schrödinger equation for an atom–diatom system in Delves-hyperspherical coordinates [54] which include three internal coordinates (hyperradius ρ and two hyperangles) and three external coordinates (three Euler angles). An adiabatic approach is adopted for the solution of the Schrödinger equation by discretizing the hyperradius into a large number of sectors. The solution of the reactive scattering problem involves two key steps: (i) eigenvalue problem involving the surface function hamiltonian (overall hamiltonian for the atom–diatom system excluding the radial kinetic energy operator in ρ) at each sector to evaluate the hyperspherical surface functions and the corresponding adiabatic energies, and (ii) propagation of the radial Schrödinger equations from small ρ within the classically forbidden region to a large asymptotic value in ρ. The reactive scattering boundary conditions are then applied at the last sector in ρ to evaluate the reactance matrix ${{{\bf K}}^{JP}}$ and the parity adapted S-matrix, ${\bf S}_{{{v}_{f}}{{j}_{f}}{{k}_{f}}\leftarrow {{v}_{i}}{{j}_{i}}{{k}_{i}}}^{JP}$. For an incident collision energy E_c , separate calculations are needed for each {J, P, q} triple, where J is the quantum number denoting the total angular momentum (vector sum of the orbital angular momentum and rotational angular momentum of the diatomic molecule in a given arrangement channel), P = (−1)^J is the triatomic inversion symmetry (sign change under co-ordinate inversion) and q = (−1)^j is the diatomic parity (identical particle permutation symmetry in the reactant and product diatomic molecule). The quantum numbers v and j represent the vibrational and rotational quantum numbers of diatomic species and the helicity quantum number k denotes the projection of J on the body fixed z-axis which is directed along the atom–molecule center-of-mass vector. Once the parity-adapted S-matrix is converted into standard helicity representation [53], ${\bf S}_{{{v}_{f}}{{j}_{f}}{{k}_{f}}\leftarrow {{v}_{i}}{{j}_{i}}{{k}_{i}}}^{J}$, the initial-state selected reaction cross sections are obtained according to

$$\begin{eqnarray}\begin{array}{ccccccccccccccc} {{\sigma }_{{{v}_{i}}{{j}_{i}}}} & = & \frac{\pi }{k_{{{v}_{i}}{{j}_{i}}}^{2}(2j+1)}\mathop{\sum }\limits_{J=0}^{{{J}_{{\rm max} }}}(2J+1)\mathop{\sum }\limits_{{{v}_{f}}{{j}_{f}}{{k}_{f}}{{k}_{i}}}|{\bf S}_{{{v}_{f}}{{j}_{f}}{{k}_{f}}\leftarrow {{v}_{i}}{{j}_{i}}{{k}_{i}}}^{J}({{E}_{c}}){{|}^{2}}, \\ \end{array}\end{eqnarray} \tag{ 1 }$$

where ${{k}_{{{v}_{i}}{{j}_{i}}}}=\sqrt{\frac{2\mu {{E}_{c}}}{{{\hbar }^{2}}}}$ is the wave vector in the incident channel and μ is the three-body reduced mass given by $\mu =\left[ \frac{{{m}_{{\rm Li}}}\times ({{m}_{{\rm H}}}+{{m}_{{\rm F}}})}{{{m}_{{\rm Li}}}+{{m}_{{\rm H}}}+{{m}_{{\rm F}}}} \right]$. The helicity quantum numbers, k_i and k_f for the initial and final states are restricted to $0\leqslant {{k}_{i}}\leqslant {\rm min} (J,{{j}_{i}})$ and $0\leqslant {{k}_{f}}\leqslant {\rm min} (J,{{j}_{f}})$. In the present study, the initial rotational quantum number of the HF molecule is restricted to j_i = 0 and only one triatomic parity P = (−1)^J needs to be considered for each J. For the heteronuclear HF molecule, the diatomic parity, q is irrelevant.

Extensive convergence studies of J = 0 initial-state-selected and state-to-state reaction probabilities have previously been performed by Weck and Balakrishnan [48] with respect to various parameters of the ABC code, including, maximum rotational quantum number ${{j}_{{\rm max} }}$, energy cut-off, emax, the maximum value of the hyperradius ${{\rho }_{{\rm max} }}$, and spacings between adjacent sectors in hyperradius, $\Delta \rho$. Converged results were obtained for ${{\rho }_{{\rm max} }}=50.0\;{{a}_{0}}$, j_max = 20, emax = 2.9 eV and $\Delta \rho =0.005\;{{a}_{0}}$ for the PES of Aguado et al [37]. The choice of emax = 2.9 eV and j_max = 20 in the basis sets led to 771 solutions for J = 0. Note that emax and j_max parameters in the ABC code do not discriminate reactant and product molecules and they apply to both reactant and product channels. This may lead to undesirably large number of coupled channel equations especially when the initial state of the reactant molecule is restricted to the ground state and/or the reactant and product molecules have widely different rotational and vibrational energy level spacings. Also, for $J\ne 0$ calculations, the number of coupled-channel equations will increase rapidly with k-values. Thus, we have found it convenient to introduce separate cut-off energies for the reactant and product molecules to keep the basis set optimum. In figure 1, we present results of convergence tests with respect to the energy cut-offs, emax1 (HF) and emax2 (LiF), respectively, for the reactant and product molecules for J = 0 (the LiH+F channel is not energetically accessible in the energy range considered here and is not included in the calculations). To enable this, we have appropriately modified the ABC code to accommodate the two cut-off energies. Figure 1 shows that emax1 = 2.25 eV and emax2 = 1.25 eV yield nearly the same results as emax = 2.9 eV adopted in [48]. This choice results in only 297 solutions for J = 0 compared to 771 in [48]. The computational savings are much more dramatic for higher J. All other parameters, j_max, ρ_max and $\Delta \rho$ are kept the same as in [48]. The convergence tests with respect to k_max are presented for J = 5 and 6 in figure 2 for the PES of Aguado et al [37]. The upper panel corresponds to J = 5 and the lower panel corresponds to J = 6. For both cases, significant deviations are found by restricting k_max = 4. However, k_max = 5, yields nearly the same results as k_max = 6 (lower panel) for J = 6. Hence for the production calculations with this potential for $J\geqslant 6$ we have restricted k_max = 5. This leads to a total of 1557 coupled-channel solutions for $J\geqslant 6$.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zanchet et al [17] reported state-to-state reaction probabilities and cross sections for the Li + HF $({{v}_{i}}=0,{{j}_{i}}=0)\to$ LiF + H reaction on the APW PES [38] using the time-dependent wave packet method. Figure 3 displays a comparison of reaction probability for this process from the present study and that of Zanchet et al [17] for J = 0. The black curve denotes the wave packet results and the various colored curves correspond to the present results. The agreement is excellent and it shows the convergence of our results with respect to ρ_max and $\Delta \rho$. However, similar comparisons for J = 10 presented in figure 4 show rather surprising trend with respect to k_max for the present results. The black curve denotes the wave packet calculations of Zanchet et al [17] whereas green, blue and red curves denote our results. They correspond to k_max = 5 (green), 7 (blue) and 10 (red), respectively. Note that k_max = 10 yields exact results for J = 10 in the ABC code. It is seen that only k_max = 10 yields good agreement with the wave packet results, indicating that, truncating the basis set with a smaller value of k_max leads to significant error in the reaction probabilities. Interestingly, Zanchet et al found that ${{k}_{{\rm max} }}$ (${{\Omega }_{{\rm max} }}$ in their notation) can be restricted to 7 in their wave packet calculations employing reactant Jacobi coordinates and they obtained nearly indistinguishable results for k_max = 7, 15 and 31 for J as high as 45. The observed convergence pattern on the APW PES is somewhat surprising considering the good convergence obtained for the PES of Aguado et al [37] for k_max = 5 for the J = 6 cross sections shown in figure 2. However, it is possible, even in this case, for higher values of J, restriction on k_max may lead to larger errors. Since the Li + HF reaction has a bent transition state, many values of the projection quantum number contribute to the reaction probability and a truncation of the basis set with a smaller ${{k}_{{\rm max} }}$ may lead to errors unlike reactions with linear geometry for the transition state. Although a value of ${{k}_{{\rm max} }}=7$ was found to be adequate for $J\gt 7$ in the wave packet calculations of Zanchet et al [17] to yield vibrationally resolved cross sections (a larger values of k_max = 15 was required for differential cross sections) this value may also depend on the coordinate system adopted (Jacobi versus hyperspherical) and the choice of body-fixed z-axis.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In light of the convergence pattern observed on the APW PES with respect to k_max, cross sections and rate coefficients reported here on the APW PES did not involve any truncation of the projection quantum number. Thus, the results reported on the APW PES are numerically exact for the basis set adopted. Since the maximum value of the projection quantum number is restricted by ${\rm min} (J,{{j}_{{\rm max} }})$, we have k_max = j_max for $J\gt 20$. The choice of j_max = 20 in the basis set led to a total of 3144 coupled-channel equations for $J\geqslant 20$. We have restricted j_max = 20 based on convergence tests for J = 0. It is possible that for higher J values, a larger basis set with j_max = 30 or higher may be required to accurately describe rotational populations of the heavy LiF molecule formed. However, choice of j_max = 30 will lead to a total of 5888 coupled-channel equations for $J\gt 30$ when no truncation in k_max is applied. Such calculations, are beyond the scope of this paper. Since our main goal is to investigate the dynamics of the reaction in the low energy regime, we believe, meaningful predictions of the low energy behavior of the reaction can still be made with the basis set adopted here. In addition to k_max, our convergence tests on the APW PES [38] with respect to various basis set parameters led to the choice of emax1 = 1.75 eV and emax2 = 1.25 eV. The somewhat smaller basis set for this case may be attributed to the slightly different exoergicity and diatomic potential energy functions. The number of coupled-channel equations specified above are based on these energy cutoff values.

In figure 5 we present total and partial-wave resolved reaction cross sections for Li + HF(${{v}_{i}}=0,{{j}_{i}}=0)\;\to \;$ LiF(v_f , j_f ) + H collisions for incident energies ranging from 10⁻⁷ to 0.4 eV. The partial cross sections for J = 0 − 7 (l_i = J for j_i = 0 initial state) are depicted by different colored curves. The total cross section that includes contributions from J = 0−40 is shown by the solid black curve. These results correspond to the PES of Aguado et al [37]. Note that the calculations based on this PES are numerically exact for energies up to 10⁻³ eV. However, the results may not be fully converged with respect to k_max (in this case k_max = 5, based on results from figure 2) in the higher energy range. The three different peaks labeled A, B and C, respectively, show three resonances originating from different partial waves. For peak A, the dominant contribution arises from J = 0 (red curve), and it persists in the total cross section. The other two resonances, B and C originate from l_i = 5 partial wave. Interestingly, shape resonance B centered at 2.5 × 10⁻⁵ eV, appears at a much lower energy compared to A.

Zoom In Zoom Out Reset image size

We have performed similar computations on the APW PES [38]. A comparison of reaction cross sections for LiF formation in Li + HF(${{v}_{i}}=0,{{j}_{i}}=0)\;\to \;$ LiF(v_f , j_f ) + H collisions evaluated using the APW PES (red curve) and the PES of [37] (black curve) is shown in figure 6. Both sets of results include contributions from J = 0 − 40. The results show dramatic sensitivity to the choice of the PES at energies below 1 meV. Results on the PES of Aguado et al [37] are significantly suppressed for energies in the range ${{10}^{-7}}\to {{10}^{-3}}$ eV compared to that of the APW PES. The difference is more than seven orders of magnitude for energies below 10⁻⁵ eV. This is attributed to the larger energy barrier for the PES of Aguado et al [37] which lead to less efficient tunneling as the collision energy is decreased. The smaller energy barrier and larger exoergicity for the APW PES also lead to diminished resonance structures in the energy dependence of the cross section compared to the PES of [37]. However, the discrepency between the two results tends to disappear as the collision energy is increased and for energies above 0.1 eV both PESs yield comparable results.

Zoom In Zoom Out Reset image size

Figure 7 displays the state-to-state reaction cross section for LiF(v_f ) formation in different vibrational levels as a function of the incident kinetic energy. The upper panel corresponds to results obtained from the PES of Aguado et al [37] whereas the results in the lower panel are based on the APW PES [38]. The different curves in both panels correspond to different vibrational populations of the LiF molecule. The dominant population of v_f = 0 vibrational level for both PESs indicates its wide accessibility in the entire energy range. Vibrational levels (v_f = 1 − 3) are populated as they become accessible at higher collision energies. Only v_f = 0 − 3 are accessible in the collision energy explored here for the PES of Aguado et al [37] whereas v_f = 0 − 4 are accessible on the APW PES. This is attributed to the higher exoergicity of the APW PES resulting from the choice of diatomic potentials for HF and LiF that differ from the PES of Aguado et al [37]. This is clearly visible in the threshold energies for the population of different vibrational levels of the LiF molecule from the two surfaces.

Zoom In Zoom Out Reset image size

A comparison of rotationally resolved cross sections between the two PESs is presented in figure 8 as a function of the rotational quantum number of the LiF molecule for several collisional energies between 10⁻⁷ and 0.1 eV. Only contributions from v_f = 0 are included. The left panel corresponds to calculations based on the PES of [37] while the right panel presents results on the APW PES [38]. Due to the smaller exoergicity, only rotational levels j_f = 0−7 are accessible on the PES of [37] for energies below 10⁻³ eV. The full spectrum of rotational levels (j_f = 0−20) is accessible at energies above 10⁻¹ eV. On the other hand, due to the higher exoergicity, rotational levels j_f = 0−20 are populated on the APW PES at all values of collision energies. The requirement of higher k_max to yield converged results for the APW PES, may in part, be attributed to the population of high rotational levels of LiF. On the other hand, the availability of fewer rotational levels (from ${{10}^{-7}}\to {{10}^{-2}}$ eV) on the PES of [37], indicates that a smaller value of k_max may be adopted for the basis set. Indeed, the rapid convergence with respect k_max for J = 6 shown in figure 2 is consistent with this result. At low energies the reaction appears to selectively populate a few rotational levels on the PES of [37], while a bimodal distribution is observed on the APW PES. At a collision energy of 0.1 eV, the populations approach a Boltzmann-type distribution on both PES.

Zoom In Zoom Out Reset image size

In the upper panel of figure 9, we compare reaction cross sections obtained from our work with the recent experimental results of Bobbenkamp et al [16]. The result on the PES of [37] is shown by the green curve whereas that on the APW PES is depicted by the red curve. The experimental results are denoted by various symbols (triangle up, triangle down, circles and squares) as reported by Bobbenkamp et al [16] which include results of their prior measurements [25]. Unfortunately, the experimental cross sections were reported in arbitrary units and for comparison purpose we have scaled the computed results by a factor to match the experimental data at a chosen collision energy. Since the computed cross sections are less oscillatory in the energy range of 0.1–0.4 eV, we chose an energy value of 0.32 eV for which experimental result is available. These factors are 3.43 and 1.32, respectively for the PES of [37] and APW PES [38]. The computed cross sections are scaled by this factor in the entire energy range shown in the upper panel of figure 9. It can be seen that the results are in overall good agreement with the experimental data except that the computed results show a series of resonances below 0.15 eV which are absent in the experiments. It should be mentioned that due to the elevated nozzle temperature the mean rotational level of the HF molecule in the experiment is 1.4 [16] while our results is for j_i = 0.

Zoom In Zoom Out Reset image size

The bottom panel of figure 9 shows a comparison of our results (green and red curves) with those of Zanchet et al [17] on the APW PES evaluated using a wave packet approach (blue curve). Interestingly, the results of Zanchet et al exhibit a different threshold behavior compared to the present results and also the experimental data. While we do not know the reason for the difference in threshold behavior of the cross sections reported by Zanchet et al [17], it is well known that the time-dependent wave packet calculations are less reliable at low energies and at energies close to the reaction threshold. However, the wave packet calculations and the present results (red curve) on the APW PES show good agreement in the 0.01–0.05 eV energy range. At higher energies, our calculations consistently underestimate the cross sections compared to the wave packet results. This may be attributed to the choice of j_max = 20 in our calculations compared to j_max = 31 (based on the value of Ω_max reported) by Zanchett et al [17]. Thus, it appears that for energies above 0.05 eV, the choice of ${{j}_{{\rm max} }}=20$ is not adequate to yield converged results. The larger difference between the results obtained from the PES of [37] (green curve) and the those of Zanchett et al may also be attributed to the smaller basis set and ${{k}_{{\rm max} }}=5$ restriction for the projection quantum numbers.

Finally, rate coefficients for the Li + HF(v = 0, j = 0) reaction derived from the two PESs are shown in figure 10 as a function of the temperature. The resonance structures are washed out during the Boltzmann average over the collision energies. We are not aware of any experimental measurements of thermal rate coefficients for this reaction. The rate coefficients on the PES of Aguado et al [37] exhibit a sharp threshold below 10 K reflecting the energy dependence of the corresponding reaction cross section. In contrast the rate coefficients on the APW PES show no threshold and continue to increase with decreasing collision energy. Indeed, the behavior of the rate coefficient on the APW PES is very similar to the experimental cross sections depicted in figure 9.

Zoom In Zoom Out Reset image size

The hyperspherical coupled-channel formalism adopted for the Li + HF system is computationally demanding and not practical for the vast majority of systems of interest to current cooling and trapping experiments. As a result, over the last several years, a number of approximate methods based on long-range theories and the MQDT formalisms have been developed to estimate the total depletion rate coefficients for ultracold chemical reactions [55–58]. These formalisms provide simple expressions for the overall rate coefficients and require only the knowledge of the analytical behavior of the long-range interaction potential. The short-range interaction potential, which is important for describing chemical reaction, is not explicitly considered in these models and they do not yield product resolved reaction rate coefficients [55, 57]. To address this issue, we have recently formulated an alternate approach that combines quantum CC with the MQDT formalism. The CC approach is based on the formalism given in the ABC code [53] and the MQDT formalism is based on the approach developed by Ruzic et al [59]. The CC calculations are restricted to the region of the strong interaction potential where chemical transformation occurs and the MQDT formalism is used outside this region. Matching the solutions of the CC equations to the MQDT reference functions at a suitable matching distance outside the region of reactive coupling yields a short-range K-matrix, ${{K}^{{\rm sr}}}$, that is weakly energy-dependent in the threshold regime. Details of the formalism, including, transformation of the short-range K-matrix to physical K and S matrices are given in [18, 60]. The method yields full product state-resolved cross sections as in full CC calculations but at a much reduced computational effort as demonstrated for D + ${{{\rm H}}_{2}}(v,j)\to$ HD(v', j') + H reaction [18].

While it would be desirable to apply this approach for the Li + HF reaction and compare the results with the full CC calculations presented here we have not been able to implement it for the PESs adopted in this study as no analytical long-range form is included in either of the two PESs. Construction of a full analytical long-range interaction that will smoothly fit to the short-range interaction potential is beyond the scope of this work. Instead, we highlight usefulness of this approach by presenting illustrative results for the benchmark D + H₂ and H + D₂ reactions focusing on the initial vibrational level v_i = 3 for both molecules. In left panels ((1a) and (1b)) of figure 11, we show the total and vibrational-level resolved reaction cross sections for D +${{{\rm H}}_{2}}({{v}_{i}}=3,{{j}_{i}}=0)$ → H + HD $({{v}_{f}},{{j}_{f}})$ reaction for J = 0. The total cross section for this vibrational level has been presented in Hazra et al [18] but not the product resolved cross section. Similar results are presented in the right panels ((2a) and (2b)) for H + ${{{\rm D}}_{2}}({{v}_{i}}=3,{{j}_{i}}=0)\to$ D + HD $({{v}_{f}},{{j}_{f}})$ reaction. For ${{v}_{i}}=3,{{j}_{i}}=0$ initial state, the number of coupled-channels is restricted by E_max = 4.25 eV (this includes v_max = 10 for H₂, v_max = 15 for D₂ and v_max = 12 for HD in the basis sets) and j_max = 8 for rotational levels. This leads to a total of 165 and 189 coupled-channels for D + H₂ and H + D₂ reactions, respectively. In the figure the solid black curves represent the results obtained from full CC calculations whereas the results from the CC–MQDT approach are shown by the red curves. The full CC results are evaluated by matching the log-derivative matrix with the free particle solutions at 100 ${{a}_{0}}$ whereas a short range matching distance of 15 ${{a}_{0}}$ has been chosen for the MQDT part. The agreement between full CC and MQDT methods is excellent and for both reactions the percentage error is less than 5%. Instead of performing the calculation at each energy grid which is essential in full CC, an interpolation of the short-range K-matrix is employed. This results in a significant saving in computation time. Due to the weak energy dependence of the ${{K}^{{\rm sr}}}$ matrix in the Wigner threshold regime, it is only evaluated at 1 μK and 1 mK and the same is used at all energies in this region. Between 100 mK and 1 K a few more ${{K}^{{\rm sr}}}$ matrices are evaluated with an energy spacing of 200 mK. Thus, the entire energy region shown in figure 11 includes actual CC calculations only at seven energies for the CC–MQDT approach and ${{K}^{{\rm sr}}}$ matrix at other energies in this regime is obtained from an interpolation scheme. Figure 11 shows that for D + H₂, vibrational levels of HD up to v_f = 3 are populated ((1b) of left panel) whereas for H + D₂, only v_f = 0 − 2 of HD are accessible ((2b) of right panel). This is due to the difference in the zero-point energies of H₂ and D₂ that leads to different exoergicities for the reaction starting from the same initial vibrational level.

Zoom In Zoom Out Reset image size

A comparison between full CC and CC–MQDT is presented for rotationally resolved cross sections in figure 12 for D + H₂ and figure 13 for H + D₂. As before the black curves correspond to full CC and the red curves correspond to CC–MQDT results. The different panels depict rotational-level resolved cross sections for different vibrational levels of the HD molecule. For the most part, the agreement is excellent as in the case of vibrational-resolved and total cross sections. Small deviations observed for a few highest populated rotational levels of the least populated HD(v_f = 2) state may be attributed to neglect of the anisotropy of the interaction potential in the construction of MQDT reference functions and parameters.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In this paper a time-independent quantum mechanical investigation of the Li + HF reaction is presented focusing on its dynamics at low collision energies to the thermal energy regime. The computations were performed using the ABC code [53] and by introducing separate energy cutoffs for reactant and product rovibrational energy levels a significant reduction of the basis set size is achieved without compromising on the accuracy of numerical results. A key motivation for this study is the recent experimental work of Bobbenkamp et al [16] that revealed a vanishing threshold for the reaction compared to a recent quantum mechanical study based on the time-dependent wave packet method. Our results are consistent with the experimental data indicating a vanishing threshold. Two PESs for the Li + HF reaction have been adopted in our studies to explore sensitivity of the reaction dynamics to details of the interaction potential. Both depict numerous narrow resonances supported by the van deer Waals interaction potential. However, results on the two PESs differ dramatically at energies below 1 meV where the PES of Aguado et al [37] yields cross sections that are three to seven orders of magnitude smaller than that obtained from the presumably more accurate APW PES. Results illustrate high sensitivity of the reaction to slightly different barrier heights and depths of the van der Waals interaction potentials for the two PESs. The reaction occurs primarily by tunneling at energies less than 100 meV and our results illustrate the sensitivity of tunneling dominated-reactions to small changes in the barrier height. Temperature dependence of thermal rate coefficients evaluated using the APW PES seems to reflect the energy dependence of the cross section measured by Bobbenkamp et al.

We have also demonstrated the usefulness of a recently introduced method for cold and ultracold chemical reactions that combines quantum CC at short-range with the MQDT method at long-range. The success of the method is illustrated for rotationally and vibrationally resolved cross sections for D + H₂ and H + D₂ reactions with three quanta of vibrational excitation for both molecules.

This work was supported in part by NSF grant PHY-1205838 and ARO MURI grant No. W911NF-12-1-0476. Authors are grateful to Prof John Bohn and Brandon Ruzic for collaboration on the development of the CC–MQDT method. NB is grateful to Profs Marcel Mudrich and Hansjürgen Loesch for providing the experimental data on the Li + HF reaction and Prof Octavio Roncero for helpful discussions and for providing the results of quantum wave packet calculations reported in [17].