Dissociative charge transfer through a conical intersection: quantum thermal rate constants up to 1000 K for the He⁺ + H₂ → He + H + H⁺ reaction

Abstract

In this work, a previous quantum wave-packet time-dependent study (D. De Fazio, A. Aguado and C. Petrongolo, Non-adiabatic quantum dynamics of the dissociative charge transfer He⁺ + H₂ → He + H + H⁺. Front. Chem., 2019, 7, A249) on the charge-transfer dissociation of the hydrogen molecule by the helium cation is continued, extending the calculation to eight rotational H₂ reactant states. New data are required to obtain convergent (within about 1 percent) Boltzmann-averaged thermal rate constants up to 1000 Kelvin, which are necessary to provide reliable reaction yields for astrophysical and cosmological computational models. To the best of our knowledge, these are the first quantum mechanical thermal rate constant calculations for a dissociative charge transfer process. As shown in the previous work for the effect of vibrations, a relevant role of rotations is found, and the use of roto-vibrational ground state rate constant only, as often employed in astrophysical model is a poor approximation, especially above room temperature. The large computational scaling of wave-packet methods for excited rotational reactant states, due to the linear increase of the number of projections of the diatomic rotational angular momentum along the internuclear axis, is handled efficiently in the calculations, providing an accurate documentation for developing approximations that are needed to extend the calculations to higher temperatures. A particular effort has been made to better clarify the reaction mechanism, which indicated the key roles of the geometrical phase and non-adiabatic effects. Additionally, the rich resonance pattern, which is the principal reaction mechanism at these temperatures and represents the main source of the computational load, is analyzed and rationalized in detail according to the resonance analysis.

---

1. Introduction

The reactivity of the helium atom has been a topic often neglected because its inertness as a noble gas restricts its involvement to only a few interesting chemical processes. Nowadays, the relevance of helium chemistry in astrochemistry and cosmology¹ has radically changed this view since the reactions of helium and those of its cations have arguably played crucial roles in the evolution of our universe, soon after the primordial Big-Bang nucleosynthesis.²

Known as a rare gas, helium is actually one of the most abundant elements, second only to hydrogen. Recently, renewed interest in the scientific community in its reactions came from the possibility of obtaining important data from orbiting³ or flying infrared telescopes⁴ that can observe starlight emissions outside the Earth's atmosphere. This permitted the observation of unequivocal fingerprints of the long searched⁵ HeH⁺ ion in the planetary nebula NGC 7027,⁶ confirming the hypothesis that the helium hydride cation was the first molecule formed in the Early Universe, and participates in a chemical network including the reactions He⁺ + H → HeH⁺ + hν and HeH⁺ + H → He + H₂⁺ whose theoretical quantum mechanical (QM) cross sections⁷ and rate constants⁸ were used for the data analysis.⁹ This synergy between observations and calculations points out the necessity of accurate QM theories of the reaction dynamics and of extensive calculations to account for astro-chemical findings.

The reaction dynamics of the HeH₂⁺ system have been employed as benchmarks in several experimental and computational studies (for a recent comprehensive review on the subject, see ref. 10). From a theoretical point of view, the presence of few electrons and low mass of this system enables highly accurate QM calculations; accordingly, several potential energy surfaces (PESs) for the electronic ground and excited states have been reported at different levels of theory.^11–26 Additionally, reactive scattering calculations of the direct and inverse reactions on the ground PES and for some isotopic variants have been obtained with a variety of reactive theories using classical (see e.g. ref. 27–29) or quantum mechanical time-independent^30–33 and wavepacket^34–38 methods. This theoretical information is also important because of the difficulty in experimentally producing well-collimated molecular beams of highly reactive He cations to perform high-resolution scattering measurements.^39–44

However, comprehensive dynamical studies of the reactions on the excited PESs are rare and mostly focused on the properties of the stable bound or quasi-bound triatomic complex in given electronic states of the system.^45–53 The first large-scale study of X̃²A′, òA′, and B̃²A′ adiabatic PESs and their non-adiabatic (NA) couplings was published in the pioneering work of ref. 22, fitted later by the analytical functionalities in ref. 24. More recently, an impressive work reporting the first four adiabatic PESs and their NA couplings⁵⁴ was presented in ref. 26, and the NA quantum dynamics on the lowest two ²A′ electronic states was calculated.⁵³ Large NA effects were also found in the dynamics of the ground-state proton transfer reaction He + H₂⁺ despite the large energy gap between the first and second adiabatic PESs for most of the physically relevant triatomic geometries. A closer agreement with the experimental results for the ground reactant roto-vibrational state (ref. 39 and 40) was found, supporting the accuracy of the presented results. However, no analytically fitted PES functions were supplied, and no reliable interatomic interaction was given in the long-range region, which are highly relevant for ionic systems (see e.g. ref. 20). Therefore, the convergence of the dynamics (see Fig. SM8 of the SI of ref. 26) needs to be confirmed to fully assess their reliability.

Another important well-known obstacle to the comparison with experiments and/or to use QM data in astrophysical models, is that the reactions of this system are highly sensitive to the reactants’ initial states. Large vibrational effects in the direct reaction of the ground electronic state were already evident in the first experimental investigations,^39–42 and the relevant rotational effects were confirmed in ref. 32. A rotational state resolution is therefore required for comparison with theoretical data, which is very demanding in molecular beam experiments of ionic species.⁴⁴

In ref. 52 we have tackled the first rigorous QM calculation of the dissociative charge transfer (DCT) process including the conical intersection (CI) between the first and the second electronic adiabatic à and B̃ states of the HeH₂⁺ system (see a schematic energy diagram in Fig. 1): the reaction He⁺ + H₂ → He + H + H⁺, all chemical species being in their ground electronic state. This reaction can occur in He⁺-rich interstellar molecular clouds, such as NGC7027⁶ and has recently been found to be of relevance to interpret astronomic observations outside of the solar system by modeling escape from the atmosphere of warm Neptune-sized exoplanets⁵⁵ on account of the effect of vibrations.⁵⁶

Fig. 1 Correlation diagram (in eV) of the adiabatic òA′ (red) and B̃²A′ (blue) and diabatic (1)²A₁ and (2)²A₁ PESs. All reactants and products are in their ground electronic states, except the unbound H₂⁺(A²Σ⁺_u). The energies of the minima of the van der Waals well and conical intersection seam are also indicated.

In particular, we compared the reaction dynamics from the lowest H₂ rotational level j₀ = 0 in the ground and first excited vibrational states ν₀ = 0 and 1 using the analytical PESs of ref. 24 (A′ PESs).

A main result from this study is that DCT can also occur at low collision energies notwithstanding the high energy location of the CI minimum (about 1.34 eV over that of reactants). This is possible because of the extremely long-lived resonances trapped in the van der Waals (vdW) well of the surface. They are particularly demanding computationally but still affordable with the present methodology. The geometry of the well (see details in Section 3) is very similar to the minimum CI. This intense resonance pattern survives in the kinetic observables, giving a broad maximum in the rate constants of around 100 K. A strong vibrational effect was also characterized, so the reactivity increased by more than an order of magnitude for the state of the first excited vibrational reactants.

In Fig. 2, we show the thermal populations (following the Boltzmann distribution) of the normal hydrogen molecules up to 1000 K.

Fig. 2 Thermal populations of a normal (75% ortho and 25% para) hydrogen molecule as a function of temperature.

As illustrated in Fig. 2, where several excited rotational states are significantly populated in the range of temperature (T) considered, restriction to the reactant ground states to the reaction kinetics appears as a poor approximation already from one hundred Kelvin.

To meet these requirements, in this work, we extend the calculations of ref. 52 by reporting initial-state-resolved integral cross sections and rate constants for H₂ in the ground vibrational state and rotational levels of j₀ ≤ 7, thus obtaining highly converged thermal rate constants for para-, ortho-, and normal-H₂ up to 1000 K. To the best of our knowledge, this is the first time that convergent quantum thermal rate constants have been presented for a CI DCT process.

The rest of this article is organized as follows. In Section 2, a brief account of the reactive computational method adopted (essentially the same as that used previously) is given together with a discussion of the convergence input parameters employed. In Section 3, we show the results obtained, including those from two useful approximations, and offer a convergence analysis for extending the calculations to higher temperatures. In Section 4, we analyze the obtained results to extract physical insights into the main reaction mechanism. In particular, a resonance analysis is performed to understand the complex resonance spectrum presented in the reactive observables and its evolution as the reactant rotational quantum number varies. Finally, in Section 5, we present the conclusions and perspectives.

2. Theory and calculations

The PESs and CI coupling are discussed in detail in previous works,^24,52 so here we present only a summary following the schematic correlation diagram of the electronic states plotted in Fig. 1, showing the symmetry, CI, and energetics of the (Ã, B̃) adiabatic and (1, 2) diabatic states. Note that the DCT channel is strongly exoergic by 6.25 eV but the CI barrier inhibits the reactivity. A small vdW well (with a depth of about 80 meV) near the CI configuration is also indicated. As we reported in detail later, although energetically small, this feature deeply affects the reactivity of the process. The (1, 2) diabatic PESs and the associated electronic coupling calculated in ref. 22 were used.

