Kinetic Isotope Effect in Low-Energy Collisions between Hydrogen Isotopologues and Metastable Helium Atoms: Theoretical Calculations Including the Vibrational Excitation of the Molecule

We present very accurate theoretical results of Penning ionization rate coefficients of the excited metastable helium atoms (4He(23S) and 3He(23S)) colliding with the hydrogen isotopologues (H2, HD, D2) in the ground and first excited rotational and vibrational states at subkelvin regime. The calculations are performed using the current best ab initio interaction energy surface, which takes into account the nonrigidity effects of the molecule. The results confirm a recently observed substantial quantum kinetic isotope effect (Nat. Chem. 2014, 6, 332–335) and reveal that the change of the rotational or vibrational state of the molecule can strongly enhance or suppress the reaction. Moreover, we demonstrate the mechanism of the appearance and disappearance of resonances in Penning ionization. The additional model computations, with the morphed interaction energy surface and mass, give better insight into the behavior of the resonances and thereby the reaction dynamics under study. Our theoretical findings are compared with all available measurements, and comprehensive data for prospective experiments are provided.


INTRODUCTION
A chemical reaction is a fundamental process in chemistry. External control of its dynamics is crucial to understand many physical and chemical phenomena, especially at the quantum level. Recent sophisticated merged beams experiments for lowenergy atom−molecule collisions revealed the possibility of controlling molecular reactions. 1−10 These experiments were groundbreaking achievements in cold chemistry. 11−15 The first merged neutral supersonic beams technique was demonstrated in 2012 by the group of Narevicius, 1 approaching a collision temperature below 10 mK. Narevicius and co-workers investigated the Penning ionization (PI) reactions by colliding excited metastable helium atoms 4 He* (≡ 4 He(2 3 S)), contained in one beam, with ground-state molecular hydrogen H 2 , contained in another beam He H e 2 2 They observed a few shape resonances by measuring reaction rate coefficients, systematically decreasing collision energy from a few tens of kelvins down to millikelvins. In follow-up studies, they focused on the isotope effect, 2 where the hydrogen molecule was substituted by hydrogen deuteride (HD) and deuterium (D 2 ). This was the first observation of such an effect in low-energy collisions. It turned out that the positions of the resonances were different for the hydrogen isotopologues, far from intuitive expectations. The substitution strongly affected the PI process with an order of magnitude increase in the reaction rate coefficient at certain collision energies. This effect may significantly matter to the isotopic composition and formation of interstellar clouds and structure. 16,17 Hydrogen and helium are the lightest and most abundant chemical elements in the Universe, while molecular hydrogen is the lightest and the most common molecule. 18,19 For instance, molecular hydrogen and helium are the prime constituents of gas giants. 20 New ground-and space-based high-resolution telescopes enable the observation of light elements, which was previously inaccessible. Very recently, scientists have made the first-ever detection of 2 3 S helium in the atmosphere of one of the planets (WASP-107b) beyond the solar system. 19 Soon after, excited metastable helium was confirmed on other exoplanets. 21−29 The measurements of the helium absorption signal in exoplanet atmospheres are now important markers of a variety of planetary processes. They have become a unique probe for assessing exoplanetary compositions and understanding the formation of exoplanetary systems. The observation of long-lived metastable helium atoms has changed the way astronomers look at exoplanets and has raised the significance of theoretical chemistry for astrophysical domains. Undoubtedly, the collisions of 2 3 S helium with various atomic and molecular partners, especially with hydrogen isotopologues, are of primary importance for astrochemistry. It is also worth noting that there are a considerable number of experimental studies for He* in collisions with more complex molecules (for example, hydrogen sulfide, 30 nitrous oxide, 31 benzene, 32 naphthalene, 33 anthracene, 33 [2,2]-paracyclophane, 34 pyrene, 35 chrysene, 35 coronene, 35 glycine, 36 anisole, 37 thioanisole, 37 and selenoanisole 37 ).
As was shown in ref 2, the kinetic isotope effect in the collisions between light subsystems (such as 4 He* with H 2 , HD, and D 2 ) is dramatic, fundamentally changing the reactivity of the PI process. A high level of theory is required to evaluate the experimental findings. The main motivation of our present studies was to confirm the strong isotope effect in cold chemistry reactions based on our recently developed, very accurate interaction energy surface. 38 This ab initio surface includes nonrigidity effects of the molecule, which play an important role in low-energy molecular anisotropic collisions. 38 Additionally, this surface also allows us to predict resonance structures for experimentally uncharted cases. Since we can account for nonrigidity effects, complexes with the vibrationally excited hydrogen molecule or its isotopologues can be considered. Such complexes have become easier to access experimentally due to recent progress in the production and the control of the excited hydrogen isotopologues. 39−41 Thus, in the present work, we compute reaction rate coefficients of the H 2 /HD/D 2 molecules in the ground and excited rotational (j = 0/1) and vibrational (v = 0/1) states with 4 He*, and additionally with 3 He*. We demonstrate that by changing the quantum state of the molecule, one can turn on or off particular low-energy resonances, thereby strongly modify the PI. It enables us better control over reactions in the subkelvin regime.

