Non-adiabatic Quantum Dynamics of the Dissociative Charge Transfer He++H2 → He+H+H+

We present the non-adiabatic, conical-intersection quantum dynamics of the title collision where reactants and products are in the ground electronic states. Initial-state-resolved reaction probabilities, total integral cross sections, and rate constants of two H2 vibrational states, v0 = 0 and 1, in the ground rotational state (j0 = 0) are obtained at collision energies Ecoll ≤ 3 eV. We employ the lowest two excited diabatic electronic states of HeH2+ and their electronic coupling, a coupled-channel time-dependent real wavepacket method, and a flux analysis. Both probabilities and cross sections present a few groups of resonances at low Ecoll, whose amplitudes decrease with the energy, due to an ion-induced dipole interaction in the entrance channel. At higher Ecoll, reaction probabilities and cross sections increase monotonically up to 3 eV, remaining however quite small. When H2 is in the v0 = 1 state, the reactivity increases by ~2 orders of magnitude at the lowest energies and by ~1 order at the highest ones. Initial-state resolved rate constants at room temperature are equal to 1.74 × 10−14 and to 1.98 × 10−12 cm3s−1 at v0 = 0 and 1, respectively. Test calculations for H2 at j0 = 1 show that the probabilities can be enhanced by a factor of ~1/3, that is ortho-H2 seems ~4 times more reactive than para-H2.