The reaction dynamics are performed by applying a well-consolidated quantum, time-dependent, real wavepacket (WP) method^57,58 implemented and parallelized to treat efficiently the main NA couplings (Renner–Teller, CI and spin–orbit effects) in different triatomic processes and adopted to study a large variety of reactions.^59,60 It is based on a general theory⁶¹ of NA effects in triatomic species; many details of the theory and an account of the numerical code can be found in the literature cited above. The parameters used in this specific calculation are also provided, and the physical quantities are defined.

We use atomic units (a.u.) unless otherwise specified. We describe the He⁺ + H₂ configurations by reactant Jacobi coordinates R, r, and γ in a body-fixed reference frame with the z axis along R and by a HeH₂⁺ spinless rovibronic Hamiltonian Ĥ, which contains the electronic Ĥ^el and the total angular momentum Ĵ operators. Ĥ is represented in an orthonormal basis of diabatic electronic states (1)²A₁ and (2)²B₂, radial grid |Rr〉, associated Legendre |jK〉, and symmetry Wigner states |K + p〉. Here (1)²A₁ and (2)²B₂ are coupled by Ĥ^el owing to the CI, and K is the Ĵ_z eigenvalue. Omitted are J and its space-fixed Z component in |K + p〉, where the total parity inversion parity is p = (−)^J+K_min with K_min = 0 or 1. The 2J + 1 values of K are thus factorized into two noninteracting groups, with K_min ≤ K ≤ J, of dimensions J + 1 or J according to K_min = 0 or 1, respectively.

The wavepacket Chebyshev propagation is performed using the diabatic PESs directly and therefore automatically contains (see e.g. ref. 62) the geometric phase (GP) effects,^63–65 which, as discussed at the end of Section 3, play a key role in the dynamics of this process.

Initial-state-resolved reaction probabilities are computed through the quantum, real WP formalism of Gray and Balint-Kurti,⁵⁷ which is essentially equivalent to the Chebyshev approach of Guo.⁶⁶ Briefly, an arccos mapping of the HeH₂⁺ time-dependent Schrödinger equation is solved recursively using a scaled and shifted Hamiltonian Ĥ_s and starting from an initial complex WP⁵⁷ |ψ₀〉 = |a₀〉 + i|b₀〉, which describes the entrance channel He⁺(²S) + H⁺₂(X̃¹Σ⁺_g).

We have

|ψ₀〉 = |a₀〉 + i|b₀〉 = |1〉|g₀(R)〉|0j₀(r)〉|j₀K₀〉, (1)

where |1〉 is the first diabatic electronic state;

g₀(R) = π^−1/4α^−1/2 exp[−(R − R₀)²/2α²]exp[−i(2μ_RE₀)^1/2(R − R₀)]; (2)

and |0j₀(r)〉|j₀K₀〉 are the H₂⁺(X̃¹Σ⁺_g) vibrational and rotational states, respectively. α is the Gaussian width parameter, and μ_R is the He⁺ + H₂ reduced mass,

The recursions are

|a₁〉 = Ĥ_s|a₀〉 − (1 − Ĥ_s²)^1/2|b₀〉, which is the first complex propagation, (3)

|a_n+2〉 = 2Ĥ_s|a_n+1〉 − |a_n〉, which is theother real Chebyshev propagation, (4)

where the square root in eqn (3) is evaluated as a Chebyshev expansion and eqn (4) is the standard Chebyshev propagation of a real WP.

The WP is absorbed at R > R_abs and r > r_abs using the Gaussians exp[−C^R_abs(R − R_abs)²] and exp[−C^r_abs(r − r_abs)²], in order to terminate propagation in the reactant Jacobi coordinates and in time: with the parameters in Table 1, the corresponding basis states are 7 053 300 basis states at K = 0, including two coupled electronic states. These parameters span the whole range of E_coll from 0.0002 to 1 eV, with ΔE_coll = 0.0001 or 0.001 eV below or above 0.2 eV, respectively.

Table 1 General parameters for the calculations. Distances are in bohr

Gaussian α, R0, and E0	0.2, 20, and 1.5 eV
R range and number of grid points	0.5–45 and 461
r range and number of grid points	0.5–15 and 153
Number of associated Legendre functions |jK〉	50 (j = even)
R and r absorption start at	30 and 11
R and r absorption strength	0.0005 and 0.005
Flux analysis at rf	10

---

The nonadiabatic initial-state-resolved reaction probabilities are calculated at the end through a time-to-energy Fourier transform and a flux analysis⁵⁸ on the |1〉 PES at r_f.

The reaction probabilities are obtained by the formula

where (1) and (2) are the diabatic electronic states; v′, j′, and K′ label open vibrational, rotational, and helicity states of the products, and S^Jp is the state-to-state symmetry adapted S-matrix at total angular momentum J and inversion parity p.

Coupled-channel calculations are performed up to J = 164, j₀ = 7, and E_coll = 1 eV, propagating WPs by kilo-steps (k-steps). After extensive calculations and checks, converged k-steps were found at the values shown in Table 2 and K_max to j₀ + 3.

Table 2 Parameters for the calculations of j_< = min(j₀, J) and K_max = j₀ + 3

j 0	K 0	J	k-Steps
0	0	0–24	480
		25–52	40
		53–84	20
		85–150	10
1	≤j<	0–22	480
		23–32	120
		33–42	70
		43–164	10
2	≤j<	0–19	200
		20–35	100
		36–51	50
		52–99	10
3	≤j<	0–20	160
		21–30	80
		31–40	40
		41–149	10
4	≤j<	0–20	140
		21–30	80
		31–50	40
		51–100	10
5	≤j<	0–18	120
		19–27	80
		28–38	40
		39–106	10
6	0	0–17	170
		18–27	90
		28–37	50
		38–102	10
7	0	0–18	170
		19–28	90
		29–36	50
		37–104	10

---

Note that in Table 2, a large number (up to k-step = 480) of propagation steps is required to obtain convergent (within a few percent) reaction probabilities at low collision energies (below 0.2 eV) for the first 20 J. As discussed in detail in subsequent sections, this is needed because of the presence of very sharp resonance peaks with unusually large lifetimes. They markedly increase the computational load, requiring wavepacket propagation for a large number of time steps to obtain the high accuracy of the resonance poles. In particular, the real part of the resonance energy pole (the energy peak position) converges reasonably fast, whereas the imaginary part (related to the heights and widths of the peaks) converges slowly. After J = 20, the resonance features sharply disappear, and the reaction probabilities converge in few k-steps (up to 10 for the larger partial waves). At a high collision energy, only background direct reactivity occurs, and convergence is achieved by a few propagation steps.

For the calculation of the integral cross sections (ICS) and more markedly the reactive rate constants that are the main target of the present work, convergence requirements are less strict, essentially because the resonance dynamics affects less these observables. This is, of course, due to the Boltzmann averaging that smears out the resonance contribution, delocalizing the convergence error, and to the higher relevance of excited reactant rotational states that, as shown in the next sections, are less sensitive to resonance dynamics.

Another important parameter that requires control to permit the calculations of highly reactant rotational states j₀ is the number of helicity states (K) included in WP. As shown in ref. 52, a nearly centrifugal sudden behavior is presented by this reaction, so for the ground rotational state, only K_max = 3 is sufficient to obtain convergence in the quantum rate constants. Increasing j₀ increases the number of reactant's helicity projections K₀ so that a larger number of K is coupled. However, our convergence test shows that the near conservation of the helicity quantum number also applies to higher rotational states, so K_max = j₀ + 3 is sufficient for the convergence of the ICS.

The total integral cross section is defined by

Here, the sum over J is cut at a convergence value J_max (dependent on j₀ and E_coll, see Table 2), while the sum over K₀ truncates to min(j₀, J). An additional considerable increase in the computational load occurs when high initial rotational states are calculated, increasing linearly with the number of WP to be propagated. Note that the sum over K₀ and p is a simple averaging of the helicity reaction probabilities . However, as shown in detail in the next sections for high rotational states there is a small dependence on K₀ so that high convergent rate constants can be obtained, eliminating these sums and the (2j₀ + 1) term in eqn (6) with a large reduction in the computational load. For this reason, calculations for the highest initial states (j₀ = 6 and 7) are performed only for K₀ = 0, as reported in Table 2.

The initial-state-resolved ICSs are then translationally averaged to obtain the initial-state-resolved rate constant k_j0(T)

where T is the temperature and k_B is the Boltzmann constant. Again, in numerical implementation, the integral limit must be cut at a finite numerical value, which is obtained by numerical convergence, and is dependent on j₀ and T. In these calculations, we found that an energy limit of 1 eV is sufficient to obtain initial state specific rate constants up to 1000 K at the claimed accuracy of one percent.

Initial state specific k_j0(T) are Boltzmann averaged to obtain para- and ortho-thermal rate constants and combined weighting by the nuclear spin multiplicity to obtain normal-H₂ rates

