Nonrigidity effects -- a missing puzzle piece in the description of low-energy anisotropic molecular collisions

TL;DR

This study demonstrates that accounting for the nonrigidity of H₂ molecules is essential for accurately predicting resonance positions and reaction rates in low-energy anisotropic molecular collisions, significantly improving theoretical-experimental agreement.

Contribution

The paper introduces the inclusion of molecular nonrigidity into interaction potential models, revealing its critical role in accurately describing cold molecular collision resonances.

Findings

Nonrigidity effects are crucial for correct resonance predictions.

Inclusion of molecular flexibility improves agreement with experimental reaction rates.

Neglecting nonrigidity leads to significant discrepancies in cold collision modeling.

Abstract

Cold collisions serve as a very sensitive probe of the interaction potential. In the recent study of Klein et al. (Nature Phys. 13, 35-38 (2017)) 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 H$_2$ molecule on the resonant structure, absent in the previous study, is critical to predict 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.

Equations

$$V_{\mathrm{TE}}(R,\theta ,r)=f_0(R,\theta )+f_1(R,\theta )(r-r_c)+\frac{1}{2}{f}_2(R,\theta )(r-r_c{)}^2,$$

$$⟨V_{\mathrm{TE}}⟩(R,\theta )$$

$$⟨V⟩(R,\theta )=\sum _{\eta }{V}_{\eta }(R)P_{\eta }(\cos \theta ).$$

Full text

