Evidence of Nonrigidity Effects in the Description of Low-Energy Anisotropic Molecular Collisions of Hydrogen Molecules with Excited Metastable Helium Atoms

Cold collisions serve as a sensitive probe of the interaction potential. In the recent study of Klein et al. (Nature Phys.2017, 13, 35–38), the one-parameter scaling of the interaction potential was necessary to obtain agreement between theoretical and observed patterns of the orbiting resonances for excited metastable helium atoms colliding with hydrogen molecules. Here, we show that the effect of nonrigidity of the H2 molecule on the resonant structure, absent in the previous study, is critical to predict the correct positions of the resonances in that case. We have complemented the theoretical description of the interaction potential and revised reaction rate coefficients by proper inclusion of the flexibility of the molecule. The calculated reaction rate coefficients are in remarkable agreement with the experimental data without empirical adjustment of the interaction potential. We have shown that even state-of-the-art calculations of the interaction energy cannot ensure agreement with the experiment if such an important physical effect as flexibility of the interacting molecule is neglected. Our findings about the significance of the nonrigidity effects can be especially crucial in cold chemistry, where the quantum nature of molecules is pronounced.


INTRODUCTION
The breakthrough in controlling the movement and internal degrees of molecules with external fields, 1 which started about 20 years ago, currently allows one to study the collisions, reactivity, and properties of molecules in cold regime. Such molecules provide new and propitious prospects in precision spectroscopy, fundamental physics, astrochemistry, and quantum engineering. Understanding of quantum effects, resonance phenomena, and reaction dynamics in a low-energy range opens the gate to design and create materials with unusual functionality and elements of quantum computers. 2,3 Cold collision experiments realized by merging two supersonic beams have become an important technique for studying chemical reactions in temperatures near 1 K 4−9 in which unseen earlier quantum features, such as resonances or interference, are revealed. The group of Narevicius performed the first experiment in that field focused mainly on the Penning ionization (PI) process of colliding hydrogen isotopologues with excited metastable helium atoms. 4,5,8,9 In such a reaction, an electron is moved from the molecule to the only partially occupied orbital of an excited atom; then, the initially excited electron of the atom is kicked out of the system. Three products occur: the atom in the ground state, the molecular ion in the ground state, and a free electron. These colliding systems are of great interest to astrophysicists and astrochemists studying conditions and reactions in outer space. Hydrogen and helium are the most abundant elements, whereas molecular hydrogen is the most common molecular species in the Universe. 8,10 The first observation of a metastable helium atom in the atmosphere of one of exoplanets 11 has boosted its importance for astrochemistry, and one can expect that its interaction with the omnipresent hydrogen molecule will be carefully investigated.
Recently, Narevicius and coworkers directly probed the anisotropy in atom−molecule interactions through orbiting resonances by changing the rotational state of the molecule. 9 That work reveals a crucial role of the anisotropy of the interaction energy, due to various orientations of the H 2 molecule in the complex of He(1s2s, 3 S 1 ) (≡ He*) with H 2 , in the dynamics in the subkelvin regime. To elucidate physical phenomena presented in their novel experiment, the authors of ref 9 used state-of-the-art first-principles calculations. The interaction energy of the complex composed of the metastable helium atom and the hydrogen molecule was first calculated using the supermolecular approach employing the coupled cluster method with singles, doubles, and perturbative triples (CCSD(T)) and is further denoted as E int CCSD(T) . The CCSD(T) method is known as the gold standard in quantum chemistry and in many cases provides values of properties accurate enough to predict experimental results. 12 However, in the case of He* + H 2 , the values of reaction rate coefficients calculated with the potential based on E int CCSD(T) did not agree with the experiment even qualitatively. Thus, the correction to the interaction energy, denoted here by δE int FCI , was calculated from the full configuration interaction (FCI) method and added to the E int CCSD(T) energy. The agreement of the rate coefficients calculated from the E int CCSD(T) + δE int FCI energies with the experimental values was improved but still qualitatively incorrect, since one additional resonance, not present in the experiment, was predicted for low collision energies. To obtain quantitative agreement of the calculated rate coefficients with the experiment, the authors scaled the correlation part of the interaction energy by a factor of 1.004, which can be viewed as adding the correction 0.4%E int corr to the E int CCSD(T) + δE int FCI energy, suggesting that the basis set incompleteness was responsible for the inaccuracy of the rigorous ab initio interaction energy surface.
The motivation of the present studies was to find a reason why the ab initio interaction energy surface, even obtained at the FCI level of theory, was not accurate enough to precisely reproduce the experimental results. Here, we show that the missing piece of the puzzle is the nonrigidity of the H 2 molecule, neglected in the theoretical study of ref 9. It has been recently shown that taking into account the monomer nonrigidity effects is necessary to obtain precise agreement with the spectroscopic or scattering experiments. 13−17 In the present Article, we demonstrate the striking importance of the monomer flexibility in low-energy molecular anisotropic collisions. Only by inclusion of vibrations of the molecule in description of the complex are we able to correctly predict results of subtle scattering experiments with no fine-tuning to the experimental data whatsoever. Moreover, we present how to incorporate, in a simple and effective way, the nonrigidity effects into theoretical studies for collisions of excited species. Our approach involves the calculation of derivatives of the interaction energy with respect to a varied molecular geometry. This is nontrivial in the case of such system, and we present how symmetry-adapted perturbation theory (SAPT) can be used to obtain stable numerical derivatives.