where ev/od is labels for even/odd j₀ and p^B is the H₂ Boltzmann population (see Fig. 2).

3. Results

3.1 Integral cross sections and initial state selected rate constants

In Fig. 3, initial-state-resolved ICSs σ_j0 as functions of the collision energy E_coll (see eqn (6)) at j₀ = 0–7 for the para- and ortho-molecular hydrogen cases are shown in panels (a) and (b), respectively.

Fig. 3 Initial-state-resolved ICS (σ_j0) as a function of collision energy (E_coll) (see eqn (6)) for the values of j₀ = 0–7 for the para-(a) and ortho-(b) hydrogens.

All j₀ curves have an overall similar behavior, with intense resonance features at low collision energies (under 200 meV) and a flat minimum of around 0.4 eV. Essentially, at low collision energies, the reactivity proceeds nearly exclusively by resonance dynamics; then, the background reactivity starts to increase moderately around about 0.5 eV. The curves show patterns of several resonance series, gradually becoming less intense, in larger numbers at lower values of j₀.

The important features exhibited by the ICS are as follows.

(1) The rotational energy of the reactants increases the background reactivity, so the ICSs of the largest ortho- and para-rotational states are about one order of magnitude larger than that of the ground state ICS. This means that the yield of this process is not only sensitive to the reactant vibrational state (see ref. 52) but also is sensitive to prominent rotational effects.

(2) The role of resonances in reactivity is present for all j₀ curves, but relevant differences appear. The most intense and narrowest resonance features are in the ground states of the reaction (j₀ = 0 and 1), while the resonance pattern for the excited rotational states is broader, and the resonance series moves at higher energies.

These aspects significantly affect the convergence and requirements for the calculation of the thermal rate constants shown by the parameters in Table 2. To illustrate the shift with j₀ of the resonance series, in Fig. 4 and 5, a blowup of up to 200 meV is shown in Fig. 3. For clarity, we report the resonance energy pattern of the initial-helicity-state-resolved ICS σ_j0K0 = 0 (see discussion of eqn (6)) using an energy grid of 0.1 meV, which can resolve most of the resonance peaks. As discussed below, for the K₀ helicity-averaged ICS σ_j0 (such as the ones shown in Fig. 3), some additional peaks occur essentially because of a slight dependence of the resonance energies on the parity inversion p, making the resonance patterns in the figure less resolved.

Fig. 4 Initial-helicity-state-resolved ICS (σ_j0K0) for K₀ = 0 as a function of collision energy (E_coll) (see eqn (6)) for j₀ = 0–3 at collision energies below 0.2 eV.

Fig. 5 Initial-helicity-state-resolved ICS (σ_j0K0) for K₀ = 0 as a function of collision energy (E_coll) (see eqn (6)) for j₀ = 4–7 at collision energies below 0.2 eV.

In this enlarged energy scale, we can observe that the resonance features below 50 meV do not appear to shift with the reactant rotational energy, while the other series move to higher collision energies. In Section 4, this resonance behavior is analyzed and rationalized. In the j₀ = 0 case, the low energy peaks and the resonance of the first group overlap, giving a considerable amplification of reactivity. For j₀ = 1 the splitting of the two series is incipient, fully developed from j₀ = 2, where the second shifting series is outside of the energy range shown.

As mentioned in Section 2, the resonance pattern plays a profound role in the convergence of the rate constants. The bottleneck for the calculations is that the resonances involved in this process have unusual long decaying lifetime, so that they are very difficult to be pinpointed by the wavepacket propagation methods, requiring to follow the reaction's evolution along very large collision times. Fortunately, translational averaged quantities, such as rate constants, have a slightly less strict convergence requirement because isolated resonances merge by Boltzmann averaging. Moreover, most of the resonance features occur in the low collision energy range, affecting the initial state selected rate constants k_j0 only below 300 K, where most of the excited reactants rotational states are not significantly populated (see Fig. 2). To obtain highly converged thermal rates of up to 1000 K (one of the goals of this work), the use of long wave-packet propagations is therefore required only for the ground rotational states of the ortho- and para-reactants (see Table 2). In addition, as shown in Fig. 4 and 5, the resonance lifetimes reduce significantly, increasing the energy of the rotational reactants so that the resonance peaks converge with a smaller number of time steps.

To show this behaviour, we compare the convergence with the number of time steps of the state specific rate constants k_j0 for j₀ = 0, 1, 3 and 5 in Fig. 6.

Fig. 6 State-specific rate constants (k_j0) for j₀ = 0, 1, 3 and 5 using 200, 300, 400 and 480 k-steps.

As shown, the convergence is slow for j₀ = 0 but rapidly increases, enhancing the rotational quantum number. For j₀ = 0, 480 k-steps are required to obtain a converged maximum at around 1%, but only 200 k-steps suffice to well converge at the j₀ = 5 rate. The more computationally demanding aspect of the whole calculation is therefore less so for excited rotational states.

Another important drawback of the proposed calculation is the reactive dependence on the initial helicity K₀, which makes the calculation very demanding for high rotational states, where 2j₀ + 1 initial WP must be propagated.

In Fig. 7, we show the helicity-resolved rate constants k_j0K0 for the j₀ = 5 case. We observe that the results from different helicity quantum numbers and/or parity inversions yield differences within 40%. The largest difference is for K₀ = 1, while the others are within 10–20%.

Fig. 7 State-specific helicity-resolved rate constants (k_j0K0) for j₀ = 5. (a) Curves for different K₀; dashed lines (with the respective color) are the results for the negative values of the inversion parity p. WPs were propagated for 200 k-steps. (b) Curves for K₀ = 0 at 120, 200 and 300 k-steps.

In Fig. 8, we compare K₀ = 0 helicity results with the K-averaged value used in ICS calculations for six values of j₀ (0–5). We observe that restricting the ICS calculation to only K₀ = 0, as discussed in the comment of eqn (6), gives accurate results for the kinetics of the reaction for the highly excited initial rotational states of the reactants, with a substantial reduction in computational cost. The result is poor for j₀ = 1 (about 15%), but the error decreases rapidly and becomes near the numerical accuracy of the calculation for j₀ = 5. This behavior is very important for extending the calculations to higher temperatures, which is more relevant in the study of chemical processes in the sun and in nebulae and exoplanet atmospheres.

Fig. 8 State-specific rate constants (k_j0) for j₀ = 0–5. Dashed lines are the corresponding curves obtained using only K₀ = 0.

Considering the high precision and low population of the rotational states involved (see Fig. 2) in the production runs, we used all the K₀ components until j₀ = 5. For j₀ = 6 and 7, only the component K₀ = 0 was included, producing thermal rate constants within our claimed numerical accuracy (1%).

3.2 Thermal rate constants

Fig. 9 reports all the eight rotationally selected rate constants k_j0 considered. Weighted by the Boltzmann populations illustrated in Fig. 2, they are employed to calculate the thermal rate constants up to 1000 K, as shown in Fig. 10. We can observe that the shape and intensity of the curves change significantly when the reactant rotational number varies.

Fig. 9 State-specific rate constants (k_j0) for j₀ = 0–7 as a function of temperature. Different colors label different reactant rotational states, as indicated in the graph. In the inset, an enlarged rate constant range is adopted to show the reactivity of the largest initial rotational states (j₀ = 6 and 7) considered in the calculations of thermal rate constants.

Fig. 10 Thermal rate constants (k) for the reaction of He⁺ with para (blue curve), ortho (red) and normal (black) H₂ as a function of temperature. For comparison, the reported initial state-selected rate constants k_j0 = 0 (blue) and 1 (red) curves are given, which are often used in applications.

The j₀ = 0 case is the only one that presents a maximum where the resonance reactivity is clearly demonstrated in the kinetics of the reaction. Already, for j₀ = 1, the maximum is very flat, and at j₀ = 2, it disappears. The overall reactivity initially decreases until it reaches the minimum at j₀ = 2, then steeply increases, so that the yield of the largest two j₀ (6 and 7 in the inset of Fig. 9) are about an order of magnitude larger. In other words, the rotational effects of the reactants are markedly relevant in the kinetics of this reaction, like the vibrational degree of freedom reported in ref. 52.

In Fig. 10, the thermal para, ortho, and normal rate constants (see eqn (8)) are compared with their ground states (j₀ = 0 and 1) k_j0, which are often used in applications. As observed by comparison with the black k_norm curve, the ground state approximation for the para-hydrogen is very poor, giving about double the yield and showing a decreasing behavior with T, while considering rotational effects, the reactivity increases at high temperatures. Some traces of the maximum around 150 K present in the ground state kinetic still survive in the black curve but appear as a flat shoulder until about 500 K when the rates start to increase linearly with T. However, the approximation works better using j₀ = 1 (dashed red curve), at least up to 400 K. Therefore, differences start to increase linearly so that the corresponding rates differ by more than 20% at 1000 K. The normal-hydrogen curve is very similar to that of ortho, except around 150 K, where a clear shoulder survives. We note that, although rotational effects increase moderately (about 30%), the rates at 1000 K, extrapolating the linear increase, could change the yield even by an order of magnitude at the high temperature (up to 5000 K), relevant in the model for exoplanet atmospheres as used in ref. 55 and 56.