THEORY
2.1. Interaction Energy Surface. Our theoretical approach consists of two parts: the first is to prepare the best possible interaction energy surface, and the second is to perform quantum-dynamical calculations of the PI rate coefficients. To parameterize the He* + H 2 complex, we use Jacobi coordinates: the distance R between He* and center of mass (COM) of H 2 , the angle θ between the H 2 bond and the COM−He* direction, and the distance r between the hydrogen nuclei. Thus, the full-dimensional energy surface can be written as a function of these coordinates, V(R, θ, r). In our previous paper, 38 we have shown for 4 He* + H 2 that if the widely used rigid-rotor approximation is applied, one cannot achieve a quantitative agreement with experiment, even if stateof-the-art ab initio methods are employed to calculate the interaction energy surface. The rigid-rotor approximation procedure involves two stages: the construction of the interaction energy surface for rigid molecules, and scattering calculations limited to only intermolecular coordinates. We have demonstrated 38 that if the interaction energy surface averaged over the vibrations of the hydrogen molecule is used instead of the rigid-molecule surface, an excellent agreement with the experiment can be obtained, despite the rigid-rotor approximation in the scattering calculations. This is because though the vibrationally averaged surface depends only on the intermolecular coordinates, to calculate it one needs to know the dependence of the interaction energy also on the intramolecular coordinates. Thus, the surface includes effective information of the nonrigidity effects and is not equivalent to the rigid-monomer one. In principle, to calculate the vibrationally averaged surface, one needs to know a fulldimensional surface for a given complex, 42 however, the surface averaged over the vibrational state v can be obtained approximated 43 using a Taylor expansion of the interaction potential V(R, θ, r) with respect to the intramonomer coordinate r, around some reference value r c where f 0 (R, θ) = V(R, θ, r c ), f 1 (R, θ) = ∂V(R, θ, r)/∂r at r = r c , and f 2 (R, θ) = ∂ 2 V(R, θ, r)/∂r 2 at r = r c . The ⟨r⟩ j,v and ⟨r 2 ⟩ j,v symbols denote the averaged values of r and r 2 for the rovibrational states of H 2 defined by the quantum numbers j and v. In our investigations, we took r c equal to the H−H distance averaged over the rovibrational ground state, i.e., r c = ⟨r⟩ 0,0 = 1.4487 Bohr. The values of f 1 and f 2 can be calculated on the grid of the intermolecular coordinates with the finitedifference method. The vibrationally averaged surfaces obtained within this approximation were used in bound-state calculations for the H 2 −CO, 44−46 Ar−HF, 47 and N 2 −HF 48 complexes, and an excellent, or at least, very good agreement of the theoretical spectra with experimental ones was observed.
For the H 2 −CO complex, scattering calculations were also performed 49 with the surface defined through eq 1, and an excellent agreement with the experimental data 50,51 and the full-dimensional results was achieved.
In the present study, we have used the values of the f i expansion functions, where i = 0, 1, 2, calculated in ref 38. The details can be found in that paper, but here we would like to mention that the leading f 0 term was obtained using the coupled cluster method with singles, doubles, and perturbative triples (CCSD(T)), and improved with a full configuration interaction (FCI) correction. The f 1 and f 2 derivatives were calculated with the symmetry-adapted perturbation theory (SAPT) method. 53 As was demonstrated for the 4 He* + H 2 complex, the diagonal Born−Oppenheimer (DBO) correction 54 is significantly smaller than other contributions to f 0 , so we neglected this correction in that work. Moreover, for the 4 He* + HD and 4 He* + D 2 complexes, the DBO correction can be expected to be even smaller than for 4 He* + H 2 due to the larger masses of HD and D 2 in comparison to H 2 . Thus, we can still neglect the DBO correction, and the f i functions therefore do not depend on masses of helium or hydrogen isotopologues.
According to eq 1, to obtain the vibrationally averaged surface ⟨V⟩ j,v (R, θ), we need to know the values of ⟨r⟩ j,v and ⟨r 2 ⟩ j,v for H 2 , HD, and D 2 , for the required vibrational states v of these molecules. Moreover, for the ortho-H 2 and para-D 2 cases, when the molecules are in the j = 1 rotational state, the averaging was done with the (j = 1, v) rovibrational wave functions of H 2 and D 2 , respectively. Although it was shown in ref 49 for H 2 −CO that the cross sections calculated for the (j = 1, v) averaging are only slightly different from those obtained with the (j = 0, v) one, we have used the appropriate versions of ⟨r⟩ j,v and ⟨r 2 ⟩ j,v depending on the system, for the sake of the overall accuracy. The whole set of values of ⟨r⟩ j,v and ⟨r 2 ⟩ j,v are given in Table 1.
To demonstrate how the vibrational averaging affects the interaction energy surface for various isotopologues and vibrational states of H 2 , we present, in Figure 1, a comparison of the interaction energy surfaces ⟨V⟩ j,v (R, θ), defined in eq 1, Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article for two angular orientations, θ = 0 and π/2, as a function of the intermolecular separation R. The corresponding numerical data are provided in the Supporting Information. For v = 0, the curves representing different isotopologues are rather close to each other. However for v = 1, the differences are much larger; for example, the minimum of the interaction energy for He* + H 2 is equal to −18.1 cm −1 , compared to −16.7 cm −1 for the isotopically substituted He* + D 2 . Even larger differences are due to the vibrational excitation of the molecules; for example, the minimum energy for He* + H 2 (v = 0) is equal to −15.1 cm −1 , which is 3 cm −1 higher than the v = 1 minimum, a 20% change. The position of the minimum is also shifted toward smaller values of R. Such large differences in the potentials used in the scattering calculations will be seen to result in substantial differences in the values of the rate coefficients to be discussed later.