THEORY WITH RESULTS
Theoretical investigations of rotationally and vibrationally inelastic scattering in atom−diatom and next diatom−diatom systems, based on a full-dimensional treatment of the problem, have been intensively developed since the late 1990s (for a comprehensive review, see ref 18). The research was stimulated mainly by astrochemical observations and cold chemistry. The full-dimensional approach requires reliable potential energy surfaces, an advanced theoretical description, and finally high-end computing resources. These are the reasons why the full-dimensional scattering calculations are still limited to the relatively small interaction partners. Such calculations were started by Balakrishnan and coworkers 19 to study the quenching of H 2 vibrations in collisions with He. Since then, more ab initio calculations of cross sections and rate coefficients for vibrational relaxation in cold inelastic collisions of atoms with diatomic molecules were published. 20−22 Using the full-dimensional approach to describe processes that involve a change of vibrational state of an interacting molecule, like relaxation of vibrationally exited molecules in collisions with atoms, seems to be natural, yet it can be very expensive. However, one can think to avoid a full-dimensional description when only rotational excitations or deexcitations take place. In such cases, the rigid-rotor approximation has been widely employed. The direct comparison of the two theoretical models, one treating the molecule as a rigid rotor and the other with vibrations included for colliding O 2 with He in the subkelvin regime, was presented by Volpi and Bohn. 23 All these with other theoretical studies on reactive scattering 24−26 initiated by Balakrishnan and Dalgarno 27 made a seminal contribution to the emerging field of cold controlled chemistry, where unexpected quantum effects occur. Until recently, the full dimensional calculations for the van der Waals complexes containing two diatomic molecules have been limited to the simplest case, i.e., H 2 + H 2 . 28,29 Lately, the larger complexes involving other atoms have been considered. The systems that are of interest to the astrochemical community, such as H 2 + CO, 16 35 and H 2 + SO, 36  In our case of scattering H 2 and He*, the helium atom is in the excited state. On the other hand, there is no vibrational excitation of the hydrogen molecule, which allows us to employ the rigid-rotor formalism in the scattering calculations. In principle, theoretical consideration of our problem can be divided into two steps: First is the preparation of the most reliable interaction energy surface, and second is the scattering calculations. Since the complex comprises the metastable atom He* and the diatomic molecule H 2 , the positions of the nuclei can be described by three 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, one has to use the three variables (R, θ, r) to parametrize the interaction energy surface and next to perform the scattering calculations. However, even for relatively small atom−diatom systems, the full-dimensional treatment is still rather rare. In most applications, the quantum scattering calculations are performed within the rigid-rotor approximation, i.e., assuming that the molecules in the complex are rigid. This widely used approximation is well physically motivated, since the internal vibrations of interacting molecules are much more energetic than the intermolecular modes. 13 Usually, in the rigid-rotor calculations, the rigid monomer interaction energy surfaces are used, obtained from ab initio calculations for the monomers with fixed geometries. Such calculations have been employed also in ref 9 to study the PI reaction of He* + H 2 . However, it has been shown very recently 16 that if in the rigid-rotor scattering calculations one uses the interaction energy averaged over the vibrations of the monomers then the results are closer to the full-dimensional calculations and experimental data. It has been demonstrated that the vibrationally averaged surfaces perform better than the rigid-monomer ones also in predictions of other physical properties, like rovibrational spectra [13][14][15]38,39 or virial coefficients. 40,41 A main drawback of such an approximation is that, in principle, to obtain the vibrationally averaged surface one has to know the corresponding full-dimensional one. To avoid construction of a full-dimensional surface and minimize the number of geometries for which ab initio calculations have to be performed, we use the method developed in ref 42. The details of this method, referred to as the Taylor-expansion method, are presented in the Computational Details section. Thus, in the current study, we try to capture the nonrigidity Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article effect by using the interaction energy surface averaged over the vibrations of H 2 and rigid-rotor scattering calculations. We have to emphasize that although the averaged surface depends only on the intermolecular coordinates (R, θ) it cannot be regarded as the rigid one, since to obtain it one has to perform ab initio calculations for various values of the internuclear distance r. Such a surface includes in an effective way information about the nonrigidity effects of the complex.
To use the Taylor-expansion method, one has to calculate the derivatives of the interaction energy with respect to the internal coordinates of the molecule. We have employed the SAPT method 43 to obtain these derivatives, instead of the supermolecular approach based on the CCSD(T) method used in the previous applications. The complex including the metastable helium is much more demanding with respect to calculations of the interaction energy than the complexes composed of the electronic ground state molecules previously studied. The combination of the Taylor-expansion method and SAPT allowed us to calculate the interaction energy with a precision sufficient to calculate numerical derivatives.
The most common technique to obtain PI rate coefficients is to use the complex potential in which the imaginary part describes the losses due to the ionization process. 44,45 With such a potential, one can solve the Schrodinger equation for the nuclear coordinates, for instance, using the close-coupling scattering method. More recently, two of us developed a new approach based on adiabatic theory and scattering theory for cold collision experiments 46,47 dubbed as adiabatic variational theory (AVT). AVT together with the dedicated new closedform expression for PI rate coefficients 48 provides a relatively simple method to implement and was used in our investigation. This technique, where the diatom is treated as a rigid rotor, allows one to uncouple the rotations of the diatom and the complex from the atom−molecule separation.
In Figure 1, we present the reaction rate coefficients of He* with para-H 2 in the ground rotational state (j = 0) and ortho-H 2 in the first excited rotational state (j = 1). The strong effect of anisotropy on the resonant structure was discussed in details in refs 9 and 47; however, it is worth mentioning that in the interaction between He* and para-H 2 (j = 0) only the isotropic part of the potential is probed because the wave function of the molecular hydrogen in the lowest rotational state is spherically symmetric. When the interacting H 2 molecule is in the j = 1 rotational state, the leading term of anisotropy of the potential contributes directly into the effective interaction and can firmly affect the positions of resonances. As demonstrated in ref 47, by excluding from calculations the orientation-dependent part, the low-temperature resonance at collision energy of k B × 0.27 K is missed in the theoretical results. It shows that this peak, as opposed to the peak at 2.4 K, arises totally from the anisotropic interaction. It takes place when the excited helium atom collides with ortho-H 2 (j = 1). Hence, the rate structures are entirely different for different rotational states of the hydrogen molecule.
The evidence of monomer nonrigidity effects in cold anisotropic collisions is demonstrated in Figure 1. It is clearly seen that the results with the interaction potential from the CCSD(T) method supplemented by the FCI correction are still not satisfactory. We found out that the discrepancy between the reported theoretical and experimental data is not due to the incompleteness of the used basis set as the authors of ref 9 stated but due to the assumed stiffness of the molecule. By adding the correction corresponding to the flexibility of the diatom, termed as δE int flex , where the vibrations are averaged, we obtained an excellent agreement with the measurements over the whole range of temperatures. The nonrigidity correction shifts the energy of the resonance at 0.27 K to the position matching the experimental data. Also, the magnitude of the rate coefficient curve is about 25% larger than the one calculated without the nonrigidity correction. These results are very similar to those obtained in ref 9 with the artificial 0.4% increase in the correlation part of the interaction energy. Note that in our entire calculations we did not apply any scaling or fitting parameters as well as we did not shift the final results to adjust to the experiment. Our results are slightly below the experimental ones, but the latter have been normalized to the absolute scale according to thermal rate observations at 300 K (see the Methods section in refs 8 and 9), and this procedure introduced a systematic error to the experimental data much larger than the vertical discrepancy. In the Supporting Information, we provide a figure corresponding to Figure 1, with the rate coefficient curves shifted by a constant value to match the normalized experimental data at the collision energy around k B × 2.4 K, as was done in ref 9. After such additional Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article "normalization", the theoretical resonance structure around 0.27 K perfectly agrees with the experiment. In the experimental data for para-H 2 , given in the upper panel of Figure 1, one can see a flat bump around the collision energy of k B × 0.3 K not discussed in ref 9. If a part of the plot with the bump is magnified, as in the inset of Figure 2, one can see that the experimental points forming this feature have error bars smaller than a magnitude of the bump; thus, one can suspect that it is not an accidental effect. On the basis of our theoretical considerations, we can try to explain an origin of that feature. One of the beams used in the Narevicius' experiment is formed of the para-H 2 molecules. However, a purity of para-H 2 was limited to 98%. 9 Therefore, we have added 2% of the reaction rate coefficients of ortho-H 2 (j = 1) given in the lower panel of Figure 1 to 98% of the reaction rate coefficients of para-H 2 (j = 0) given in the upper panel of Figure 1; the resulting rates are presented in Figure 2. Now one can see, in the inset of this figure, that the bump appears on the theoretical curve, properly predicting the position and the shape of its experimental counterpart. The explanation presented above can be verified in an experiment similar to that of Klein et al. 9 with amount of ortho-H 2 gradually increasing in the beam to observe if the magnitude of the bump changes accordingly.