4. Results

4.1 Reaction mechanism

As the main result of this work, resonances dominate the reaction mechanism, especially at low temperatures. To understand why and how it works, in Table 3, we report the geometries of the vdW and CI minima of the òA′ PES and the corresponding energies.

Table 3 vdW and CI minimum geometries. The zero for energy is taken at the minimum of the reactant's valley

	vdW	CI
R/a.u.	4.45	4.89
r/a.u.	1.42	2.18
γ/deg	90	90
Energy/eV	−0.082	1.344

---

Note that R_vdW and R_CI are very similar, as well as r_eq and r_vdW, while r_CI is larger. Moreover, γ_vdW and γ_CI coincide. The geometry of the vdW well and of the CI minimum is very similar, but the diatom internuclear distance must stretch to reach the CI minimum geometry. These features, as will be shown in the following, considerably affect deeply the reaction dynamics.

In Fig. 11, we show one dimensional cut in the C_2v configuration for the à and B̃ adiabatic PESs, with different values of r (namely at r_eq, r_CI and r_vdW of Table 3). From the figure, it can be observed that at r_eq of H₂, the CI seam is more than five eV higher in energy with respect to the collision energies involved. However, as r increases up to r_CI, the CI seam moves rapidly to lower energies until it reaches its minimum at about 1.3 eV. The cuts at r_vdW are nearly indistinguishable from the ones at r_eq even in the inset where an enlarged energy scale is adopted to exhibit the shape of the vdW well. This vdW well, localized in the reactant valley, is easily accessible to the system. The inset also shows that the vdW well is deep enough (see Table 1) to support several bound and quasi-bound states, arguably responsible for giving the resonance series observed.

Fig. 11 Cuts of the adiabatic potential energy surfaces of à (black lines) and B̃ (red lines) as a function of R taken at the values of r equal to the vdW minimum, equilibrium geometry of H₂, and CI seam minimum (solid, dashed and dashed-dotted lines, respectively). The zero for energy is taken at the minimum of the reactant's valley. In the inset, the energy scale is enlarged to show the vdW well.

In panels (a) and (b) of Fig. 12, 1D cuts over the Jacobi angle γ at CI and vdW geometries reported in Table 3 are plotted to show the anisotropy of the features (the energy scales are the same in both panels).

Fig. 12 Cuts of à (black lines) and B̃ (red lines) adiabatic PESs of ref. 24: in panels (a) and (b) as function of γ for r equal to that of the VdW minimum, and to the CI minimum (see Table 3), solid and dashed-dotted lines respectively. In panel (c) a cut of the adiabatic PESs in C_2v configuration vs. r is plotted. Note the different energy scale in panel (c). The energy zero is at the minimum of the reactant's valley.

Both the features present a marked anisotropy (note that the same energy scale is used in both panels), suggesting that the vdW well plays an important role in driving the system to the C_2v configuration where the avoided crossing between à and B̃ reaches its minimum. A similar stereodynamic effect of the vdW well is also observed in the reaction mechanism of the Cl + HD⁶⁷ and in the F + HD⁶⁸ reactions.

In panel (c), a cut along r in the C_2v configuration at R_CI and R_vdW is plotted. The figure shows that the CI between the two PESs acts as a high barrier in r that needs to be crossed by tunneling to have the H₂ dissociation. To have a non-adiabatic crossing, a stretch of the hydrogen molecular bond is required to favor charge transfer and diatomic dissociation by tunneling.

According to this reaction mechanism, although the crossing between the two electronic PESs is lower at large r, it remains at higher energies with respect to those of the system, especially in the low collision energy range. Therefore, to have efficient tunneling, sufficient time is required to convert the translational collision energy into the vibrational diatomic degree of freedom. This is favored by the formation of resonance metastable states, especially if they are populated in high vibrational diatomic states.

As is well known (see, e.g., ref. 69–71), CI seams can influence dissociation dynamics even when the dissociating species is well below the CI. Already since the dawn of the quantum theory, von Neumann and Wigner⁷² discovered potentials able to support spatially discrete states with energies within the continuum. In 2003, Cederbaum, Moseyev et al. published a nonadiabatic bi-dimensional analytic model⁷³ where bound states embedded in the continuum were found. Note the clear similarity between Fig. 1 of ref. 73 and panel (c) in Fig. 12. As discussed by Althorpe,⁷⁴ this happens because of the destructive interference between two types of quantum ‘trajectories’ circling around the CI on the opposite sides. Because these two ‘trajectories’ have opposite geometrical phase effects, they interfere destructively, leading to a node in the wavefunction. Of course, bound states in the continuum cannot be observed in scattering experiments. However, if the system is perturbed and its symmetry is lowered, the bound states interact with the continuum and turn into extraordinarily long-lived resonances,⁷³ as observed in this dissociative process. As discussed in Section 2, GP effects are automatically included when the wavefunction propagation is performed in diabatic representation, as in the calculations presented here.

4.2 Resonance analysis

Because of the high relevance of resonances in the reaction mechanism of the title reaction, we further analyzed the resonance pattern of the reaction probabilities for several of the reactant rotational states investigated, mostly along the lines followed previously by us.

In Fig. 13, J = 0 reaction probabilities for the initial reactant rotational states j₀ of the ortho-hydrogen molecule are plotted as a function of the total energy E_tot in the energy range 0.25–0.82 eV in logarithmic scale. As can be observed, all curves exhibit a similar behavior with narrow oscillations near the reaction threshold and two strong resonance peaks (more than one order of magnitude larger than the scattering background) just below the opening of each reactant channel. The features near the reaction thresholds have no significant impact on the reactive observable (ICS and rates), so in this work their origins are not further investigated (arguably not associated to resonances). The strong resonance peaks just below the opening of a new reactive channel are instead responsible for the resonance features discussed in Section 3; their evolution with J is subsequently examined.

Fig. 13 J = 0 reaction probability for the initial reactant rotational states (j₀) of the ortho-H₂ as a function of total energy (E_tot) in the energy range of 0.25–0.83 eV. Resonance features analyzed in this work are labeled as R_1–6. Note that the logarithmic scale is used for reaction probabilities.

All curves show resonance peaks at exactly the same total energy, whereas when plotted as a function of the collision energy, the peaks shift systematically to higher energies as j₀ increases (see ICS discussion in Subsection 3.1). This is, of course, due to the quadratic dependence on j₀ of the rotational energy that moves the energy threshold of the next rotational excited level; a metastable state is thus trapped at a higher collision energy.

Deep investigations of vdW resonance effects in the reaction dynamics for the F + H₂ and F + HD reactions were presented in ref. 75–77. Therefore, accurate energy poles obtained by a systematic theory of resonance lifetimes using the Felix Smith⁷⁸ method were successfully compared with simple adiabatic models already used in ref. 79 (and reported in detail in ref. 80) for the assignment of quasi molecular quantum numbers to the resonance peaks found in reaction probabilities. The fast vibrational motion of the product diatomic molecule HF was decoupled and integrated to obtain potential curves as a function of the product Jacobi arrangement R at different values of the diatomic product rotational quantum number j′. These adiabatic potentials show wells at the vdW geometry supporting (quasi-) bound states, where incoming atom trapping generates long-lived resonance states. The resonance energies (uniquely related to the real part of the poles in the complex energy plane) directly reflect the roto-vibrational spectrum of the metastable states involved. The resonance width is related to their lifetime and to the imaginary part of the resonance poles. Resonance energies are fixed at the same total energy and that when they match the reactants’ energy, a resonance feature can appear (or not) in the reaction probabilities if reactant and product eigenstates are significantly coupled (or uncoupled) with the resonance wavefunction: mathematically, when the residue of the resonance pole is sufficiently large (small) for the specific transition under study.

Observing the resonance patterns shown in Fig. 13 and with the adiabatic model in mind, we can conclude (and confirm next) that the labeled doublet resonance states are vdW quasi-bound triatomic states trapped in the adiabatic curves of the higher excited rotational states of the reactants. The first peak of the doublets corresponds to ground vibrational asymmetric stretching, whereas the second is the first excited vibrational level, which is very close to escaping from the adiabatic well. However, the coupling of the resonance (the resonance residues) is relevant only for the first two or three near neighborhood excited rotational states, as can be clearly observed by comparing the intensities of the peaks for j₀ = 1 and 7 at the highest total energy. A similar behavior is also observed for the para-H₂J = 0 reaction probabilities.

In the case of the FH₂ system, the extraction of the resonance parameters from the quantum dynamical calculation is complicated by coupling with the background direct scattering; sophisticated Q-matrix⁷⁸ or Regge pole⁸¹ analysis (as the ones presented in ref. 82 and 83) was required to clearly characterize the resonance effects in the reaction observables. However, as shown for the cold and ultra-cold regimes,^80,84 when the resonance mechanism is dominant with respect to the background direct scattering and the resonance widths are small enough to avoid interaction between different metastable states, we are under isolated narrow resonance conditions.⁸⁵ In this regime, the energy resonance poles can be extracted directly by reaction probabilities, obeying the basic Lorentzian functionality:

where A and bk are energy independent parameters accounting for the coupling with the reactant state (the resonance residue) and for the direct background scattering, respectively, supposed slowing varying along the resonance width Γ. E_res and Γ are respectively the energy position and the half width at half maximum of the distribution, which are directly related to the resonance complex energy pole (E_r = E_res + iΓ/2) and to its lifetime (h/2πΓ), respectively.

To check the working of the formula in eqn (9) for the resonances of this reaction, we plot in Fig. 14, as an example, the results of the Lorentzian fits for the four main resonance peaks for the J = 0 reaction probabilities in the j₀ = 0 case in the neighborhood of the first four resonance features. As can be clearly observed in all cases, the numerical results are very well described by the analytical Lorentzian functionals, and this holds for most of the j₀ and J calculated. Both double peaks show a similar spacing of about 26 meV, suggesting that changes in the adiabatic wells of different reactant rotational states are small. In both cases, the ground resonance state of the two doublets has a larger residue (such as in the doublets of Fig. 13), but the excited resonance state is significantly narrower (longer living). The lifetimes of all resonances are of the order of magnitude of picoseconds (about 2.5 for R₁) even though the lifetime resonances of the second doublet (R₃ and R₄) are about one order of magnitude smaller and continue to decrease gradually, increasing j₀. As expected, the second doublet also shows a larger direct (no resonant) scattering, but the nearly constant background condition within 2–3 times the resonance width is confirmed to be a good approximation in a narrow energy range. We stress here again that this order of magnitude for the resonance lifetimes is unusual in reactive scattering and comes from the extraordinarily long living properties of the bound molecular states embedded in the continuum generated by the conical intersection.⁷³

Fig. 14 J = 0 Reaction probabilities for the ground reactant rotational state j₀ = 0 in the vicinity of the resonance energies. Solid lines are the results of the fittings by eqn (9) with the parameters reported in each panel, full black circles are the results of the dynamical calculations. Curves of different colors indicate different doublets (see text for details).

Reaction probabilities at J = 0 are exemplary for clearly showing the position of the triatomic vibrational states of the metastable triatomic state without congestion due to the rotational manifold. Because of the higher multiplicity of the larger J (see eqn (6)), they affect the reactive observables in a minor way. To gain a deeper understanding of the resonance effects on reactivity, it is therefore important to study the evolution of the patterns with the total angular momentum J.

In Fig. 15, we show the opacity functions (2J + 1)P^J of j₀ = 0 for the first four J in the collision energy range of R₃ and R₄. As can be observed, the resonance peaks of J = 0 shift to higher collision energies as J increases. However, each resonance is split into multiplets so that the number of peaks for J different from zero rapidly grows. The splitting of J = 0 resonance poles into multiplets has been observed and analyzed in ref. 76 and 77, and confirmed by Regge pole analysis in ref. 83 for the F + H₂ reaction and its isotopic variants⁸⁶. In ref. 76, the components of each multiplet were shown to correspond to the rotational multiplicity of the adiabatic well supporting the quasi-bound state and are due to the dependence on the helicity states K of the adiabatic wells. Further increasing the J values, each component of the multiplet shifts to J(J + 1). Shifting growth is proportional to the rotational constant of the triatomic resonance state and is specific to each helicity component of the multiplet.⁷⁶

Fig. 15 Opacity functions (2J + 1)P^J for the ground reactant rotational state j₀ = 0 in the vicinity of the resonance energies R₃ and R₄ (see Fig. 14) for the first four values of J (indicated by different colors specified in the plot) showing the split of the resonance poles.

In Fig. 16, all the resonance positions up to 200 meV for the j₀ = 0 reaction are plotted as a function of J(J + 1). Each component of the multiplet is observed to follow rigorously a linear behavior. The angular coefficient of each regression line gives the rotational constants of the resonance.

Fig. 16 Energy positions (open circles) of the resonance peaks for j₀ = 0 as a function of J(J + 1). Dashed lines are regression lines obtained by the linear fit of data. Different colors indicate different branches of the multiplets (see text for details).

At a J value specific to each pole branch, the residues become very small (when the quasi-bound state escapes from the adiabatic well), and no significant peaks appear in the reaction probabilities. Therefore, the number of resonance peaks in the reaction probabilities has a maximum near J = 5 and then decreases until around J = 20, where the reaction probabilities appear smooth and the peaks disappear. Therefore, only the components of the resonances that survive at high J significantly affect the ICS.

The figure also shows that at large values of J, new progressions appear. These “new” resonances (reported as black and red in the figure) are bound states for J = 0, but they escape from the adiabatic well, becoming resonances, when the rotational energy of the triatomic complex increases. Although the resonance progressions have similar behavior for different resonances, they split into different numbers of branches. The resonances R₁ and R₂ (Fig. 14) become triplets, while the J = 0 bound states do not split and become resonances escaping from the well. Otherwise, the higher energy resonances (R₃ and R₄) show a larger number of branches. These differences are due to the different rotational adiabatic wells trapping the doublets in the three cases (bound states, R_1–2 and R_3–4). Additionally, some other regularity behavior emerges regarding the persistence with J for the different resonances. The vibrational excited peaks (R₂ and R₄ in Fig. 14) disappear rapidly with J and have minor effects on the reactive observables. Instead, the lowest resonance of each doublet is the component that survives at the largest J, producing the largest impact on reactive observables.

Analyzing the angular coefficient of the different regression lines in Fig. 16, the first branch is always found to correspond to the lowest angular coefficient (higher inertia moment), as observed for the vdW resonances of the F + H₂ reaction.⁷⁶ For this reason, if the two resonances survive enough, it is possible that the straight line of the higher helicity component of the ground vibrational resonance crosses the lower component of the first vibrational multiplet. From Fig. 16, we can observe that in the j₀ = 0 case, it happens three times, both for the low energy (R_1–2) and high energy (R_3–4) multiplets. A similar crossing also occurs in the ground rotational resonance of the F + H₂ reaction, and as analyzed in ref. 87, it is responsible for experimentally detected resonance features⁸⁸ in the angular distribution⁸⁹ and intense oscillations in ICS.⁹⁰ As reported in ref. 87, when these crossings occur, the Riemann surfaces of the two resonances interact in the complex plane, providing a condition for a double pole. The lifetime of both resonances steeply increases, strongly enhancing reactivity.

To illustrate the effect on the reaction probabilities when such an interaction of resonances occurs, in Fig. 17, we plot the opacity functions (2J + 1) P^J from J = 9 to 14 at energies close to the first resonance crossing in Fig. 16 (yellow and brown lines) involving two branches of R₁ and R₂.

Fig. 17 Opacity functions (2J + 1) P^J for j₀ = 0 and J = 9–14 at the collision energies in the vicinity of the crossing between yellow and brown regression lines depicted in Fig. 16. Note the sudden decrease in the resonance width at J = 12, where crossing occurs.

As discussed in the comments in Fig. 16, both R₁ and R₂ are triplets with different vibrational levels in the same adiabatic well (j = 2). However, the highest branch of R₂ survives only for a few J, so at J of the first panel (J = 9), it is actually a doublet (indicated by arrows in the panel). At J = 10, the highest component of R₁ also starts to disappear, clearly decreasing its lifetime, and the second component of R₂ is barely visible at J = 11. However, the lowest component of R₂ still survives at J = 11 although it starts its decay. At this J, the highest component of R₁ is just the small peak a few meV below the first peak of R₂. At J = 12, the two resonances cross, producing the single narrow peak shown in the left panel of Fig. 17. At a higher J, just one significant peak (although not Lorentzian in shape) survives, following the regression line of R₁. Such crossing, therefore, significantly enhances the resonance reactivity, and it is one of the reasons that makes resonance effects larger for the ground rotational state in the reactive observables, as discussed in Section 3.

To observe how the pattern of resonance series varies for the excited rotational states of the reactants, in Fig. 18, we plot the example of all the resonance peaks in j₀ = 1 (K₀ = 0, 1; p = 1 and −1) reaction probabilities in the collision energy range for the first two resonances R₁ and R₂ (see Fig. 14). In agreement with the adiabatic model, the two resonances split now into quadruplets (quasi-bound states supported by the j₀ = 3 adiabatic well), while the resonances as bound states for J = 0 split in doublets. The collision energies of the appearance of resonances increase, as discussed before, because of the widening of the energy gap of the rotational states that supports the resonances. The figure clearly shows that the position of the resonances is independent of the reactant helicity number K₀ (open and full circles). However, the resonance residues can vary, so for some K₀, the resonance peaks cannot be distinguished by scanning the reaction probabilities, and some integer branches disappear; this is the case of K₀ = 0 (empty circles), where the lowest branch is not present. Additionally, some branches of different helicity or parity survive better than others at higher J, but the differences are small. However, it is confirmed again, as depicted in Fig. 16, that the excited vibrational states (R₂ and R₄) survive for only a few J. We note that the parity inversion p significantly changes the position of resonances, so for K₀ larger than zero, all branches with K differ from the zero shift for the odd parity (open triangles). The ICS for j₀ different from zero (see Fig. 3) shows a larger number of resonance peaks; however, the parity split can be barely observed in ICS because even parity is dominant, amplified also by the K₀ = 0 probabilities involving, of course, only even parity.