INTRODUCTION
Atomic Hydrogen and Helium are the dominant chemical species of the early Universe (Galli and Palla, 2013) and of the interstellar medium, and are easily ionized by cosmic rays. Therefore, these atoms and their ions are the astrochemical fundamental reactants (Lepp et al., 2002), together with ubiquitous photons, giving first simple diatoms as H 2 , H + 2 (Stancil et al., 1993), and HeH + (Zygelman et al., 1998) and then atom+diatom bimolecular collisions as He+H + 2 (De Fazio et al., 2012), H+HeH + (De Fazio, 2014; Gamallo et al., 2015), and He + +H 2 .
When all chemical species are in the ground electronic states, the He+H + 2 reaction is endothermic and rather slow, but the other two are exothermic by about ∼0.7 and 6.2 eV, respectively, as many ion+neutral astrochemical reactions (Herbst and Klemperer, 1973). H+HeH + is barrierless and gives quickly the He+H + 2 products, but the collision-induced dissociative charge transfer (DCT) He + ( 2 S)+H 2 (X 1 + g )→ He( 1 S)+H( 2 S)+H + is very slow at low collision energy E coll and at room temperature, if H 2 is in the ground vibro-rotational state. In fact, H+HeH + ↔He + +H 2 occurs on the ground potential energy surface (PES)X 2 A' of HeH + 2 , which is well-separated from the excited electronic species, but the low-E coll DCT involves the first two excited adiabatic electronic states A 2 A' andB 2 A' of HeH + 2 , which differ by a two-electron excitation and are coupled by a C 2v conical intersection (CI) (Preston et al., 1978;McLaughlin and Thompson, 1979). In the latter case the non-adiabatic coupling is weak and the lower, dissociative cone of the intersection seam gives rise to an adiabatic barrier with diabatic character, that is the DCT tends to follow the diabatic PESs without changing the electronic configuration. This strongly inhibits the reactivity that is associated with the tunneling through the barrier (Preston et al., 1978). This is schematically shown in the correlation diagram of Figure 1 for theÃ 2 A' andB 2 A' adiabatic and (1) 2 A 1 and (2) 2 B 2 diabatic PESs V of HeH + 2 , where all chemical species are in the ground electronic states, save H + 2 (A 2 + u ) that is unbound, and the energy is referred to the reactant minimum. This diagram is obtained from the ab initio Multi-Reference Configuration-Interaction results of McLaughlin and Thompson (1979), changing the label of the reactant diabatic state 2 A 1 from (5) to (1), and from the analytical fits of the associated diabatic PESs V 11 and V 22 by Aguado et al. (1993). In the scheme of Figure 1 we omit the ground adiabatic PES, well below the excited PESs and with too small non-adiabatic couplings. However, all three adiabatic PESs are strongly coupled by intense laser pulses (Szidarovszky and Yamanouchi, 2016) when the ground PES becomes populated (Schauer et al., 1989). A more complete description of the reaction dynamics in presence of electric and magnetic fields is out of the scope of the present article. Note in the figure the small ion-induced dipole minimum in the C 2v entrance channel at R(He-H 2 ) = 4.45 a 0 , r(H-H) = 1.42 a 0 , and V = −0.08 eV and the CI C 2v minimum at R = 4.89 a 0 , r = 2.18 a 0 , and V = 1.34 eV, owing to the intersection (Preston et al., 1978) between the H 2 (X 1 + g ) and H + 2 (A 2 + u ) curves at r = 2.19 a 0 and V = 1.28 eV.
Experimental studies on the He + +H 2 → He+H+H + DCT dynamics date back to 1955 (Stedeford and Hasted, 1955) and 1961(Giese and Maier, 1961, when the integral cross sections (ICSs) were measured at E coll = 4 eV and found <0.05 Å 2 . ICSs as functions of E coll were then measured in many works up to 1996 (Dhuicq et al., 1996), finding values from 0.01 up to ∼2 Å 2 , for E coll between 3 and 100 eV, by considering H 2 in the ground vibro-rotational state and all ground and excited open states of H. Below 3 eV, the H product is in the ground 1s state and the ICSs are so small and so difficult to measure that experimental values present large discrepancies (Reinig et al., 1992). However, the ICSs increase by one/two orders of magnitude if H 2 is excited by one vibrational quantum (Preston et al., 1978) or E coll grows up to ∼100 eV. Differential cross sections were also measured (Dhuicq et al., 1998) at E coll > 9 eV, that is for the formation of the H * excited product. Accordingly, small rate constants were observed (Johnsen et al., 1980), with values of (1.5 ± 0.15) × 10 −13 and (1.1 ± 0.1) × 10 −13 cm 3 s −1 at 78 and 330 K, respectively, that is nearly four order of magnitude lower than the Langevin (Gioumousis and Stevenson, 1958) estimates for ion+neutral barrierless and exothermic reactions. FIGURE 1 | Dissociative charge transfer (DCT) and collision induced dissociation (CID) of He + +H 2 . Schematic correlation diagram in eV of the adiabaticÃ 2 A ′ (red) andB 2 A ′ (blue) and of the diabatic (1) 2 A 1 and (2) 2 B 2 PESs. All reactants and products are in the ground electronic state, save H + 2 (A 2 + u , unbound).
A few works have also theoretically investigated the dynamics of the DCT collision since 1994, when Aguillon (1994) employed a semiclassical coupled wavepacket (WP) method and the analytical diabatic PESs of Aguado et al. (1993) for computing ICSs for ground and excited vibrational states v 0 of H 2 , up to v 0 = 4 and in the E coll range from 2 to 10 eV. He found total-ICS values from ∼0.001 Å 2 (v 0 = 0, E coll = 2 eV) to ∼1 Å 2 (v 0 = 3, E coll = 10 eV), in agreement with the most recent measurements (Dhuicq et al., 1996). In a subsequent work (Aguillon, 1998), an improved version of this semiclassical method was used for obtaining new ICS values and differential cross sections. Approximated theoretical models were also employed (Dhuicq et al., 1996) for explaining observed cross sections above 9 eV, when excited H * (n = 2) was produced. Finally, quantum infiniteorder-sudden cross sections were calculated (Martínez et al., 2000) for He + +H 2 → He+H * (n ≥ 2)+H + at E coll ≥ 10 eV, using the accurate diabatic representation of Aguado et al. (1993) and six more approximated diabatic electronic states (Sidis, 1996).
As far as we know, no further studies of the title reaction have been published and its collision dynamics below ∼2 eV is unknown. In particular, only semiclassical or approximated quantum theoretical studies were carried out, although both conical intersection and barrier tunneling are purely quantum effects. We thus here report a rigorous time-dependent quantum study of the DCT He + +H 2 (v 0 = 0,1)→ He+H+H + with all species in the ground electronic state, at thermal and hyperthermal collision energy up to 3 eV, using the diabatic PESs of Aguado et al. (1993) and WP and flux formalisms. In section Theory and Calculations we present the theoretical method and its numerical details. Section Collision Results reports initial-state-resolved total reaction probabilities, ICSs, and thermal rate constants. Finally, we present our conclusions in section Conclusions.

Potential Energy Surfaces and Non-adiabatic Coupling
As we said in the Introduction, the HeH + 2 adiabatic and diabatic electronic states relevant in the present work were obtained ab initio in McLaughlin and Thompson (1979), fitted analytically in Aguado et al. (1993), and they are schematically plotted in Figure 1. We label the adiabatic species and PESs byÃ 2 A' andB 2 A' and by V A and V B , respectively. As already discussed (Aguado et al., 1993), these states belong to the fully symmetric irreducible representation for linear (C ∞v ) and non-symmetric (C S ) geometries, while Figure 1 shows that they transform as A 1 /B 2 and B 2 /A 1 for perpendicular geometries (C 2v ), before/after the CI, respectively, which rules the title reaction. On the other hand, we label the associated diabatic electronic states and PESs by (1)A 1 and (2)B 2 and by V 11 , V 22 , and V 12 , respectively, where the third surface describes the CI non-adiabatic coupling in the diabatic representation.
In order to simultaneously describe these electronic species and take into account the CI, a fit of the diabatic PESs and coupling V 11 , V 22 , and V 12 was made in Aguado et al. (1993). The description of the adiabatic states is obtained as the eigenvalues , in which the interaction term V 12 must have the correct symmetry, being anti-symmetric with respect to the permutation of the H atoms and thus vanishing identically for equal He-H distances. This coupling term V 12 was fitted in Aguado et al. (1993) to reproduce the CI between theÃ 2 A' andB 2 A' adiabatic states.
The diabatic surfaces where fitted in Aguado et al. (1993) using the Aguado-Paniagua functional form (Aguado and Paniagua, 1992;Aguado et al., 1998) that expands the energy as a multidimensional permutationally invariant polynomial (Aguado et al., 1994(Aguado et al., , 2001 in Rydberg type variables (Rydberg, 1931) of the form ρ AB = R AB exp (−β AB R AB ), where A and B are two nuclei and R AB is their distance. For the interaction term, a simple expansion in Rydberg functions, The DCT is produced through the CI, as shown in Figure 2 using nuclear Jacobi coordinates R, r = R HH ′ , and γ , at R = 4 a 0 and γ = 0. The top panel shows the adiabaticÃ 2 A' andB 2 A' PESs of HeH + 2 in the reactant channel, obtained from the diabatic V ij surfaces, whereÃ 2 A' correlates with He + ( 2 S)+H 2 (X 1 + g ) for r < 2 a 0 and with He( 1 S)+H + 2 (A 2 + u , unbound) for r > 2 a 0 . To analyze the accuracy of the fit, we compare the analytical nonadiabatic coupling matrix elements (NACMEs) in the adiabatic representation, obtained from the fitted diabatic energies and coupling, with the ab initio results calculated using the MOLPRO program package (Werner et al., 2018). As done previously for H + 4 and H + 5 in Sanz-Sanz et al. (2015), the analytical NACMEs can be calculated from the generalized Hellmann-Feynman theorem, whereĤ el is the electronic Hamiltonian and the rhsm is obtained from the derivatives of the diabatic energies V ij (Sanz-Sanz et al., 2015). Ab initio calculations have been performed using the Multi-Reference Configuration-Interaction method, with the aug-cc-pVTZ basis set of Dunning (1989) and Woon and Dunning (1994). The ab initio NACMEs are obtained using a first order difference method with an interval of 0.01 a 0 , and that for R AB = r is compared with the fitted one in the bottom panel of Figure 2. The agreement, in form and in position, between both is excellent, indicating that the equation used for the interaction term V 12 is also appropriate to reproduce the NACMEs. The probability density of the first two vibrational states of H 2 (X 1 + g ) has been included in the top panel of Figure 2, in order to analyze the reason why the reaction is faster when H 2 is vibrationally excited (Aguillon, 1994(Aguillon, , 1998. While the maximum of the probability density of the vibrational state v 0 = 0 is found at 1.4 a 0 , those in the first excited vibrational state v 0 = 1 are at 1.25 and 1.73 a 0 . The latter is close to the region in which the NACME has is maximum, that explains the enhancement of the reactivity as we shall see in section Reaction Probabilities. In Figure 3 the ratio between the interaction term and the energy difference of the diabatic states, V 12 /(V 22 -V 11 ), is plotted as a function of the R HeH and R HH ′ distances, for several values of the angle θ ≤ 180 • between these distances. This ratio is important for the calculation of the adiabatic PESs from the fitted diabatic ones, according to: As expected, V 12 /(V 22 -V 11 ) changes sign and its absolute value is maximum in the region of the diabatic crossing line shown in Figures 2, 3 of Aguado et al. (1993).
In conclusion, these results show that the collisional dynamics of DCT He + ( 2 S)+H 2 (X 1 + g )→ He( 1 S)+H( 2 S)+H + can be investigated with high accuracy in the present diabatic electronic representation, using the V 11 , V 22 , and V 12 PESs.

Collision Dynamics
Since many years we are presenting our quantum theory (Petrongolo, 1988) and results of non-adiabatic effects in spectroscopy of triatomics and dynamics of atom+diatom collisions and we here report a brief summary, following our work on the CI dynamics of the OH(A 2 + )+H( 2 S) reaction (Gamallo et al., 2013).
The He + +H 2 collision is described by reactant Jacobi coordinates R, r, and γ , by a body-fixed reference frame with the z axis along R, and by a HeH + 2 spinless rovibronic Hamiltonian H, which contains the electronicĤ el and the total angular momentumĴ operators.Ĥ is represented in an orthonormal basis of diabatic electronic states (1) 2 A 1 and (2) 2 B 2 , radial grid |Rr>, associated Legendre |jK>, and symmetry Wigner states |K+p>. Here (1) 2 A 1 and (2) 2 B 2 are coupled byĤ el owing to the CI,h K is theĴ z eigenvalue, and we omit J and its space-fixed Z component in |K+p>, where the total parity is p = (−) J+K min with K min = 0 or 1. The 2J+1 values of K are thus factorized in two noninteracting groups, with K min ≤ K≤ J, of dimensions J+1 or J according to K min = 0 or 1, respectively.
Initial-state-resolved reaction probabilities are computed through the quantum, real WP formalism of Gray and Balint-Kurti (1998) and Meijer et al. (1998), essentially equal to the Chebyshev approach of Guo (2012). Shortly, an arccos mapping of the HeH + 2 time-dependent Schrödinger equation is solved recursively, using a scaled and shifted HamiltonianĤ s and starting from an initial and complex WP |ψ 0 >= |a 0 > +i|b 0 > . This initial WP describes the entrance channel He + ( 2 S)+H 2 (X 1 + g ), with the diabatic electronic state (1) 2 A 1 and the R-dependent term where µ R is the He + +H 2 reduced mass. The r and angular terms of |ψ 0 > are the vibrational |v 0 > and rotational |j 0 K 0 > states of H 2 (X 1 + g ), and finally |K 0 +p> is the initial overall rotational species. The recursions are |a n+2 = 2Ĥ s |a n+1 − |a n , other real propagations, where the square root in Equation (5) is evaluated with a Chebyshev expansion, and Equation (6) is a standard Chebyshev propagation of just a real WP, which is also absorbed at R > R abs and r > r abs by the Gaussians exp[-C R abs (R-R abs ) 2 ] and exp[-C r abs (r-r abs ) 2 ], respectively. At the end of the propagation, we obtain the probability via a time-to-energy Fourier transform and a flux analysis (Meijer et al., 1998) on the (2) 2 B 2 PES. We compute non-adiabatic initial-state-resolved reaction probabilities P J v 0 (E coll ), with the initial WP on the HeH + 2 (1) 2 A 1 diabatic PES and H 2 (X 1 + g ) in the ground and first excited vibrational state, v 0 = 0 and 1, and in the ground rotational state j 0 = 0. Then the initial K 0 is equal to K min = 0 and the total parity p is (-) J . Here we have where (1) and (2) are the diabatic electronic states, v ′ , j ′ , and K ′ label open vibrational, rotational, and helicity states of the products, respectively, and S J is the state-to-state parity adapted S-matrix at J. The total ICS is then defined by and the initial-state-resolved rate constant is where T is the temperature and k B is the Boltzmann constant. The calculations are done up to J = 150, using the coupledchannel formalism with Coriolis couplings among the K values, up to K max , that are necessary for converging the probabilities. Considering the H-H permutation symmetry, the numerical parameters of the WP propagations are listed in Table 1: they correspond to 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 3 eV, with E coll = 0.0001 or 0.001 eV below or above 0.2 eV, respectively. The most important parameters that affect the convergence of the calculations are the number of the propagation kilo-steps, kstep, the maximum K value, K max , and the dimension of the R grid, nR. We shall see in the next sections some convergence results with respect to kstep and K max . The probabilities obtained with nR = 461 of Table 1 are practically indistinguishable by those corresponding to nR = 559. All calculations are done with our J-K-parallelized Open MPI codes.

Reaction Probabilities
We plot in Figure 4 the opacities functions (2J+1)P J v 0 (E coll ) at J = 20 and 40, v 0 = 0 and 1, E coll ≤ 3 eV, and K max = 0, 1, and 2, with kstep = 5, which converges the probabilities within 0.1% above ∼0.5 eV. The curves at v 0 = 0 and 1 present a similar behavior, namely narrow resonance features at J = 20 and E coll ≤ 0.5 eV and a smooth, fast increase above. From the comparison among lower and upper panels we see that v 0 = 1 is more reactive  than v 0 = 0 by ∼one order of magnitude. Of course, this finding reflects the HeH + 2 (Ã 2 A ′ ) early potential barrier of 1.34 eV and the v 0 = 1 density maximum at 1.73 a 0 (Figure 2), close to the NACME maximum. Comparing the curves inside each panel, we observe a very fast convergence in K max , so that the Coupled-State approximation (McGuire and Kouri, 1974;Pack, 1974) at K max = 0 gives already reasonable results. Just two Ks converge the plots within the graphical accuracy, implying that the Coriolis couplings between K and K ± 1 vanish quickly when the initial K 0 is equal to zero. To reach the high accuracy claimed above at all the partial waves required to give convergent ICSs in all the collision energy range considered, we employ K max = 2 and 3 for v 0 = 0 and 1, respectively. In the J = 20 left panels, the dashed lines show the convergent results below 0.5 eV and the weight of the resonances with respect to the background. The resonance features are probably due to rotational metastable states of the C 2v ion-induced dipole minimum in the entrance channel. These states, embedded in the collisional continuum, are particular Feshbach resonances induced by the CI (Cederbaum and Friedman, 2003), with lifetimes of the order of the nanosecond. The v 0 = 0 opacities then increase in a monotonic way above 0.6 eV up to a maximum value of ∼0.06 at 3 eV, that is onetwo orders of magnitude larger than the strongest resonance. The v 0 = 1 opacities show the same behavior, but the value at 3 eV is just two times larger than the resonance maximum. However, to converge the resonance energy pattern very different values of the convergent parameters K max and kstep are required.
We show in Figure 5 the opacities at J = 10 and 20, for v 0 = 0, E coll ≤ 0.3 V, K max = 0, 3, and 4, and kstep = 120. In the upper panels a blow up of 0.05 eV also shows the details of the resonance features. From this figure we can observe that the convergence in K max is slower than for high energies. The curves at K max = 0 exhibit less peaks and the resonances shift in energy by increasing this parameter. Three main groups of resonances could be identified, decreasing and then disappearing at high E coll . This result suggests that the different peaks inside each group correspond to different bending energies of the collision complex (Aquilanti et al., 2005) while different groups correspond to different energies of the symmetric stretching motion. The v 0 = 1 curves have a similar behavior, with the K max convergence slightly faster and the resonances at lower energies.
In Figure 6 the resonance energy patterns at J = 5, 10, 15, and 20, v 0 = 1, and K max = 3 are plotted for E coll ≤ 0.11 eV at different kstep values. The K max value gives a graphical convergence to the curves and the logarithmic scale of the energy points out the resonance details. Here kstep is much greater than for the background and changes markedly with J and E coll . In fact, kstep = 350 is still not enough to converge all the J = 5 and 10 features, especially below 0.01 eV, but kstep = 160 perfectly converges the resonances at J = 15 and 20. As expected, the slowest convergent resonances are the narrowest ones, with largest lifetimes (Aquilanti et al., 2004). Moreover, the lifetimes of the collision complexes decrease at high J and E coll , because the centrifugal barrier increases and the resonances are less trapped inside the shallower well. Above J = 25 the well is so shallow that it does not support any resonance and the features disappear. Of course, increasing kstep means increase the collision time that must be of the same order of magnitude of the lifetime of the resonance intermediate to give convergent results.

Integral Cross Sections and Rate Constants
We plot in Figure 7 the total ICS at v 0 = 0 and 1, with J ≤ 150 to obtain convergent results up to 3 eV, noting that the convergent requirements change markedly with E coll and J. To minimize the computational effort, the total number of partial waves was therefore shared in different groups, and different values of K max and kstep were employed in each group, as shown in Table 2.
Notwithstanding the large kstep used for the lowest partial waves, test calculations show that some narrow ICS peak slightly increases with more steps. With these input data, a convergence within 0.1% is reached only by the resonance peaks above     (Aguillon, 1998).
v 0 E coll /eV present work (Aguillon, 1998) in two main groups, the first and stronger below 0.1 eV, and the latter and weaker from 0.1 to 0.3 eV. At v 0 = 0, the strongest resonance is equal to 0.003 Å 2 at 0.0116 eV and the maximum at 3 eV has nearly the same value. At v 0 = 1, the strongest resonance is larger than the value at 3 eV by ∼one order of magnitude, and is equal to 0.455 Å 2 at 0.0096 eV. In conclusion, the one-quantum vibrational excitation of H 2 (X 1 + g ) increases the cross section by more than two orders of magnitude at ∼0.01-0.02 eV and by ∼20 times at 3 eV.
As said in the Introduction, only the semiclassical theoretical works by Aguillon (1994Aguillon ( , 1998 have obtained ICSs at or below 3 eV. We present in Table 3 a comparison between our quantum ICSs and those estimated from Figure 4 of Aguillon (1998), showing the good agreement between the results that differ at the most by ∼24% at v 0 = 0 and E coll = 2 eV. Taking into account the different theoretical treatment and the small reactivity at these conditions, this implies that the Aguillon semiclassical treatment of the DCT reaction works remarkably well at this collision energies, where the ICSs do not present any quantum resonance.
Finally, we report in Table 4 and Figure 8 the initial-stateresolved rate constants k v 0 (T), at v 0 = 0 and 1, as functions of the temperature T up to 2,000 K. Their accuracy is ∼2%, considering the slow convergence of the resonances, and the rates are stable with respect to changes of the parameters in Tables 1, 2. Both rates increase quickly by a factor of ∼30 from 20 to 200 K, with a maximum at ∼250 K equal to 1.75 10 −14 and 1.97 10 −12 cm 3 sec −1 for v 0 = 0 and 1, respectively. This behavior is associated with the sharp ICS resonances above ∼0.006 eV and is similar to that of H+HeH + (De Fazio, 2014). Note that a rateconstant maximum was also found in the adiabatic ground PES dynamics of H+HeH + (Esposito et al., 2015) using the PES of FIGURE 8 | Initial-state-resolved rate constants at v 0 = 0 (red) and 1 (blue). Ramachandran et al. (2009), but it was at ∼10,000 K and due to a complete different mechanism. Then the rates slowly decrease with T and become nearly constant above 1,500 K, where they differ by two orders of magnitude. On the overall, this behavior is due to the ICS resonances shown in Figure 7 overlapped to a nearly Langevin decrease of the associated cross sections up to ∼0.5 eV.
Because we have considered just the ground rotational state of H 2 (X 1 + g ), j 0 = 0, the present rate constant at v 0 = 0 and 300 K underestimates by a factor of ∼6 the thermal experimental k(330) = (1.1 ± 0.1) × 10 −13 cm 3 s −1 (Johnsen et al., 1980), which is however 39 years old. The agreement increases if we consider H 2 in the excited rotational states that are open at 300 K, because test calculations suggest that the reactivity can be enhanced by ∼33% when j 0 = 1. Taking into account the H 2 (X 1 + g ) Fermi-Dirac nuclear spin statistics, the room-temperature populations of j 0 = 0 and 1 are ∼0.13 and 0.66, respectively, and this implies a rotational enhancement of the rate constant by a factor of ∼7. Owing to the room-temperature branching ratio of paraand ortho-H 2 and the rotational effects on the reactivity, we can roughly estimate than ortho-H 2 is ∼4 times more reactive than the para species. In closing this section, we have also done some test calculations of the reaction probabilities in the adiabatic approximation, on theÃ 2 A' PES of HeH + 2 (see Figure 1). Without reaching the accuracy and the stability of the results in section Reaction Probabilities, we present an example at v 0 = 0, J = 0, and kstep = 120, contrasting in Figure 9 adiabatic A 2 A' and non-adiabatic CI (1) 2 A 1 /(2) 2 B 2 results. Although both probabilities present a resonance structure at low E coll and increase above ∼1 eV, the reactivity is dramatically different, with theÃ 2 A ′ probability hugely larger than the CI one, from 6 orders of magnitude at 0.001 eV to 2 orders at 3 eV. This finding shows that theÃ 2 A ′ adiabatic PES drives the WP into the exothermic product channel, following a C s pathway that avoids the C 2v barrier of the CI. On the contrary, the non-adiabatic CI WP remains essentially on the initial and repulsive (1) 2 A 1 diabatic PES V 11 , owing to the weak non-adiabatic interaction. As expected, only at energies larger than 3 eV the two probabilities seem to merge and the adiabatic approximation is probably less worse. We roughly estimate that the adiabaticÃ 2 A ′ rate constants at 330 K is ∼1,000 times larger than the experimental value (Johnsen et al., 1980).

CONCLUSIONS
In this article we have presented quantum non-adiabatic DCT results of the He + +H 2 collision, employing an electronic diabatic representation, previously computed by Aguado et al. (1993), a WP time dependent formalism, and a flux analysis. Specifically, we have taken into account the non-adiabatic CI coupling between the first two excited diabatic PESs of HeH + 2 . Dynamical calculations are performed for the ground and first excited vibrational states of H 2 , for investigating vibrational effects on the DCT dynamics, and collision energies up to 3 eV are considered. Reaction probabilities and ICSs present strong and narrow resonance features up to 0.5 eV, due to quasi-bound molecular states embedded in the continuum and trapped in the ion-induced dipole minimum of the reactant channel, near the CI. These features are very hard to converge and affect all the DCT dynamics. These intense resonance features have a determinant role in the radiative chargetransfer formation (Mrugala et al., 2013) of stable HeH + 2 , probably present in the interstellar medium (Tennyson and Miller, 1987). At higher collision energies the computational load reduces drastically owing to a near conservation of the helicity quantum number. Our results are in reasonable agreement with previous experimental and theoretical studies of this reaction, confirming the strong vibrational enhancement previously founded. In fact, rate constants increase of about two orders of magnitude just adding one vibrational quantum to the H 2 reactants.
At the best of our knowledge, these are the first rigorous quantum DCT calculations presented for a chemical system, made possible by the joint implementation of time-dependent WP, time-to-energy Fourier transform, and flux methods. On the other hand, a more rigorous approach, as time independent close coupling calculations, cannot be attempted at the present state-of-the-art of the reaction dynamics theories, because of the difficulty of these methods to treat the three-bodies breakup (see e.g., (Pack et al., 1998), and references therein). The results achieved could have relevant consequences in astrophysics and in particular in the early Universe evolution models (Bovino et al., 2011). Although the reaction is probably negligible when the molecular hydrogen is relaxed in its ground vibrational state, the strong increase with the vibrational energy suggests that it could have a role during the early stages of the adiabatic expansion at high redshift, destroying the molecular hydrogen formed and slowing down further cooling of the primordial gas. Of course, to determinate its role, evolution Universe models with at least a vibrational resolution of the chemical network (Coppola et al., 2011(Coppola et al., , 2012 are required. Moreover, this reaction could be also relevant in other different astrophysical environments, as for example in the study of the Sun atmosphere where the temperature is very high and hydrogen and helium are very abundant species.

AUTHOR CONTRIBUTIONS
All authors listed have made a substantial, direct and intellectual contribution to the work, and approved it for publication.

ACKNOWLEDGMENTS
We are very grateful to P. Gamallo and to the Dipartimento di Biotecnologie, Chimica, e Farmacia, Universita' di Siena, Italy, for the availability of essential computer resources. We want also to thank O. Roncero for fruitful discussions. The research leading to these result has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007(FP/ -2013/ERC Grant Agreement n. 610256 (NANOCOSMOS), the HPC Europa3 Project (INFRAIA-2013-1-730897) with the support of the EC Research Innovation Action under H2020 Programme and the COST Action CM1401 (Our Astrochemical History). Also, we acknowledge support by the Ministerio de Economía y Competitividad under grants No.