Surface Expansion.
The intermolecular potential given in eq 1 may be expanded in Legendre polynomials as follows 38 The first term, V 0 (R), is isotropic. The odd terms vanish for symmetric molecules, thereby the first anisotropic term for He* colliding with H 2 and D 2 is of the form V 2 (R)(3 cos 2 θ − 1)/2. However, in the case of He* + HD, the center of mass of HD is not in its geometric center, and if the above expansion is carried out in the appropriate Jacobi coordinates, the V 1 term does not vanish. 55,56 In Figure 2, we display the isotropic, V 0 , and anisotropic, V 1 and V 2 , radial interaction potential terms of He* + H 2 /HD/D 2 .
As one could expect, the leading term V 0 is almost 3 times deeper than the corresponding V 2 . The V 1 term, present for the He* + HD case, has a magnitude smaller than V 2 for R > 9.45 Bohr and larger for shorter intermolecular distances. This means that the repulsive wall of the V 0 potential is modified more by V 1 than by V 2 . Our main interest is the dependence of the V η terms on the masses of the molecules and their vibrational states. One can see that V 0 becomes shallower when the mass of the molecule increases. There is no such simple correlation of the depth of V 2 with the mass of the molecule. Figure 2 presents also the radial potentials when the molecule in the complex is excited to the first vibrational v = 1 state. The isotropic term is significantly shifted down compared to the v = 0 state. It should be emphasized that the anisotropy here is doubtless smaller compared to the rotational constant (B) of the colliding molecule.
2.3. Adiabatic Variational Theory (AVT) and Adiabats. The widely used method for quantum scattering calculations is close coupling. This method requires a properly prepared interaction energy surface. For instance, Balakrishnan and coworkers 57 performed close-coupling computations for cold collisions of nonexcited 4 He(1 1 S) and 3 He(1 1 S) with H 2 including vibrations of the molecule. However, in our studies, we have applied the recently developed adiabatic variational theory (AVT) for cold atom−diatom collision experiments. 58−60 It involves the adiabatic approximation and holds for any collision where the anisotropy of the interaction energy is small compared to the molecular rotational constant. In general, this theory allows us to reduce the multidimensional problem to a set of more easily solvable subproblems, without losing physical essence, by separating the radial from the angular motion. Thus, the interaction potential has to be expanded in terms of products of radial and angular functions, as, for example, in eq 2. Then, taking more than one nonzero term of the expansion guarantees that the anisotropy is included. AVT works well with complex potentials, where the imaginary part is related to the exponential decay rate of the entire system being in a resonance state. 61−63 In our    Within the framework of AVT, we construct the potential in eq 2 with the rotation terms (l̂2/(2μ A+M R 2 ) and Bĵ2) in the matrix representation in the basis set of products of two spherical harmonics. Here, μ A+M is the reduced mass of the complex, l̂2 is the operator for the end-over-end rotation of the complex, and ĵ2 stands for the angular momentum operator of the molecule. Each matrix built for a given R (where R is a parameter) has to be diagonalized. By repeating this procedure for different values of R and ordering the eigenvalues, we obtain one-dimensional effective adiabatic potentials, known as adiabats. The technical details are described in ref 59. The adiabats asymptotically coincide with a specific rotational state of the molecule (j) as well as an end-over-end angular momentum referred to as a partial wave (l). Finally, we solve the radial Schrodinger equation separately for every adiabat. The eigensolutions serve to determine reaction rate coefficients as is shown in ref 60.