Fig. 18 Energy positions of the resonance peaks for j₀ = 1 as a function of J(J + 1) up to 0.13 eV. Meaning of symbols regarding different values of K₀ and parity inversions p are specified in the graph. Note that the larger number of peaks with respect to Fig. 16 is due to the higher resonance splitting involved.

A similar detailed analysis of the resonance energies is also carried out for other reactant rotational states, and all the features and trends reported above are confirmed. All the opacity function data are available for further analysis upon request from the authors.

5. Summary, conclusions and perspectives

In this paper, we reported an extensive study of the kinetics of the dissociative charge transfer process of the helium cation with the hydrogen molecule, from few up to 10³ Kelvin.

Quantum thermal rate constants are calculated for the ortho, para and normal H₂ cases (see Fig. 10). Rotational effects show moderate variation in the reaction thermal kinetics, within a factor of 2 compared with j₀ = 0, resulting in the range of temperatures considered. However, the large kinetic enhancement for high j₀ shown in Fig. 9 supports changes of about one order of magnitude at higher temperatures, where the highest rotational states are more populated.

Approximations limiting the reactivity of this system to the ground state are therefore not adequate for estimating its relevance in astrophysical models, which also assume a different (decreasing) behavior as a function of T. At least in the investigated range better should be to focus on the rates for j₀ = 1 (the ground state of ortho hydrogen). The large rotational effect could also be the reason for the unexpected discrepancies (by one order of magnitude) with the available experimental data in ref. 91 possibly sensitive to the thermalization of efficient reactants. An additional reason for the discrepancies could be that the non-adiabatic couplings considered, based on the old ab initio data of ref. 22, are too weak. As cited in Section 1, new excited PESs and couplings have been published²⁶ recently, so it should be interesting to check whether updated PESs can fill this gap.

The dominant reaction mechanism at the temperatures investigated is deeply influenced by a small van der Waals well in the entrance channel of the reaction (Table 3). This feature produces extraordinarily long-lived quasi-bound states embedded in the continuum,⁷³ generated by the conical intersection when geometric phase effects are considered. This trapping favors the energy transfer in the vibrational diatomic mode, which permits the elongation required to reach the configuration of the conical intersection minimum and to penetrate by tunneling the thin and high barrier (see Fig. 12, panel c). The van der Waals well also influences the direct background reactivity by driving the process toward the perpendicular configuration, where the CI shows its minimum (see Fig. 12 panels a and b).

Exploiting the extremely long resonance lifetime, of the order of pico-seconds, we obtained the real part of the resonance poles with good accuracy, better characterizing the prominent resonance series shown in the ICS of Fig. 13, simply by monitoring the positions of the peaks in the reaction probabilities as a function of the collision energy. As shown in Fig. 14, the J = 0 resonance pattern is essentially composed of couples of vibrational states trapped in the adiabatic wells of the rotational states of the reactants, which split into multiplets,⁷⁶ increasing the total angular momentum following the multiplicity of the excited rotational states supporting the resonance (see Fig. 16). The excited vibrational states are closer to escaping from the wells but, in some cases, survive enough with J, so lower branches can cross the trajectories of the vibrational ground resonance states, giving a strong enhancement⁸⁷ of the resonance features (see Fig. 17).

Although traces of the complete rotational excited resonance progressions are presented in Fig. 13, only the first two or three series have significant pole residues to emerge in the integral cross sections. The quadratic increase in the reactant rotational energy shifts the position of the collision energies, where these resonances occur to higher values, increasing the collision energy gap between two consecutive doublets. Remarkably, in the j₀ = 0 case, the two resonance patterns overlap, producing a clear maximum in the rates (Fig. 9). For the other initial rotational states, the resonance contribution is smeared out in a larger range of T and progressively less intense. Therefore, resonance dynamics gradually become less relevant in the kinetics of excited reactant rotational states. This may permit feasible extensions of the present calculations at a higher range of T, as demanded in the atmosphere modeling of the Neptune-like exoplanets.⁵⁵ In Fig. 8, we exhibited complications for the wave packet methods in requiring the repetition of the propagation for each K₀, and the small dependence by K₀ of the rate for high excited reactant rotational levels found here may favor accurate approximations of the rate constant behavior. Consideration of at least one vibrational excited reactant's state because of the enhancement of one order of magnitude of the reaction rates shown in ref. 52 should suffice to account for the increase in temperature. Such extension and other applications are in progress.

The relevance of this system as a benchmark in modern chemical kinetics and dynamics continues to be a main concern of our research. Recent advances⁹² of a basic renormalization approach have shown perspectives for rationalizing key fundamental parameters controlling elementary chemical reactions, in particular, for ionic processes.

Conflicts of interest

There is no conflict to declare.

Data availability

State specific rate constants for  j₀ = 0–7 as a function of temperature, are given in Tabular form in the supplementary information (SI). Supplementary information is available. See DOI: https://doi.org/10.1039/d5cp04927j.

Acknowledgements

This work is a tribute to our eminent colleague and friend N. Sathyamurthy, to whom this special issue is dedicated on his 75th birthday. We are very grateful to A. Garcìa Muñoz (Paris-Saclay) for scientific discussions and for revising a draft of the manuscript. We acknowledge the CINECA award under the ISCRA initiative for ensuring the availability of high-performance computing resources and support.

References