The breakthrough in controlling the movement and internal degrees of molecules with external fields van de Meerakker et al. (2008), which started about 20 years ago, currently allows to study the collisions, reactivity and properties of molecules in very 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 Côté (2014); Bohn et al. (2017). Cold collision experiments realized by merging two supersonic beams have become an important technique for studying chemical reactions in temperatures near 1 K Henson et al. (2012); Lavert-Ofir et al. (2014); Jankunas et al. (2014, 2015); Shagam et al. (2015); Klein et al. (2017) 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 Henson et al. (2012); Lavert-Ofir et al. (2014); Shagam et al. (2015); Klein et al. (2017). In such a reaction, an electron is moved from the molecule to the only partially-occupied orbital of 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 Dalgarno (2000); Shagam et al. (2015). The first observation of metastable helium atom in the atmosphere of one of exoplanets Spake et al. (2018) 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 co-workers directly probed the anisotropy in atom–molecule interactions through orbiting resonances by changing the rotational state of the molecule Klein et al. (2017). That work reveals a crucial role of the anisotropy of the interaction energy, due to various orientations of the H2 molecule in the complex of He(1s2s,3S1) ($\equiv$ He*∗) with H2, in the dynamics in sub-kelvin regime. To elucidate physical phenomena presented in their novel experiment, the authors of Ref. Klein et al. (2017) used the state-of-the-art first-principles calculations. The interaction energy of the complex comprised 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_{\rm int}^{\rm 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 Vogels et al. (2018). However, in the case of He∗*+H2, the values of reaction rate coefficients calculated with the potential based on $E_{\rm int}^{\rm CCSD(T)}$ did not agree with experiment even qualitatively. Thus, the correction to the interaction energy, denoted here by $\delta E_{\rm int}^{\rm FCI}$, was calculated from the full configuration interaction (FCI) method and added to the $E_{\rm int}^{\rm CCSD(T)}$ energy. The agreement of the rate coefficients calculated from the $E_{\rm int}^{\rm CCSD(T)}$+$\delta E_{\rm int}^{\rm FCI}$ energies with the experimental values was improved but still qualitatively incorrect, since one additional resonance, not present in the experiment, was predicted for very 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, that can be viewed as adding the correction $0.4\%E_{\rm int}^{\rm corr}$ to the $E_{\rm int}^{\rm CCSD(T)}$+$\delta E_{\rm int}^{\rm 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 H2 molecule, neglected in the theoretical study of Ref. Klein et al. (2017). It has been recently shown that taking into account the monomer nonrigidity effects is necessary to obtain very precise agreement with the spectroscopic or scattering experiments Jeziorska et al. (2000); Jankowski et al. (2012, 2013); Faure et al. (2016); Stoecklin et al. (2017). In the present paper, we demonstrate for the first time the importance of the monomer flexibility in low-energy molecular anisotropic collisions. Only by inclusion of vibrations of the molecule in description of the complex, we are able to correctly predict results of very 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 any colliding molecules.

Results and discussion

Theoretical consideration of our problem can be divided into two steps: first, the preparation of the most reliable interaction energy surface and second, scattering calculations. Since the complex comprises the metastable atom He*∗* and the diatomic molecule H2, the positions of the nuclei can be described by three coordinates: the distance $R$ between He*∗* and center of mass (COM) of H2, the angle $\theta$ between the H2 bond and the COM–He*∗* direction, and the distance $r$ between the hydrogen nuclei. Thus, in principle one has to use the three variables $(R,\theta,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 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 Jeziorska et al. (2000). 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. Klein et al. (2017) to study the PI reaction of He∗+H2. However, it has been shown very recently Faure et al. (2016) that if in the rigid-rotor scattering calculations one uses the interaction energy averaged over the vibrations of the monomers, 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 Jeziorska et al. (2000); Jankowski (2008); Jankowski and Szalewicz (2008); Jankowski et al. (2012, 2013) or virial coefficients Garberoglio et al. (2014, 2017). A main drawback of such 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. Jankowski (2004). The details of this method are presented in the Methods section. Thus, in the current study we try to capture the nonrigidity effect by using the interaction energy surface averaged over vibrations of H2 and rigid-rotor scattering calculations. We have to emphasize that although the averaged surface depends only on the intermolecular coordinates $(R,\theta)$, it cannot be regarded as the rigid one, since to obtain it one has to perform ab initio calculations for various values of the intramolecular distance $r$. Such a surface includes in an effective way information about the nonrigidity effects of the complex.

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 Miller et al. (1972); Siska (1993). With such potential one can solve the Schrödinger 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 Pawlak et al. (2015, 2017) dubbed as adiabatic variational theory (AVT). AVT together with the dedicated new closed-form expression for PI rate coefficients Pawlak et al. (2018) 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 to uncouple the rotations of the diatom and the complex from the atom–molecule separation.

In Fig. 1 we present the reaction rate coefficients of He∗ with para-H2 in the ground rotational state ($j=0$) and ortho-H2 in the first excited rotational state ($j=1$). The strong effect of anisotropy on the resonant structure was discussed in details in Refs. Klein et al. (2017); Pawlak et al. (2017), however it is worth mentioning that in the interaction between He∗ and para-H2 ($j=0$) only the isotropic part of the potential is probed, because the wavefunction of the molecular hydrogen in the lowest rotational state is spherically symmetric. When the interacting H2 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. Pawlak et al. (2017), by excluding from calculations the orientation-dependent part, the low-temperature resonance at collision energy of $k_{B}\times 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-H2 ($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 Fig. 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. Klein et al. (2017) stated, but due to the assumed stiffness of the molecule. By adding the correction corresponding to the flexibility of the diatom, termed as $\delta E^{\rm flex}_{\rm int}$, 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. Klein et al. (2017) with the artificial $0.4\%$ increase of 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 in Refs. Shagam et al. (2015); Klein et al. (2017), and this procedure introduced a systematic error to the experimental data much larger than the vertical discrepancy. In the Supplementary Information we provide a figure corresponding to Fig. 1, with the rate coefficient curves shifted by a constant value to match the normalized experimental data at the collision energy around $k_{B}\times 2.4$ K, as it was done in Ref. Klein et al. (2017). After such additional “normalization”, the theoretical resonance structure around 0.27 K perfectly agrees with the experiment.

In the experimental data for para-H2, given in the upper panel of Fig. 1, one can see a very flat bump around the collision energy of $k_{B}\times 0.3$ K, not discussed in Ref. Klein et al. (2017). On the basis of our theoretical considerations we can explain an origin of that feature. One of the beams used in the Narevicius’ experiment is formed of the para-H2 molecules. However, a purity of para-H2 was limited to $98\%$ Klein et al. (2017). Therefore, we have added $2\%$ of the reaction rate coefficients of ortho-H2 ($j=1$) given in the lower panel of Fig. 1 to $98\%$ of the reaction rate coefficients of para-H2 ($j=0$) given in the upper panel of Fig. 1; the resulting rates are presented in Fig. 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.

Conclusions

Low-energy collision experiments allow to directly and clearly reveal the true nature of interaction. Only at a high level of theory we are able to predict or confirm experimental data as well as understand quantum phenomena in chemical reactions. In this paper, 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 to 1 millikelvin. We have provided the most accurate interaction energy surface that can be treated as the reference one for novel quantum chemistry methods for molecular systems in the resonance state. The calculated rate coefficients are in remarkable agreement with measurements Klein et al. (2017). We have demonstrated that the discrepancy between the experimental and theoretical results discussed in Ref. Klein et al. (2017) is due to the assumption that the H2 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. It shows 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. The correction to the interaction energy due to the nonrigidity of the monomer is a few times larger than the uncertainties caused by basis-set incompleteness or generated by the Born–Oppenheimer approximation. Our finding 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.

Methods

Vibrationally averaged surface

To construct an interaction energy surface averaged over the vibrations of the monomers, $\langle V\rangle$, one can take advantage of the fact that the molecules in the complex preserve their identity and frequencies of internal vibrations are much higher than those of the intermolecular modes, and perform averaging over the intramonomer coordinates similar to the adiabatic approximation in the electronic-structure theory Jeziorska et al. (2000). Since we are interested in the He*∗*+H2 complex, let us limit our consideration to the atom-diatom case. If the Jacobi coordinates are used, the interaction energy surface $V(R,\theta,r)$ can be represented as the truncated Taylor expansion around some reference geometry $r_{c}$ Jankowski (2004)

[TABLE]

where $f_{0}(R,\theta)=V(R,\theta,r_{c})$, f_{1}(R,\theta)=\frac{\partial V(R,\theta,r)}{\partial r}\big{|}_{r=r_{c}}, and f_{2}(R,\theta)=\frac{\partial^{2}V(R,\theta,r)}{\partial r^{2}}\big{|}_{r=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 H2 molecule. The $f_{0}$ function is in fact the rigid monomer two-dimensional surface calculated for the H–H separation equal $r_{c}$, while the remaining terms account for the nonrigidity effects, i.e., the dependence of the surface on the intramolecular coordinate. The potential $V_{\rm TE}$ from Eq. (1) can be easily averaged over the $v$ vibrational state of the monomer and the resulting formula reads

[TABLE]

We can use the $\langle V_{\rm TE}\rangle$ surface as an approximation to the $\langle V\rangle$ one, i.e., assume $\langle V\rangle\approx\langle V_{\rm TE}\rangle$. The values of $\langle r\rangle_{v}$ and $\langle r^{2}\rangle_{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 Jankowski (2008); Jankowski and Szalewicz (2008); Jankowski et al. (2012, 2013); Faure et al. (2016). In practical applications, there is no need to know the surface $\langle V_{\rm TE}\rangle$ for any values of the $(R,\theta)$ 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,\theta)$, 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 $\langle V_{\rm TE}\rangle$ for a given grid of geometries and then fit an analytical function to obtain the surface. From the expression (2), one can see a very 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_{\rm TE}$ on $r$, may be calculated at a somehow lower level of theory. Since the first and second order terms are much smaller than the leading term $f_{0}$, such an additional approximation only slightly increases the uncertainty of $V_{\rm TE}$ or $\langle V_{\rm TE}\rangle$, whereas it may significantly reduce the computational effort required. Such a strategy has been applied, for instance, to the H2–CO complex in Refs. Jankowski et al. (2012, 2013) and led to the very accurate rovibrational spectra.

To prepare the vibrationally averaged surface for the He*∗*+H2 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 rigid-monomer interaction energy of Ref. Klein et al. (2017). That interaction energy can be written as $E_{\rm int}^{\rm CCSD(T)}$+$\delta E_{\rm int}^{\rm 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 we consider the ground vibrational state of H2, $v=0$, the values of $\langle r\rangle_{0}$ and $\langle r^{2}\rangle_{0}$ for $j=0$ were set to the round-off values from Ref. Bubin and Adamowicz (2003) equal to 1.4487 bohr and 2.1270 bohr, respectively. The equality of $r_{c}$ and $\langle r\rangle_{0}$ is of course not accidental, because the interaction energy of Ref. Klein et al. (2017) was calculated for the rigid H2 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 $\langle r\rangle_{0}$ and $\langle r^{2}\rangle_{0}$ for the molecule in the first excited rotational state ($j=1$) are slightly different: 1.4509 bohr and 2.1334 bohr, respectively. Therefore, for this case the first-order term in Eq. (2) contributes to the nonrigidity correction.

In the standard, supermolecular approach the interaction potential is very 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. Hapka et al. (2013) it was shown that using a carefully tailored start guess it is possible to converge the CCSD(T) interaction potential and that also the symmetry-adapted perturbation theory (SAPT) Jeziorski et al. (1994); Hapka et al. (2012); Żuchowski et al. (2003) 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 wavefunction of monomers, the interaction energy is stable and very inexpensive. This is due to the fact that in SAPT 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 H2 and helium analogously to Ref. Hapka et al. (2013). 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,\theta)$, the interaction energy was calculated for five separations $r_{c}+kh$, $k=-2,-1,0,1,2$. The SAPT calculations were carried out with the d-aug-cc-pVTZ basis set. We estimated that uncertainty of calculation of flexibility correction in the interaction well is underestimated by about 2% for linear geometry and about 10% for T-shape geometry. In absolute numbers these uncertainties are below 0.01 cm*-1*. For discussion of uncertainty of $\delta E^{\rm flex}_{\rm int}$ see the Supplementary Information.

In Ref. Klein et al. (2017), 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 Handy et al. (1986) calculated at the minimum of the surface. Here we performed an extended analysis for two angular orientations, $\theta=0^{\circ}$ and $90^{\circ}$, several values of $R$, and the H2 geometry fixed at $r_{c}$. The DBO correction to the interaction energy, $\delta E_{\rm int}^{\rm 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 $\delta E_{\rm int}^{\rm 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 $\delta E^{\rm flex}_{\rm int}$ at this geometry. The values of $\delta E_{\rm int}^{\rm DBO}$ for some other geometries are given in the Supplementary Information. We have also added $\delta E^{\rm DBO}_{\rm int}$ to the $E^{\rm CCSD(T)}_{\rm int}+\delta E^{\rm FCI}_{\rm int}+\delta E^{\rm flex}_{\rm int}$ energy and have performed the scattering calculations. The resulting positions of resonances are almost unaffected by the addition of the $\delta E^{\rm DBO}_{\rm int}$ correction, as can be seen in the figure presented in the Supplementary Information.

Finally, let us discuss the effect of basis set incompleteness. It is very 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 of 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 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 the deeper potential. However, if one proceeds from the aug-cc-pV6Z basis set to the aug-cc-pV7Z one, the shift of the interaction energy surface by 0.02–0.03 cm*-1* at the minimum separation is an order of magnitude smaller than in the case of the $\delta E^{\rm flex}_{\rm int}$ correction. Thus, the effect of the basis set incompleteness on the positions and intensities of the resonances is negligible in comparison to the effects caused by the $\delta E^{\rm flex}_{\rm int}$ correction and is also much smaller than the effect of the $\delta E_{\rm int}^{\rm DBO}$ correction.

According to the works of Refs. Holmgren et al. (1977); Dubernet and Hutson (1994), we may expand the vibrationally averaged surface $\langle V\rangle$ in Legendre polynomials

[TABLE]

For the collision of an atom with a homonuclear diatomic molecule the index $\eta$ is even due to symmetry reasons, $\langle V\rangle(R,-\theta)=\langle V\rangle(R,\theta)$. In other words, terms for odd $\eta$ vanish. Thus, the two leading terms are given by $V_{0}(R)$ and $V_{2}(R)(3\cos^{2}\theta-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,3S1)+H2 system is embedded in the continuum of the He(1s2,1S1)+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 Fano-algebraic diagramatic construction method Yun et al. (2018), the next one uses the stabilization method with an analytical continuation via the Padé approximant Landau et al. (2016); Bhattacharya et al. (2017). In our studies, we took the imaginary part of $V_{0}$ and $V_{2}$ from Ref. Bhattacharya et al. (2017), 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 Fig. 3. One can see that adding the $\delta E_{\rm int}^{\rm FCI}$ correction to $E_{\rm int}^{\rm 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 $\delta E_{\rm int}^{\rm flex}$ to the $E_{\rm int}^{\rm CCSD(T)}$+$\delta E_{\rm int}^{\rm FCI}$ interaction energy makes the resulting $V_{0}$ even deeper, whereas for the $V_{2}$ term the effect of $\delta E_{\rm int}^{\rm flex}$ almost cancels the effect of $\delta E_{\rm int}^{\rm FCI}$. The lower panel of Fig. 3 presents how the $\delta E_{\rm int}^{\rm FCI}$ and $\delta E_{\rm int}^{\rm flex}$ corrections to the interaction energy enter the $V_{0}$ and $V_{2}$ terms of the interaction potential.

Adiabatic variational theory

To solve the Schrödinger 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 Pawlak et al. (2015, 2017). This technique has been recently successfully applied for the He(1s2p,3P2)+H2 system Bhattacharya et al. (2019). It enables one to reduce the complexity of the problem enhancing the computational performance without losing physical essence. Within AVT, we represent the potential (3) in a basis set consisting of many angular functions. The single angular function in our case is a product of two spherical harmonics — one is responsible for the description of rotations of the molecule, whereas the other one 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 Schrödinger equation many times, each time with a different adiabat treating $R$ as a variable. Here, any technique can be used, e.g., by spanning the wave function space in a basis of trial functions Pawlak and Bylicki (2011); Pawlak et al. (2014). 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 Pawlak et al. (2018). 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 end-over-end angular momenta of the complex ($l=0,1,...,20)$. The Schrödinger equation was solved by the DVR with the 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 are convoluted over the experimental collision energy spread (8 mK). Neither scaling nor fitting parameters have been used in our calculations.

Acknowledgements

We would like to thank Debarati Bhattacharya for providing the imaginary part for radial interaction potentials.

This research was financially supported by the National Science Centre, Poland (Grant No. 2016/23/D/ST4/00341 [MP], Grant No. 2015/19/B/ST4/02707 [PSŻ], and Grant No. 2017/25/B/ST4/01300 [PJ]).

Competing financial interests

The authors declare no competing financial interests.