RESULTS AND DISCUSSION
3.1. Reaction Rate Coefficients. We have computed the PI rate coefficients for the reactions between the hydrogen isotopologues H 2 , HD, and D 2 , and the excited 4 He atoms in the metastable 2 3 S state, at the wide collision temperature range, from several dozen of kelvins down to 1 mK. In a supersonic expansion, the molecular beam is cooled to the lowest rotational state with the para/ortho ratio of 1/3 for H 2 and 1/2 for D 2 . Consequently, the only two states that matter are j = 0 for para-H 2 and ortho-D 2 , and j = 1 for ortho-H 2 and para-D 2 . Note that the HD molecule does not exhibit spin isomerism. Therefore, we first performed calculations for the vibrational ground-state hydrogen and deuterium molecules in the j = 0 and 1 rotational states, and for the rovibrational ground-state hydrogen deuteride (the j = 0 state only). These results are presented in the left panels (a 0 , b 0 , and c 0 ) of Figure  3 (see the solid blue (j = 0, v = 0) and red (j = 1, v = 0) curves). For comparison, we added all available experimental findings. Black dots with error bars represent the first-ever experimental results for H 2 (mixture para-H 2 and ortho-H 2 ; the a 0 panel), HD (the b 0 panel), and D 2 (mixture ortho-D 2 and para-D 2 ; the c 0 panel). 2 One can see that the positions of the resonances agree very well with our theoretical findings, but there is a clear difference in the relative intensities. Recently, the experiment was improved and the para-H 2 sample was isolated 6 from normal-H 2 . Hence, Figure 3 also contains the reaction rate coefficients separately for 4 He* + para-H 2 and 4 He* + ortho-H 2 (purple and green dots with error bars in the top left panel, respectively). The agreement of these results with the theoretical predictions is remarkable as was discussed for this particular case in detail in our last work. 38 Unfortunately, improved PI rate coefficients for 4 He* + ortho-D 2 , and indirectly for 4 He* + para-D 2 , are not available. Figure 3 clearly confirms the strong quantum kinetic isotope effect, initially discussed in ref 2. Comparing the a 0 , b 0 , and c 0 panels of Figure 3, one can see that the positions and intensities of the peaks are different. For instance, in the case of 4 He* + D 2 (the c 0 panel), there is no strong peak below 1 K. When we perform the isotopic substitution (H 2 → HD → D 2 ) in the theoretical calculations, there are two factors that have to be changed: the interaction energy surface and the mass of the molecule. The vibrationally averaged interaction energy becomes progressively weaker in the sequence H 2 → HD → → ⟨V⟩ j=0,v=0 (He*+ HD) , as follows where λ is the morphing parameter that varies from 0 to 1 in steps of 0.1. The Jacobi coordinate system is adjusted to be consistent with the change of m morph . Figure 4 presents the morphed results of reaction rate coefficients, which demonstrate that the resonances move to the cold regime region when the H 2 → HD change is made. Thus, one can conclude that the effect of mass replacing in scattering calculations is stronger than the effect of changes of the potential caused by the molecule flexibility. A separate problem of the disappearance of the strongest resonances for low temperatures will be discussed in the next subsection. A similar analysis was done by one of us in ref 66 for a shape resonance in N + NH collisions. Generally, scaling the mass and the interaction energy surface is a valuable technique for understanding cold collision experiments and computations. In our studies, we have also determined the PI rate coefficients of metastable helium colliding with the hydrogen isotopologues in the first excited vibrational state (v = 1). These unique results are presented as dotted-dashed curves in Figure 3. One can see that by exciting the molecule from v = 0 to 1, the resonant structure is dramatically changed. The lowenergy resonance at 0.3 K for 4 He* + ortho-H 2 (j = 1, v = 0) vanishes when the orthohydrogen is excited to the v = 1 state. A similar situation occurs for 4 He* + HD(j = 0, v = 0 → v = 1). In turn, the higher energy resonances in all cases are shifted to the lower temperature range, simultaneously increasing their intensities. A mechanism of such behavior of the resonances is explained in Section 3.2.
The right panels of Figure 3 show the reaction rate coefficients, where the 4 He* atom is substituted by its stable lighter isotope, 3 He*. All of these results are significantly different in comparison to their counterparts from the left panels of Figure 3. Interestingly, for the 3 He* + ortho-H 2 (j = 1, v = 0) and 3 He* + para-D 2 (j = 1, v = 0) systems (compare the red curves in the d 0 and f 0 panels), the high-intensity resonances are accompanied by the much weaker ones clearly seen on their shoulders at the collision energy of k B × 0.8 K and k B × 1 K, respectively. The structure of the PI rate coefficients for the 3 He* + ortho/para-D 2 (j = 0/1, v = 1) collision pairs (the f 1 panel) attracted our attention, because it corresponds to the exceptional case that takes place for 4 He* + para/ortho-H 2 (j = 0/1, v = 0), which was widely discussed in ref 6. Let us shortly describe that issue. When the 2 3 S helium atom collides with the molecular hydrogen, the extra strong resonance appears or disappears depending on the rotational state of the molecule (j = 1 or 0) (see the solid red and blue curves, respectively, in the top left panel of Figure 3). The origin of this low-energy resonance at 0.3 K is entirely due to the anisotropy of the interaction between 4 He* and H 2 . 59 When the hydrogen molecule is in the ground rotational state, its wave function is spherically symmetric and only the isotropic component of the potential is probed. Within the present work, we found 3 He* + D 2 (v = 1) as another example where the anisotropy in atom−molecule interactions can be directly switched on/off by the manipulation of the rotational state of the molecule, drastically affecting the PI reaction at subkelvin temperatures. Now, we analyze the 3 He* + ortho/para-D 2 (j = 0/1) collisions extended in context of the vibrational state. In the f 0 and f 1 panels of Figure 3, we see two strong low-energy resonances for v = 0 (each for different j), whereas only one for v = 1. It means that by the vibrational excitation of the D 2 molecule in the lowest rotational state, the resonance vanishes. To better examine the resonances, first, we plotted in Figure 5 results without the convolutionjust sharp peaks. The plot reveals that the low-energy and high-intensity convoluted resonance for j = 1 and v = 0 consists, in fact, of two closely lying resonances at about 0.5 K. The third weaker peak at 1 K corresponds, after the convolution, to the barely visible one, on the shoulder of the largest one (see the f 0 panel of Figure 3). We would like to figure out the origin of these resonances. Hence, we focus on adiabats supporting bound and resonance states in the spirit of ref 6.
3.2. Analysis Based on Adiabats. Figure 6 shows the adiabats that support the resonances discussed above and presented in Figure 5, i.e., for 3 He* + D 2 . The adiabat in blue corresponds asymptotically to the rotational state j = 0, the partial wave l = 4, and the total angular momentum J = 4. The remaining three adiabats correlate asymptotically to the rotational state j = 1, the partial wave l = 4, and the total angular momenta J = 3, 4, 5 (yellow, red, and black curves, respectively). One can see how substantially the anisotropy  Figure 3] resonances, whereas the rest of the curves are to show the evolution from the former case into the latter one. The theoretical results have been convoluted with the experimental energy resolution.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article affects the effective potential. Each of the considered adiabats for v = 0 (the top panel of Figure 6) supports one resonance state whose energy level is between the dissociation threshold and the top of the barrier. When the molecule is excited to the first vibrational state (v = 1; the bottom panel of Figure 6), the quantum wells become deeper and then the energy levels are shifted down, three of them below the dissociation threshold.
Only the state with the highest energy (in red) still remains above the threshold. In other words, three resonances turned into bound states. Consequently, the resonance observed at the collision energy of k B × 0.7 K for ortho-D 2 (j = 0, v = 0), which is shown in blue in the top panel of Figure 5, disappears in the bottom one, whereas from the three resonances placed in the range of collision energies k B × (0.4, 1.1) K, corresponding to the para-D 2 (j = 1, v = 0) case, only one survives and is present at k B × 0.4 K in the bottom panel of Figure 5. Now, one can see that the latter resonance has the same quantum nature as the much weaker one observed as a feature on the shoulder of the strongest peak of the para-D 2 (j = 1, v = 0) case, in the f 0 panel of Figure 3.
To gain even better insight into a mechanism leading to the effects described above, we can perform the following numerical test. Let us define a surface where λ is again the morphing parameter, while ⟨V⟩ j=1,v=0 and ⟨V⟩ j=1,v=1 are the surfaces averaged over the (j = 1, v = 0) and (j = 1, v = 1) rovibrational states of D 2 , respectively. Changing the values of λ from 0 to 1, one can morph the interaction energy from the v = 0 form to the v = 1 one. Note that this time the mass of the molecule is not morphed. We have carried out calculations of the rate coefficients for the Ṽm orph potential gradually increasing λ in steps of 0.1. The results are presented in Figure 7, and we can see that the convoluted resonance for 3 He* + para-D 2 (j = 1, v = 0) at 1 K becomes the resonance at 0.4 K for v = 1, whereas the convoluted low-energy resonance at 0.5 K, which in fact consists of two closely lying resonances, shifts toward lower temperatures, and finally vanishes. A similar analysis (with adiabats and morphing) can be done for every v = 0 and 1 pair of panels presented in Figure 3 as well as for every pair of the j = 0 and 1 states considered in this work, since the mechanism of appearance and disappearance of resonances is the same. Moreover, it can be also performed for the isotopic substitution.