1. D. Galli and F. Palla, The dawn of chemistry, Annu. Rev. Astron. Astrophys., 2013, 51, 163.
2. S. Lepp, P. C. Stancil and A. Dalgarno, Atomic and molecular processes in the Early Universe, J. Phys. B, 2002, 35, R57.
3. E. T. Young, et al., Early science with SOFIA, the stratospheric observatory for infrared, Astrophys. J., 2012, 749, L17.
4. S. Heyminck, et al., GREAT: the SOFIA high-frequency heterodyne instrument, Astron. Astrophys., 2012, 542, L1.
5. J. M. Moorhead, R. P. Lowe, J.-P. Maillard, W. H. Wehlau and P. F. Bernath, Search for HeH⁺ in NCG7027, Astrophys. J., 1988, 326, 899.
6. R. Güsten, H. Wiesemeyer, D. Neufeld, K. M. Menten, U. U. Graf, K. Jacobs, B. Klein, O. Ricken and C. Risacher, Astrophysical detection of the helium hydride ion HeH⁺, Nature, 2019, 568, 357.
7. S. Vranckx, J. Loreau, M. Desouter-Lecomte and N. Vaeck, Determination of photodissociation and radiative association cross sections from the same time-dependent calculation, J. Phys. B., 2013, 46, 155201.
8. S. Bovino, M. Tacconi, F. A. Gianturco and D. Galli, Ion chemistry in the early Universe: revisiting the role of HeH⁺ with new quantum calculations, Astronom. Astrophys., 2011, 529, A140.
9. E. D. S. Courtney, R. C. Forrey, R. T. McArdle, P. C. Stancil and J. F. Babb, Comprensive chemistry of HeH⁺ in the Early Universe, Astrophys. J., 2021, 919, 70.
10. S. Adhikari, M. Baer and N. Sathyamurthy, HeH₂⁺: structure and dynamics, Int. Rev. Phys. Chem., 2022, 41, 49.
11. P. J. Kuntz, Use of the method of diatomics-in-molecules in fitting ab initio potential surfaces: the system HeH₂⁺, Chem. Phys. Lett., 1972, 16, 581.
12. D. G. Hopper, The electronic structure of HeH₂⁺, Int. J. Quantum Chem., 1978, 14, 305.
13. T. Joseph and N. Sathyamurthy, Three-dimensional quasiclassical trajectory study of the reaction He + H₂⁺ → HeH⁺ + H on an accurate ab initio potential-energy surface, J. Chem. Phys., 1984, 80, 5332.
14. A. Aguado and M. A. Paniagua, New functional form to obtain analytical potentials of triatomic molecules, J. Chem. Phys., 1992, 96, 1265.
15. M. F. Falcetta and P. E. Siska, The interaction between He and H₂⁺: anisotropy, bond length dependence and hydrogen bonding, Mol. Phys., 1999, 97, 117.
16. V. Špirko and W. P. Kraemer, Ab Initio Predicted Rotation-Vibration Energy Levels of HeH₂⁺, J. Mol. Spectrosc., 1995, 172, 265.
17. P. Palmieri, C. Puzzarini, V. Aquilanti, G. Capecchi, S. Cavalli, D. De Fazio, A. Aguilar, X. Giménez and J. M. Lucas, Ab initio dynamics of the He + H₂⁺ → HeH⁺ + H reaction: a new potential energy surface and quantum mechanical cross-sections, Mol. Phys., 2000, 98, 1835–1849.
18. W. Xu, X. Liu, S. Luan, Q. Zhang and P. Zhang, An ab initio potential energy surface of the He + H₂⁺ → HeH⁺ + H reaction, Chem. Phys. Lett., 2008, 464, 92.
19. C. N. Ramachandran, D. De Fazio, S. Cavalli, F. Tarantelli and V. Aquilanti, Revisiting the potential energy surface for the He + H₂⁺ → HeH⁺ + H reaction at the full configuration interaction level, Chem. Phys. Lett., 2009, 469, 26.
20. J. Sahoo, D. Bossion, T. Gonzalez-Lezana, D. Talbi and Y. Scribano, Low temperature dynamics of the H + HeH⁺ → H₂⁺ + He reaction: on the importance of the long-range interaction, J. Chem. Phys., 2024, 161, 144312.
21. F. Schneider and L. Zülicke, Approximate diatomics-in-molecules potential energy surfaces and non-adiabatic coupling for He⁺ + H₂, Chem. Phys. Lett., 1979, 67, 491.
22. D. R. Mclaughlin and D. L. Thompson, HeH₂⁺. Ground- and lower excited-state discrete ab initio electronic potential-energy surfaces for doublet, J. Chem. Phys., 1979, 70, 2748.
23. R. J. Furlan, G. Bent and A. Russek, Ab initio (HeH₂)⁺ energy surfaces and nonadiabatic couplings between them, J. Chem. Phys., 1990, 93, 6676.
24. A. Aguado, C. Suárez and M. Paniagua, Accurate fit of the two lowest excited-state potential-energy surfaces for doublet HeH₂, J. Chem. Phys., 1993, 98, 308.
25. V. Sidis, Diabatic excited states of the (HeH₂) + molecular ion for the charge exchange-excitation reaction: He⁺ + H₂ → HeH⁺ + H, Chem. Phys., 1996, 209, 313.
26. K. Naskar, S. Ghosh, S. Adhikari, M. Baer and N. Sathyamurthy, Beyond Born–Oppenheimer Constructed Diabatic Potential Energy Surfaces for HeH₂⁺, J. Phys. Chem. A, 2023, 127, 3832.
27. N. Sathyamurthy, R. Rangarajan and L. M. Raff, Reactive scattering calculations on a spline fitted ab initio surface: the He + H₂⁺ (ν = 0, 1, 2) → HeH⁺ + H reaction, J. Chem. Phys., 1976, 64, 4606.
28. T. Joseph and N. Sathyamurthy, Dynamics of the endothermic reaction He + H₂⁺ → HeH⁺ + H on an accurate ab initio potential-energy surface, J. Chem. Phys., 1987, 86, 704.
29. F. Esposito, C. M. Coppola and D. De Fazio, Complementarity between quantum and classical mechanics in chemical modeling. The H + HeH⁺→ H₂⁺ + He reaction: a rigourous test for reaction dynamics methods, J. Phys. Chem. A, 2015, 119, 12615.
30. B. Lepetit and J. M. Launay, Quantum mechanical study of the reaction He + H₂⁺ → HeH⁺ + H with hyperspherical coordinates, J. Chem. Phys., 1991, 95, 5159.
31. V. Aquilanti, G. Capecchi, S. Cavalli, D. De Fazio, P. Palmieri, C. Puzzarini, A. Aguilar, X. Giménez and J. M. Lucas, The He + H₂⁺ reaction: a dynamical test on potential energy surfaces for a system exhibiting a pronounced resonance pattern, Chem. Phys. Lett., 2000, 318, 619.
32. D. De Fazio, M. de Castro-Vitores, A. Aguado, V. Aquilanti and S. Cavalli, The He + H₂⁺ → HeH⁺ + H reaction: ab initio studies of the potential energy surface, benchmark time-independent quantum dynamics in an extended energy range and comparison with experiments, J. Chem. Phys., 2012, 137, 244306.
33. D. De Fazio, The H + HeH⁺ → He + H₂⁺ reaction from the ultra-cold regime to the three-body breakup: exact quantum mechanical integral cross sections and rate constants, Phys. Chem. Chem. Phys., 2014, 16, 11662.
34. B. Maiti, C. Kalyanaraman, A. N. Panda and N. Sathyamurthy, Reaction probabilities and reaction cross sections for three-dimensional He + H₂⁺(ν) collisions: a time-dependent quantum mechanical study, J. Chem. Phys., 2002, 117, 9719.
35. A. K. Tiwari and N. Sathyamurthy, Preferential scattering of one isotopomer over another in (He, HD⁺) collisions, Chem. Phys. Lett., 2005, 414, 509.
36. A. N. Panda and N. Sathyamurthy, Time-dependent quantum mechanical wave packet study of the He + H₂⁺(v,j) → HeH⁺ + H reaction, J. Chem. Phys., 2005, 122, 054304.
37. X. Tang, C. Houchins, K.-C. Lau, C. Y. Ng, R. A. Dressler, Y.-H. Chiu, T.-S. Chu and K.-L. Han, A time-dependent wave packet quantum scattering study of the reaction HD⁺ (v = 0–3; j₀ = 1) + He → HeH⁺ (HeD⁺) + D(H), J. Chem. Phys., 2007, 127, 164318.
38. P. Gamallo, S. Akpinar, P. Defazio and C. Petrongolo, Quantum Dynamics of the Reaction H(²S) + HeH⁺(X¹Σ⁺) → H₂⁺(X²Σ_g⁺) + He(¹S) from Cold to Hyperthermal Energies: Time-Dependent Wavepacket Study and Comparison with Time-Independent Calculation, J. Phys. Chem. A, 2014, 118, 6451.
39. W. A. Chupka and M. E. Russell, Photoionization Study of Ion–Molecule Reactions in Mixtures of Hydrogen and Rare Gases, J. Chem. Phys., 1968, 49, 5426.
40. W. A. Chupka, J. Berkowitz and M. A. Russell. A Study of Some Reactions of H₂⁺ in Selected Vibrational States. 6th International 2 Conference on the Physics of Electronic and Atomic Collisions: ICPEAC VI, 1969, p 71.
41. D. van Pijkeren, E. Boltjes, J. van Eck and A. Niehaus, Vibrational-state-selected ion-molecule reaction cross sections at thermal energies, Chem. Phys., 1984, 91, 293.
42. T. Turner, O. Dutuit and Y. T. Lee, The effects of collision energy and vibrational excitation on H₂⁺, HD⁺ + He reactions, J. Chem. Phys., 1984, 81, 3475.
43. M. Baer, S. Suzuki, K. Tanaka, I. Koyano, H. Nakamura, Z. Herman and D. J. Kouri, Phys. Rev. A: At., Mol., Opt. Phys., 1986, 34, 1748.
44. X. N. Tang, H. Xu, T. Zhang, Y. Hou, C. Chang, C. Y. Ng, Y. Chiu, R. A. Dressler and D. J. Levandier, A pulsed-field ionization photoelectron secondary ion coincidence study of the H₂⁺(X,v⁺ = 0–15, N⁺ = 1) + He proton transfer reaction, J. Chem. Phys., 2005, 122, 164301.
45. W. P. Kraemer, V. Špirko and O. Bludsky, Bound and low-lying quasi-bound rotation-vibration energy levels of the ground and first excited electronic states of HeH₂⁺, Chem. Phys., 2002, 276, 225.
46. F. Mrugala, V. Špirko and W. P. Kraemer, Radiative association of HeH₂⁺, J. Chem. Phys., 2003, 118, 10547.
47. F. Mrugala and W. P. Kraemer, Radiative charge transfer in He⁺ + H₂ collisions in the milli- to nano-electron-volt range: a theoretical study within state-to-state and optical potential approaches, J. Chem. Phys., 2013, 138, 104315.
48. T. Szidarovszky and K. Yamanouchi, Photodissociation dynamics of weakly bound HeH₂⁺ in intense light fields, Phys. Rev. A, 2016, 94, 063405.
49. S. Ravi, S. Mukherjee, B. Mukherjee, S. Adhikari, N. Sathyamurthy and M. Baer, Non-adiabatic coupling as a frictional force in (He, H, H) + dynamics and the formation of HeH₂⁺, Mol. Phys., 2021, 119, e1811907.
50. F. Aguillon, Semi-classical coupled wavepacket study of the dissociative charge exchange He⁺ + H₂ → He + H + H⁺, Chem. Phys. Lett., 1994, 222, 69.
51. F. Aguillon, A new treatment of nonadiabatic dynamics: application to the determination of the He⁺ + H₂ →He + H + H⁺ differential cross section, J. Chem. Phys., 1998, 109, 560.
52. D. De Fazio, A. Aguado and C. Petrongolo, Non-adiabatic quantum dynamics of the dissociative charge transfer He⁺ + H₂ → He + H + H⁺, Front. Chem., 2019, 7, A249.
53. K. Naskar, S. Ravi, S. Adhikari, M. Baer and N. Sathyamurthy, Coupled three-dimensional quantum mechanical wave packet study of proton transfer in H₂⁺ + He collisions on accurate ab initio two-state diabatic potential energy surfaces, J. Chem. Phys., 2023, 159, 034302.
54. A. K. Gupta, V. Dhindhwal, M. Baer, N. Sathyamurthy, S. Ravi, S. Mukherjee, B. Mukherjee and S. Adhikari, Non-adiabatic coupling and conical intersection(s) between potential energy surfaces for HeH₂⁺, Mol. Phys., 2020, 118, e1683243.
55. A. Garcìa Muñoz, Modelling helium in exoplanet atmospheres. A revised network with photoelectron-driven processes, Astronom. Astrophys., 2025, 689, A199.
56. A. Garcìa Muñoz, D. De Fazio, D. J. Wilson and K. France, Vibrationally excited H₂ mutes the He I triplet line at 1.08 µm on warm exo-Neptunes., Astronom. Astrophys., 2025, 704, L18.
57. S. K. Gray and G. G. Balint-Kurti, Quantum dynamics with real wave packets, including application to three-dimensional (J = 0) D + H2 → HD + H reactive scattering, J. Chem. Phys., 1998, 108, 950.
58. A. J. H. M. Meijer, E. M. Goldfield, S. K. Gray and G. G. Balint-Kurti, Flux analysis for calculating reaction probabilities with real wave packets, Chem. Phys. Lett., 1998, 293, 270.
59. P. Gamallo, A. Zanchet, F. J. Aoiz and C. Petrongolo, Non-adiabatic quantum dynamics of the electronic quenching OH(A²Σ⁺) + Kr, Phys. Chem. Chem. Phys., 2020, 22, 17091.
60. P. Defazio, P. Gamallo and C. Petrongolo, Nonadiabatic dynamics of O(¹D) + N₂(X¹Σ⁺_g) → O(³P) + N₂(X¹Σ⁺_g) on three coupled potential surfaces: Symmetry, Coriolis, spin–orbit, and Renner–Teller effects, J. Chem. Phys., 2012, 136, 054308.
61. C. Petrongolo, Nonadiabatic theory of triatomics: general formalism and application to Renner–Teller and conical-intersection effects, J. Chem. Phys., 1988, 89, 1297.
62. C. Xie, C. L. Malbon, D. R. Yarkony, D. Xie and H. Guo, Signatures of a Conical Intersection in Adiabatic Dissociation on the ground electronic state, J. Am. Chem. Soc., 2018, 140, 1986.
63. M. V. Berry, Quantal Phase Factors Accompanying Adiabatic Changes, Proc. R. Soc. A, 1984, 392, 45.
64. C. A. Mead, The Molecular Aharonov-Bohm Effect in Bound-States, Chem. Phys., 1980, 49, 23.
65. C. A. Mead, The geometric phase in molecular systems, Rev. Mod. Phys., 1992, 64, 51.
66. H. Guo, Quantum dynamics of complex-forming bimolecular reactions, Int. Rev. Phys. Chem., 2012, 31, 1.
67. D. Skouteris, D. Manolopoulos, W. Bian, H.-J. Werner, L.-H. Lay and K. Liu, van der Waals interactions in the Cl + HD reaction, Science, 1999, 286, 1713.
68. D. Skouteris, D. De Fazio, S. Cavalli and V. Aquilanti, Quantum Stereodynamics for the Two Product Channels of the F + HD Reaction from Complete Scattering Matrix in the Stereodirected Representation, J. Phys. Chem. A, 2009, 113, 14807.
69. S. Henshaw and A. S. Izmaylov, Topological origins of bound states in the continuum for systems with conical intersections, J. Phys. Chem. Lett., 2018, 9, 146.
70. C. Xie and H. Guo, Photodissociation of phenol via nonadiabatic tunneling: comparison of two ab initio based potential energy surfaces, Chem. Phys. Lett., 2017, 683, 222.
71. C. Xie, J. Ma, X. Zhu, D. R. Yarkony, D. Xie and H. Guo, Nonadiabatic Tunneling in Photodissociation of Phenol, J. Am. Chem. Soc., 2016, 138, 7828.
72. J. von Neumann and E. Wigner, Über Merkwürdige Diskrete Eigenwerte, Phys. Z., 1929, 30, 465.
73. L. S. Cederbaum, R. S. Friedman, V. M. Ryaboy and N. Moiseyev, Conical Intersections and Bound Molecular States Embedded in the Continuum, Phys. Rev. Lett., 2003, 90, 013001.
74. S. C. Altorphe, General explanation of geometric phase effects in reactive systems: unwinding the nuclear wave function using simple topology, J. Chem. Phys., 2006, 124, 084105.
75. V. Aquilanti, S. Cavalli, A. Simoni, A. Aguilar, J. M. Lucas and D. De Fazio, Lifetime of reactive scattering resonances: Q-matrix analysis and angular momentum dependence for the F + H₂ reaction by the hyperquantization algorithm, J. Chem. Phys., 2004, 121, 11675.
76. V. Aquilanti, S. Cavalli, D. De Fazio, A. Simoni and T. V. Tscherbul, Direct evaluation of the lifetime matrix by the hyperquantization algorithm: narrow resonances in the F + H₂ dynamics and their splitting for nonzero angular momentum, J. Chem. Phys., 2005, 123, 054314.
77. D. De Fazio, S. Cavalli, V. Aquilanti, A. A. Buchacenko and T. V. Tscherbul, On the Role of Scattering Resonances in the F + HD Reactions Dynamics, J. Phys. Chem. A, 2007, 111, 12538.
78. F. T. Smith, Lifetime matrix in Collision Theory, Phys. Rev., 1960, 118, 349; F. T. Smith, Lifetime matrix in Collision Theory, Phys. Rev., 1960, 119, 2098.
79. J. F. Castillo, D. E. Manolopoulos, K. Stark and H.-J. Werner, Quantum mechanical angular distributions for the F + H₂ reaction, J. Chem. Phys., 1996, 104, 6531.
80. D. De Fazio, S. Cavalli and V. Aquilanti, Benchmark Quantum Kinetics at Low Temperatures toward Absolute Zero and Role of Entrance Channel Wells on Tunneling, Virtual States, and Resonances: The F + HD Reaction, J. Phys. Chem. A, 2020, 124, 12.
81. T. Regge, Introduction to complex angular momenta, Il Nuovo Cimento, 1959, 14, 951.
82. C. Xiahou, J. N. L. Connor, D. De Fazio and D. Sokolovski, A single resonance Regge pole dominates the forward-angle scattering of the state-to-state F + H₂→ FH + H reaction at E trans= 62.09 meV, Phys. Chem. Chem. Phys., 2024, 26, 3647.
83. D. Sokolovski, D. De Fazio and E. Akhmatskaya, A Transition State Resonance Radically Reshapes Angular Distributions of the F + H₂ → FH(v_f = 3) + H Reaction in the 62–102 meV Energy Range, ACS Phys. Chem. Au, 2025, 5, 219.
84. D. De Fazio, S. Cavalli and V. Aquilanti, Quantum dynamics and kinetics of the F + H₂ and F + D₂ reactions at low and ultra-low temperatures, Front. Chem., 2019, 7, 328.
85. J. R. Taylor, Scattering Theory, Wiley, New York, 1972, Cap. 13.
86. D. Sokolovski, E. Akhmatskaya, C. Echeverra-Arrondo and D. De Fazio, Complex angular momentum theory of state-to state integral cross sections: resonance effects in the F + HD → HF(v′ = 3) + D reaction, Phys. Chem. Chem. Phys., 2015, 17, 18577.
87. D. Sokolovski, S. Sen, V. Aquilanti, S. Cavalli and D. De Fazio, Interacting resonances in the F + H₂ reaction revisited: complex terms, Riemann surfaces and angular distributions, J. Chem. Phys., 2007, 126, 084305.
88. M. Qiu, et al., Observation of Feshbach resonances in the F + H₂ → HF + H reaction, Science, 2006, 311, 1440.
89. D. Sokolovski, D. De Fazio, S. Cavalli and V. Aquilanti, On the origin of the forward peak and back oscillations in the F + H₂(v = 0) → HF(v′ = 2) + H reaction, Phys. Chem. Chem. Phys., 2007, 9, 5664.
90. D. Sokolovski, D. De Fazio, S. Cavalli and V. Aquilanti, Overlapping resonances and Regge oscillations in the state-to-state integral cross sections of the F + H₂ reaction, J. Chem. Phys., 2007, 126, 121101.
91. R. Johnsen, A. Chen and M. A. Biondi, Dissociative charge transfer of He⁺ ions with H₂ and D₂ molecules from 78 to 330 K, J. Chem. Phys., 1980, 72, 3085.
92. V. H. Carvalho-Silva, F. O. Sanches-Neto, G. M. Leão, S. Cavalli, D. De Fazio, M. C. Nucci and V. Aquilanti, Renormalized chemical kinetics and benchmark quantum mechanical rates: activation energies and tunnelling transitivities for the reactions of fluorine atoms with H₂ and HD, Rend. Phys. Acc. Lincei, 2023, 34, 997.

---