COMPUTATIONAL DETAILS
3.1. Vibrationally Averaged Surface. To construct an interaction energy surface averaged over the vibrations of the monomers, ⟨V⟩, one can take advantage of the fact that the molecules in the complex preserve their identity, the frequencies of internal vibrations are much higher than those of the intermolecular modes, and thus, averaging can be performed over the internuclear coordinates similarly to the adiabatic approximation in the electronic-structure theory. 13 Since we are interested in the He* + H 2 complex, let us limit our consideration to the atom−diatom case. If the Jacobi coordinates are used, the interaction energy surface V(R, θ, r) can be represented as the truncated Taylor expansion around some reference geometry r c The higher order terms can be neglected if only modest deformations of the monomer are allowed, as those corresponding to a few lowest vibrational states of the H 2 molecule. The f 0 function is in fact the rigid monomer twodimensional surface calculated for the H−H separation equal to r c , while the remaining terms account for the nonrigidity effects, i.e., the dependence of the surface on the internal coordinate of H 2 . The potential V TE from eq 1 can be easily averaged over the v vibrational state of the monomer, and the resulting formula reads We can use the ⟨V TE ⟩ surface as an approximation to the ⟨V⟩ one, i.e., assume ⟨V⟩ ≈ ⟨V TE ⟩. The values of ⟨r⟩ v and ⟨r 2 ⟩ v can be calculated from the theoretical properties of the monomer or even estimated from the empirical spectroscopic constants. The surfaces averaged according to the approximation of eq 2 and the corresponding formula for the diatom−diatom case were used in the rigid-rotor dynamical calculations, both the bound states and the scattering, and provided the results in excellent agreement with the experimental ones. [14][15][16]38,39 In practical applications, there is no need to know the surface ⟨V TE ⟩ for any values of the (R, θ) coordinates, but it would be enough to calculate it on the grid points, for instance, the ones used in the scattering calculations. For each grid point in (R, θ), we can compute the interaction energy f 0 and the numerical values of the f 1 and f 2 derivatives. Of course, one can also calculate ⟨V TE ⟩ for a given grid of geometries and then fit an analytical function to obtain the surface. From eq 2, one can see a useful feature, namely, that the values of f 1 and f 2 may be calculated on a different level of theory than f 0 , e.g., the leading term f 0 on the highest possible level, and the values of f i defining the higher order terms, specifying the dependence of V TE on r, may be calculated at a somehow lower level of theory.
Since the firstand second-order terms are much smaller than the leading term f 0 , such an additional approximation only slightly increases the uncertainty of V TE or ⟨V TE ⟩, whereas it may significantly reduce the computational effort required. Such a strategy has been applied, for instance, to the H 2 + CO complex and led to the accurate rovibrational spectra. 14,15 To prepare the vibrationally averaged surface for the He* + H 2 complex, we have used the formula of eq 2 in the following way. The leading term f 0 was set to be equal to the rigidmonomer interaction energy of ref 9. That interaction energy can be written as E int CCSD(T) + δE int FCI , using our notation, and was obtained as a sum of the interaction energy calculated at the CCSD(T) level and the FCI correction for the intramolecular distance 1.4487 bohr. Thus, we have to set r c to be equal to the same value to make our choice of f 0 consistent with eq 2. Since Figure 2. Reaction rate coefficients of He(1s2s, 3 S 1 ) with H 2 with respect to relative energy (in K) between the colliding subsystems. The theoretical rate coefficients have been calculated based on the interaction potential obtained using CCSD(T) with additional corrections resulting from FCI and inclusion of the effect of the hydrogen molecule nonrigidity. The black solid curve represents results for 100% of para-H 2 (j = 0), whereas the magenta dashed− dotted curve represents results for the mixture of 98% of para-H 2 (j = 0) and 2% of ortho-H 2 (j = 1). The inset exhibits the region with a small bump in the millikelvin regime. The experimental data are shown as black points with error bars. 9 Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article we consider the ground vibrational state of H 2 , v = 0, the values of ⟨r⟩ 0 and ⟨r 2 ⟩ 0 for j = 0 were set to the round-off values from ref 49 equal to 1.4487 and 2.1270 bohr, respectively. The equality of r c and ⟨r⟩ 0 is of course not accidental because the interaction energy of ref 9 was calculated for the rigid H 2 molecule with the vibrationally averaged geometry. With such a choice of r c , the first-order term in eq 2 vanishes, and the nonrigidity effect is fully described by the second-order term. The values of ⟨r⟩ 0 and ⟨r 2 ⟩ 0 for the molecule in the first excited rotational state (j = 1) are slightly different: 1.4509 and 2.1334 bohr, respectively. Therefore, for this case, the firstorder term in eq 2 contributes to the nonrigidity correction.
In the standard, supermolecular approach, the interaction potential is difficult to obtain due to the fact that our potential is not the ground-state one and is coupled to the scattering state of the He + H 2 + + e − system. In ref 50, it was shown that using a carefully tailored start guess it is possible to converge the CCSD(T) interaction potential and that also SAPT 43,51,52 provides a good representation of the short-range potential. Nonetheless, it is very difficult to stabilize first and second derivatives of the interaction potential with respect to the nuclear coordinate motion in the supermolecular method as it inherently relies on subtraction of big numbers, and loss of accuracy is unavoidable. We previously stated that the derivative of the interaction energy can be obtained at a lower level of theory, and it still catches the essential physics. As a matter of fact, the application of SAPT greatly facilitates the calculation of nonrigidity effects. Since in SAPT the interaction energy is obtained directly from the wave function of monomers, it is stable and inexpensive. This is due to the fact that in this method we obtain the interaction energy directly as a sum of the perturbation theory terms in which the expansion parameter is the interaction potential between monomers. Here, we use the interaction energy which is the sum of first two terms of perturbation series in the interaction operator between H 2 and helium analogously to ref 50. The values of f 1 and f 2 were calculated from the four-and five-point central-difference formula, respectively, with h = 0.05 bohr. For each grid point in (R, θ), the interaction energy was calculated for five separations, r c + kh, where k = −2, − 1, 0, 1, 2. The SAPT calculations were carried out with the d-aug-cc-pVTZ basis set.
Let us now discuss the uncertainty of the flexibility correction δE int flex . In Table 1, we gathered second derivatives of the components of the interaction energy for 10.5 bohr for a T-shape and linear geometry. As one can see, the dispersion energy by far dominates the total flexibility correction. Given how good overall performance of SAPT was for the classically allowed region, 51 one can safely assume that the effect of higher order SAPT corrections will be marginal. There are two main uncertainties related to the dispersion derivative used in this Article: basis set incompleteness and time-dependent Hartree−Fock (TDHF) approximation. 51,52 To address the first uncertainty, we performed test calculations using the daug-cc-pVQZ basis set for a few geometries around the minimum in T-shape and linear configurations. To estimate the uncertainty of the dispersion and the exchange−dispersion components due to the basis truncation, these components can be extrapolated using the standard basis set extrapolation technique, 53 while one cannot use such extrapolation for the other components present in our SAPT expansion. We found that for the linear geometry the second derivative is underestimated by about 2%, while for the T-shape by about 10%. Since in the scattering calculations we used the interaction energy surface expanded in Legendre polynomials (see Section 3.2), it is interesting how the uncertainties of the second derivative transform to the uncertainties of two leading radial isotropic (V 0 ) and anisotropic (V 2 ) terms in that expansion. It turned out that the uncertainties amount to about 6% (underestimation) and 3% (overestimation) for V 0 and V 2 , respectively. In absolute numbers, these values are well below the uncertainty of basis set incompleteness for the E int CCSD(T) part of the total interaction energy. To address the uncertainty of the TDHF method, let us note that this model only slightly overestimates the dispersion energy for the metastable helium dimer by about 3%. 52 Similarly, the comparison of the longrange isotropic C 60 coefficient (i.e., the leading term in the inverse power expansion of V 0 ) obtained with TDHF 50  is proportional to overall performance of the dispersion, one can conclude that the inaccuracy due to the TDHF method is marginal and contributes to less than 0.01 cm −1 in the minimum range.
In Table 2, we present, for selected geometries of the complex, the values of the leading part of the interaction energy E int CCSD(T) obtained at the CCSD(T) level of theory and the values of various corrections to this energy. One can see that for the geometry close to the global minimum of the interaction energy surface, for θ = 0°and R = 10.5 bohr, the value of δE int flex is equal to −0.29 cm −1 , and although it seems to be small in the absolute scale, it amounts to about 2% of the E int CCSD(T) + δE int FCI interaction energy. For the same distance and θ = 90°, δE int flex amounts to 1% of the total energy. Thus, the δE int flex correction slightly changes the anisotropy of the potential. It is also interesting that the ratio of δE int flex to δE int FCI is significantly different for θ = 0°and 90°and amounts to 0.60 and 0.15, respectively, that shows that the anisotropy of these two corrections is completely different. The ratio is almost constant for the whole range of values of R at the same values of θ. The most important comparison one can draw from Table 2 is between the δE int flex and 0.4%E int corr corrections. For θ = 90°, they are almost equal, whereas for θ = 0°the values of , and extrapolated values of E disp (2) and E exch-disp (2) . For the induction and dispersion energies, we used the TDHF model; 50−52 exchange and electrostatic interactions were obtained from the Hartree−Fock densities. The derivatives with respect to r were calculated at r = r c , in units cm −1 /(bohr) 2 .
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article δE int flex are about two times larger than the 0.4%E int corr ones. Nevertheless, there is no doubt that the 0.4%E int corr correction introduced in ref 9 to reproduce the experimental rate coefficients can be recognized, on the basis of our investigation, as the result of the nonrigidity effect of H 2 on the interaction energy.
In ref 9, the effect on the interaction energy surface beyond the Born−Oppenheimer approximation was neglected on the basis of the analysis of the value of the diagonal Born− Oppenheimer (DBO) correction 55 calculated at the minimum of the surface. Here, we performed an extended analysis for two angular orientations, θ = 0°and 90°, several values of R, and the H 2 geometry fixed at r c . The DBO correction to the interaction energy, δE int DBO , was obtained by subtracting from the DBO correction calculated for the complex the value calculated for the monomer at large separations of interacting species. The calculations were performed at the CCSD level of theory, with the aug-cc-pVTZ basis set augmented by the midbond functions. This δE int DBO correction causes a small positive shift of the interaction energy, smaller than about 0.06 cm −1 in the minimum region, i.e., five times smaller than the value of δE int flex at this geometry. The values of δE int DBO for some other geometries are given in Table 2. For the values of R smaller than 10.0 bohr for θ = 0°and 9.0 bohr for θ = 90°, we have problems to converge the calculations of δE int DBO at the CCSD level; thus, for small values of R in the scattering calculations, we have used the values of δE int DBO extrapolated from the region of R greater than or equal to 10.0 and 9.0 bohr for θ = 0°and 90°, respectively. From Table 2, one can see that for geometries close to the geometry of the minimum, θ = 0°a nd R = 10.5 bohr, the value of δE int DBO amounts to 20% of δE int flex . The δE int DBO /δE int flex ratio is similar also for the same value of R and θ = 90°. However, one can observe that if the value of R increases then the value of δE int DBO decreases to zero faster than δE int flex . For instance, already for R = 12 bohr, δE int DBO amounts to only 10% of δE int flex , whereas for 14 bohr this ratio drops below 4%. This feature of the ratio between the δE int DBO and δE int flex corrections means that one can expect that adding δE int DBO to E int CCSD(T) + δE int FCI + δE int flex should not change significantly the calculated reaction rate coefficients, since in the major part of the range of the propagation of the resonance wave function the δE int DBO correction is very small in comparison with other components of the interaction energy. Indeed, in Figure S1 of the Supporting Information, one can see that the curve representing the rate coefficient calculated with the E int CCSD(T) + δE int FCI + δE int flex + δE int DBO surface is very close to the curve representing the calculations with the E int CCSD(T) + δE int FCI + δE int flex one. Concluding, the δE int DBO correction has a tiny effect on the position and shape of the considered resonances.
Finally, let us discuss the effect of basis set incompleteness. It is difficult to estimate since the basis set extrapolation techniques are questionable to use in the present case, where the system is not in its ground state and the system concerned is essentially a resonance. The basis set used in the present calculations at the CCSD(T) level, aug-cc-pV6Z with bond functions, is already well saturated for the dispersion energy which dominates the interaction energy for the considered system. A computationally expensive increase in the basis set to aug-cc-pV7Z shifts the interaction energy in the global (linear) minimum by about −0.022 cm −1 , while for the local minimum at a T-shape geometry (for 10.85 bohr) the shift is about −0.027 cm −1 . Close to the inner turning points at 9.5 bohr, these values are −0.034 and −0.045 cm −1 , respectively. Unfortunately, one cannot perform an extrapolation to the complete basis set, since the resonant nature of the interaction affects the stability of ab initio calculations at the level of accuracy of the order of 0.01 cm −1 . Nonetheless, one should bear in mind that such an interaction energy should not be extrapolated in the usual sense but rather than that stabilized. For this particular system, the coupling between the , δE int FCI , and 0.4%E int corr were calculated for purposes of ref 9, while δE int flex and δE int DBO in this work. Two angular orientations, the linear one (θ = 0°) and the T-shape one (θ = 90°), and selected values of the intermolecular separation R are represented. Energies are given in cm −1 . b Extrapolated value.
Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article continuum and the bound state is weak; hence, the potential obtained by the standard quantum chemistry methods gives right answer for right reasons. An improvement of theory to go below the 0.01 cm −1 accuracy is a formidable task, and a new approach would be needed to address such demands. Since the real potential is dominated by the dispersion energy which is a variational quantity, the more complete basis set implies a deeper potential. However, if one proceeds from the aug-cc-pV6Z basis set to the aug-cc-pV7Z one, a shift of the interaction energy surface by 0.02−0.03 cm −1 at a minimum separation is an order of magnitude smaller than in the case of the δE int flex correction. Thus, the effect of the basis set incompleteness on the positions and intensities of the resonances is negligible in comparison to the effect caused by the δE int flex correction and is also much smaller than the effect of the δE int DBO correction. 3.2. Expanding the Surface. According to the works of refs 56 and 57, we may expand the vibrationally averaged surface ⟨V⟩ in Legendre polynomials For the collision of an atom with a homonuclear diatomic molecule, the index η is even due to symmetry reasons, ⟨V⟩(R, −θ) = ⟨V⟩(R, θ). In other words, terms for odd η vanish. Thus, the two leading terms are given by V 0 (R) and V 2 (R)(3 cos 2 θ − 1)/2, where V 0 (R) and V 2 (R) are radial isotropic and anisotropic interaction potential terms, respectively. It should be emphasized that the considered complex is not in the bound state but in the resonance one. The total potential energy surface is above the ionization threshold, thus the electronic state of the He(1s2s, 3 S 1 ) + H 2 system is embedded in the continuum of the He(1s 2 , 1 S 1 ) + H 2 + + e − system. Consequently, the total potential energy surface has to be complex where the imaginary part represents the ionization width (inverse lifetime). Two new approaches have been lately developed for PI widths: One is based on the Fano-algebraic diagramatic construction method, 58 and the next one uses the stabilization method with an analytical continuation via the Padéapproximant. 59−61 In our studies, we took the imaginary part of V 0 and V 2 from ref 60, obtained by the latter technique.
The calculated isotropic, V 0 , and anisotropic, V 2 , radial interaction potential terms obtained on three levels of theory are presented in the upper panel of Figure 3. One can see that adding the δE int FCI correction to E int CCSD(T) apparently changes both V 0 and V 2 , but in opposite directions: V 0 becomes deeper and V 2 slightly shallower. The subsequent addition of δE int flex to the E int CCSD(T) + δE int FCI interaction energy makes the resulting V 0 even deeper, whereas for the V 2 term the effect of δE int flex almost cancels the effect of δE int FCI . The lower panel of Figure 3 presents how the δE int FCI and δE int flex corrections to the interaction energy enter the V 0 and V 2 terms of the interaction potential.
3.3. Adiabatic Variational Theory. To solve the Schrodinger equation for the complex with the previously prepared interaction potential and consequently to calculate rate coefficients, we used the AVT approach that has been developed for cold atom−molecule collision experiments. 46,47 This technique has been recently successfully applied for the He(1s2p, 3 P 2 ) + H 2 system. 62 It enables one to reduce the complexity of the problem enhancing the computational performance without losing physical essence. Within AVT, we represent the potential (eq 3) in a basis set consisting of many angular functions. In our case, the single angular function is a product of two spherical harmonics: One is responsible for the description of rotations of the molecule, whereas the other one is for the description of rotations of the whole complex. Such a matrix constructed for a given intermolecular distance R needs to be diagonalized providing a set of eigenvalues. The process has to be repeated for different values of R. The obtained eigenvalues, after ordering, create so-called adiabats (effective potentials) depending on R. Then, from the practical reasons, all adiabats are shifted to get asymptotes at zero. Therefore, for each of them, the dissociation threshold is at zero. Next, we solve the one-dimensional Schrodinger equation many times, each time with a different adiabat treating R as a variable. At this point, any technique can be used, for example, by fitting an analytical function to the adiabat and spanning the wave function space in a basis of trial functions. 63 Finally, we apply the simple and easy to implement formula for reaction rate coefficients that has been derived based on AVT and non-Hermitian scattering theory. 48 Only the information about complex eigenenergies, the reduced mass of the atom−diatom system, and the rotational state of the molecule are required. In calculations, we used 21 partial waves corresponding to endover-end angular momenta of the complex (l = 0, 1, ..., 20). The Schrodinger equation was solved by the DVR with a box size of 500 bohr and with 2000 basis sine functions. The results are fully converged with respect to the number of partial waves and of basis functions. The calculated reaction rate coefficients Journal of Chemical Theory and Computation pubs.acs.org/JCTC Article are convoluted over the experimental collision energy spread (8 mK). Neither scaling nor fitting parameters have been used in our calculations.