CONCLUSIONS
This paper shows the very accurate theoretical results of PI rate coefficients for all possible combinations of 4 He and 3 He in the metastable 2 3 S state colliding with the hydrogen isotopologues (H 2 , HD, D 2 ) in the two lowest rotational and vibrational states at subkelvin temperatures. For all but the 4 He* + para/ ortho-H 2 (v = 0) cases, these are the first theoretical results of such high quality. The calculations have been carried out based on the best available ab initio interaction energy surface, developed by us, 38 which includes molecule nonrigidity effects. The results confirm the recently observed strong quantum kinetic isotope effect in cold collisions. 2 The rate of the PI process can differ by as much as an order of magnitude at certain collision energies depending on the isotope or isotopologue. We have demonstrated that the existence of resonances is determined by the position of energy levels in quantum wells formed by adiabats. Consequently, the location of major peaks in the reaction rate coefficients can be controlled by the rovibrational state of the molecule. In sum, the change of the rotational and vibrational quantum numbers of the molecule can strongly enhance or suppress the PI reaction like a quantum switch. Moreover, we have presented the origin of the resonances and the mechanism of their disappearance. This was supported by model calculations with the morphed mass and interaction potential. We are positive that our findings are an important contribution to a precise control of chemical reactions in cold chemistry. Since nowadays the vibrational excitation of the hydrogen isotopologues is feasible, 39 Values of the interaction potential of He* + H 2 /HD/D 2 for two angular orientations: the linear one (θ = 0) and the T-shaped one (θ = π/2) (PDF)