CONCLUSIONS
Low-energy collision experiments allow one to directly and clearly reveal the true nature of interactions. Only at a high level of theory are we able to predict or confirm experimental data as well as understand quantum phenomena in chemical reactions. In this Article, we have demonstrated a significant role of monomer nonrigidity effects on the position and intensity of scattering resonances in anisotropic cold molecular collisions. We have investigated excited metastable helium atoms colliding with hydrogen molecules in the temperature range from a few dozen kelvins down to 1 millikelvin. We have provided the most accurate interaction energy surface that can be treated as the reference one for all semiclassical algorithms as well as for novel quantum chemistry methods for molecular systems in the resonance state, where the Hamiltonian is real or complex. The calculated rate coefficients are in remarkable agreement with the measurements. 9 We have demonstrated that the discrepancy between the experimental and theoretical results discussed in ref 9 is due to the assumption that the H 2 molecule is rigid. We have complemented the theoretical description of the interaction energy of the complex by inclusion of its dependence on the flexibility of the molecule. For such a challenging complex as He* + H 2 , it was necessary to use the SAPT method to calculate the interaction energy with a precision sufficient to obtain numerical derivatives required by the Taylor-expansion method to describe the flexibility effects. Our results exhibit that the approach beyond the commonly used rigid-rotor approximation is indispensable even when rigorous state-of-the-art computations are performed at the FCI level of theory. By thorough analysis of the uncertainties, we have shown that the two largest of them, caused by basis-set incompleteness and generated by the Born−Oppenheimer approximation, are a few times smaller than the correction to the interaction energy due to the nonrigidity of the monomer. Thus, this is the only missing important correction we had to take into account to achieve the full agreement between the theory and experiment for the He* + H 2 collision pair. Our finding concerning the significance of the nonrigidity effects is not limited to a specific complex considered here nor restricted to the PI reaction process. We believe that it can be vitally important in precisely controlled cold chemistry, where quantum effects in chemical reactions dominate.
■ ASSOCIATED CONTENT