Vibrational, non-adiabatic and isotopic effects in the dynamics of the H2 + H2+ → H3+ + H reaction: application to plasma modelling

The title reaction is studied using a quasi-classical trajectory method for collision energies between 0.1 meV and 10 eV, considering the vibrational excitation of $ {\rm H}_2^+ $ H2+ reactant. A new potential energy surface is developed based on a Neural Network many body correction of a triatomics-in-molecules potential, which significantly improves the accuracy of the potential up to energies of 17 eV, higher than in other previous fits. The effect of the fit accuracy and the non-adiabatic transitions on the dynamics are analysed in detail. The reaction cross section for collision energies above 1 eV increases significantly with the increasing of the vibrational excitation of $ {\rm H}_2^+ $ H2+( $ v^{\prime} $ v′), for values up to $ v^{\prime} $ v′=6. The total reaction cross section (including the double fragmentation channel) obtained for $ v^{\prime} $ v′=6 matches the new experimental results obtained by Savic, Schlemmer and Gerlich [Chem. Phys. Chem. 21 (13), 1429.1435 (2020). doi:10.1002/cphc.v21.13]. The differences among several experimental setups, for collision energies above 1 eV, showing cross sections scattered/dispersed over a rather wide interval, can be explained by the differences in the vibrational excitations obtained in the formation of $ {\rm H}_2^+ $ H2+ reactants. On the contrary, for collision energies below 1 eV, the cross section is determined by the long range behaviour of the potential and do not depend strongly on the vibrational state of $ {\rm H}_2^+ $ H2+. In addition in this study, the calculated reaction cross sections are used in a plasma model and compared with previous results. We conclude that the efficiency of the formation of $ {\rm H}_3^+ $ H3+ in the plasma is affected by the potential energy surface used. GRAPHICAL ABSTRACT


Introduction
Hydrogen, as the most abundant element in the Universe, plays a fundamental role in star formation and the chemical evolution of molecular Universe.Its molecular forms previous theoretical simulations [50,56,61].Moreover, for collision energies above 1 eV, the new experimental measurements show differences with previous experimental ones [46-49, 53, 54], all of them scattered in a rather wide interval of cross sections.
There are two main goals in this work.First, we focus on the theoretical simulations of the H + 3 formation reaction, Eq. ( 1), to reproduce the new experimental data above 1 eV [1].Second, we study hydrogen or deuterium plasma models in order to analyze the effects of the calculated cross sections and rates on the H + 3 or D + 3 densities.This work is organized as follows.In section 2, we develop new potential energy surfaces (PESs), which include non-adiabatic effects and increase the accuracy of the fit.This new fit uses a Neural Network (NN) method to describe the four-body term, to improve a zero-order description using a Triatomics-in-Molecule treatment, which accurately fits very precise ab initio calculations, over a broader energy interval than previous fits, up to 17 eV.In section 3, we study the reaction dynamics using a quasiclassical trajectory (QCT) method, including transitions among different electronic states, using the fewest switches method of Tully [66].The reaction dynamics is studied for collision energies between 1 meV and 10 eV, and for several vibrational states of H + 2 reactant, as well as for the deuterated reaction, focusing on high energy reactive cross sections recently measured by Savic et al. [1].In section 4, we investigate how the calculated reactive cross sections affect the population density of H + 3 and D + 3 in a CR models of H 2 /D 2 plasma.Finally, in section 5 some conclusions are extracted.

Potential energy surfaces of H + 4
In this work several PESs are used, and are listed here to clarify the differences: • PES1: This PES developed by Sanz-Sanz et al. [67], is the most accurate developed so far for this system.It is built as a sum of a triatomics-in-molecules (TRIM) term, H T RIM , plus a four body term, H M B .The TRIM term is a generalization of the Diatomics-in-Molecules (DIM) [68][69][70][71] method, in which the electronic Hamiltonian is factorized as a sum of triatomic and diatomic fragments as [67,72] Ĥi e = n>i,o>n where Ĥ+ ip are the monoelectronic Hamiltonians of H + 2 fragments and Ĥ+ ino (n − i, o − i) are the bielectronic Hamiltonians (for n − i, o − i electrons) describing the H + 3 system for the ino nuclei.The TRIM representation consists of a 8×8 matrix, whose elements have H + 3 and H + 2 matrix elements of different electronic configurations, in which each hydrogen atom is described by a 1s function, except one corresponding to H + .In PES1 [67], the triatomic terms included in the TRIM matrix were built as the 3×3 DIM matrix for each H + 3 fragment plus a three-body term added to describe the ground state of either singlet or triplet symmetry.The non-adiabatic terms in these triatomic fragments are therefore approximated by those obtained in a DIM treatment.
PES1 is built adding a four-body correction (MB term) to this TRIM description.This MB term is expressed as a linear combination of four-body polynomials that are symmetric with respect to the permutations of identical nuclei (Permutationally invariant polynomials or PIP).These PIP are built in terms of Rydberg polynomials of the interatomic distances (p i = r i exp(−α i r i )) [73][74][75].The set of linear and non linear coefficients, α i , are optimized to minimize the difference with the ab initio energies obtained with the multi reference configuration interaction (MRCI) method using the aug-cc-pV5Z basis set [76].In PES1, the same four-body correction term is added to all the diagonal elements of the TRIM matrix.
The reactive cross sections calculated on this PES1 shows an excellent agreement with the experimental results for collision energies between 1 meV and 1 eV [1,[61][62][63].
• PESTRIM8×8 and PESTRIM1×1: This PES is based on the TRIM model explained above, and is an improvement made in this work.The main difference is that the triatomic H + 3 is represented by 3×3 diabatic matrices fitting its three lowest singlet and triplet electronic states [77] to MRCI ab initio energies extrapolated to the complete basis set (CBS) limit [78].This improvement is crucial to study the non-adiabatic dynamics of H + 4 .Analytical derivatives of the potential energy surfaces and non-adiabatic couplings are calculated based on the Hellmann-Feynmann theorem as described in [61].The full PESTRIM8×8 potential correspond to the eight adiabatic eigenvalues including the non-adiabatic coupling terms.PESTRIM1×1 only considers the ground adiabatic energy.
where H + 4 is formed, while showing excellent agreement for H + 3 + H or H 2 + H + 2 rearrangement channels.For this reason a four-body neural network term is added to PESTRIM1×1 to improve the accuracy of the system, as described below.

Many Body Neural Network term
Many body potential energy terms are widely used to represent potential energy surface of chemical systems, either in a many body expansion [74,79,80], where one up to N -body terms are summed, or as correction terms of a zero-order description of the potential, described by the TRIM method [67,72] or by a reactive force field matrix [81][82][83].
In this work a many body neural network (MB-NN) [83] potential energy term has been built as a PIP-NN [84], in which the neural network is fed with a permutational invariant polynomial representation of the molecular geometry.A PIP is constructed by projecting a polynomial of the interatomic distances into the totally symmetric irreducible representation of the desired permutation group.A generator of these PIPs is defined as: where N d is the number of interatomic distances, p i is a function of the interatomic distance r i , l n i is the exponent of the ith monomial for the nth polynomial and Ŝ is the projector to the totally symmetric irreducible representation of the permutation group.A common choice of p(r) in PIP-NN-PES is the decaying exponential p i (r) = exp(−α i r).
The set of PIP has to be carefully filtered so that it is purely composed of N -body functions.This means that any of these functions evaluated on a geometry where the N bodies are not closely interacting should be zero.In case these terms included any lower body functions, we have shown that it would add spurious interactions between fragments [85], which specially affects the long range regions.
In appendix A it is shown that a polynomial in Eq. ( 5) can be represented as a graph, and that the subset of N -body polynomials corresponds to those whose respective graph is connected, meaning that there exists a path which connects all the vertices (particles) of the graph.In this way we can guarantee that any N -body polynomial, with p(r) defined as a decaying exponential or Rydberg function, will tend to zero as any particle or set of particles moves away from the rest.This will automatically make zero a PIP-PES and provide constant descriptor for a NN-PIP-PES, that returns a constant energy out of the N -body region, which can be trained to be as close to zero as possible and that produces no net force, since its derivative with respect to the atom coordinates is zero.
The four-body term developed here is expressed as a feed forward neural net with 11 input neurons and two hidden layers with N 2 = 32 and N 3 = 77, and sigmoid non linearities σ = 1/(1 + exp(−x)): were w (l) matrix (with elements w (l) ki ) and b (l) vector (with elements b (l) i ) represent the trainable weights and bias on layer l.The 11 input neurons correspond to four-body PIPs produced by setting a maximum polynomial degree max( l i ) = 5 and maximum monomial degree max(l i )=2.All lower body polynomials that would introduce spurious energy contributions in reactants and products channels are filtered following the steps detailed in appendix A.
The four-body term is trained with NeuralPES, an in-house Python code based on PyTorch [86].New ab initio points have been used in the fit of this work, of higher accuracy of those used in PES1 [67].The energies are obtained using a two point extrapolation method to complete basis set (CBS) [78], using the results obtained with the aug-cc-pV5Z and aug-cc-pV6Z basis sets.Around 33000 ab initio points were calculated, including the geometries of Ref. [67] and new ones selected to increase the accuracy of the PES above 2 eV.These last new points were chosen from QCT trajectories at different collision energies (from 1 eV to 10 eV) taken on the TRIMPES1×1 to populate physically accessible configurations at this high energy regions.The complete set of ab initio points is randomly split into a training set (containing 80% of the data) and validation and test sets (with 10% of the points each).The training set consists on 27294 geometries, with energies up to 17 eV over the H 2 + H + 2 asymptote, mostly corresponding to four-body interactions, and elongations to lower body geometries.The training process aims to minimize the root mean squared error between the ab initio and PESTRIM1×1 + H M B energies.

Analysis of the different PESs
As all the four-body terms described above vanish as the systems tend to reactants or product asymptotes, the long-range interactions are purely described by the triatomic terms of the TRIM model.In all the triatomic fits considered here, the long range terms + H reaction at low temperatures.The measured cross section was then scaled to reproduce the cross section calculated with the PES1 potential [61] at a single collision energy.The excellent agreement between calculated and scaled experimental cross sections between 0.5 and 5 meV demonstrates the good behaviour of the long-range interactions included in PES1.The new PESs introduced in this work, describes the long-range interaction even more accurately, by using improved long-range interactions in the triatomic fragments [77].
The RMS errors for the different potential energy surfaces are presented in Table 1, in different energy intervals.The improvement of PES-NN over PESTRIM1×1 is clear in all energy ranges due to the enhancement of the H + 4 channels as presented in Table 2. PES1 and PES-NN show comparable RMSE for energies below 2 eV, but when ab initio points above 2 eV are added, PES-NN shows a much higher accuracy.
The whole configuration space is divided as follows: when all the interatomic distances are below R thres = 4 Å, the system is taken to be in H + 4 region; if all interatomic distances of a triad of hydrogens are shorter than R thres , the region is taken to be H + 3 + H; if any two pairs of atoms present interatomic distances shorter than R thres , the region is denoted as H 2 + H + 2 ; if only one internuclear distance is < R thres , the region is called H 2 + H + H + ; otherwise, if all the interatomic distances are large, the system is in fully dissociated.Following this division, the top left panel of Figure 1 shows the distribution of the data set among the different regions as defined above, as a function of the energy, taking the H 2 + H + 2 reactants asymptote as the zero of energy.Most of the data correspond to H + 4 and geometries which connect this channel to H 2 + H + 2 and H + 3 + H.As can be seen in the top right and bottom panels, the PESTRIM1×1 is highly accurate in the latter two channels, but shows deviations in the H + 4 region, up to several eV.The effect of the four-body term on PES-NN is decisive for the proper description of H + 4 region.PES1 was fitted for energies up to ≈ 2 eV, and it presents large energy deviations over this threshold.On the contrary, PES-NN keeps a high precision at high energies, which are the interest of the present work.
In Figure 2, H 2 + H + 2 approaches for different θ 1 and θ 2 angles in their equilibrium geometry are shown for the three potential energy surfaces and compared with ab initio calculations.PESTRIM1×1 tends to predict a larger energy for H + 4 geometries, while PES1 and PES-NN yield effectively the same description.As the interfragment distance increases towards H 2 + H + 2 the four-body terms in both PES1 and PES-NN go to zero and all that remains are the respective TRIM terms.The improvement of the new PES-NN is better seen when the diatomic fragments are not at the equilibrium geometry, for energies above 2 eV.
The more accurate description of the higher energy regime, energies larger than 2 eV, can be seen in Figure 3 where H approaches to a compressed H + 3 .The differences between PES1 and PES-NN are more pronounced when bonds are compressed than stretched.This is

Quasi-classical trajectory and surface hoping calculations
The quasi-classical trajectory calculations are performed with the MDwQT code [61,82,89,90].When considering several coupled adiabatic electronic states, the fewest switches approach of Tully [66] is used as described in [61]   are sampled with the usual Monte Carlo method [91].The initial conditions for the vibrational modes of H 2 and H + 2 are quantized using the adiabatic switching method [92][93][94], yielding vibrational energies within 0.3 meV with respect to the exact vibrational levels of each diatomic fragment.The rotational states of H 2 and H + 2 are set to zero in these studies, and the initial distance between the two center-of-mass is set to 105 bohr.The initial impact parameter, b, is sampled between 0 and B, according to a quadratic distribution on b, where B is determined for each energy according to a capture model [95] as B = (α/2E) 1/4 , for a charge-induced-dipole interaction, described by −α/2R 4 , where α is a constant, with dimensions [EL 4 ], which is proportional to the average polarizability of H 2 , β, at its equilibrium configuration, as α = βe/4π 0 .All trajectories are stopped when any internuclear distance becomes longer than 125 bohr, where they are analysed.
The reactive cross section for each collision energy, E, is calculated as [91] σ where N r is the number of trajectories leading to products and N tot is the total number of trajectories with initial impact parameter lower than b max , the maximum impact parameter for which reaction takes place at energy E.Here we have considered H 2 (v=0), while H + 2 vibrational level v varies between 0 and 6.For each (v, v ) couple and each energy a set of 10 5 trajectories are run, with energy error lower than 0.01 meV.
The final energy distribution of H + 3 products is also analyzed.We simply evaluate classical energies, without trying to consider the permutation symmetry of identical fermions (for H + 4 ) or bosons (for D + 4 ).This is done in three steps.First, the kinetic energy of H + 3 and H products are calculated and substracted.Second, the rotational angular momentum of H + 3 products is evaluated, and its rotational energy.By setting the origin of energy at the bottom of the H + 3 well, the remaining energy corresponds to vibrational energy.Here, we do not attempt to assign the internal vibrational modes, which deserves further development and is led for a future work.

H 2 + H + 2 (v ) collisions in PES1
The new experimental results for the title reaction of Savic et al. [1] differ from previous ones, both experimental and theoretical, above 2 eV.The difference with previous experimental data, which are scattered above 1-2 eV, may be due to different conditions in the generation of the reactants.In low temperature plasma the vibrational temperature of H 2 is of the order of 2500 K, so that the population of H 2 (v=1) is expected to be lower than 10% [31][32][33].On the contrary, H + 2 is formed by electronic impact or photoionization, which may yield to different vibrational and rotational excitations.Vibrationally excited H + 2 can partially thermalize, yielding to different initial conditions in different experimental setups.As an example, recent theoretical calculations [96,97] on the H + H + 2 charge transfer reaction, and some isotopic variants, have found that the reaction cross section highly depends on the initial vibrational Collision energy (eV)  [1], from Glenewinkel-Meyer and D. Gerlich [55] and Shao and Ng [53], and the theoretical results of Eaker and Schatz [56].
state of the diatomic ion.Following this idea, in this work we have performed QCT calculations of the cross section of the reaction for several initial vibrational states of H + 2 reactant using the global PES1 of Refs.[61,67], which are shown in Figure 4.
For exothermic reactions without a barrier, the long-range interaction between the reactants dominates the reaction dynamics.In the case of charge-induced dipole longrange interactions the cross section for exothermic reactions takes the form [95] σ This is approximately the behaviour of the cross sections for energies below 1 eV for every initial vibrational state v .Therefore, we can conclude that in this energy interval reaction dynamics is dominated by long range interactions, independently of the initial vibrational state of H + 2 (v ).However, above 1 eV the reactive cross sections present important differences among the v considered, showing that the vibrational excitation has a strong impact on the reactivity.In general, the reactive cross section increases with increasing v , which could be simply understood assuming that H + 2 can break more easily.In the left panels of Figure 5 the two main mechanisms to form H + 3 are separated as H-hop and proton-hop, corresponding to the fragmentation of H 2 or H + 2 reactants, respectively.For collision energies below 1 eV, the cross section for the proton-hop mechanism is slightly larger.However, above 3 eV the H-hop cross section is larger for v =0, and the two mechanisms tend to a rather similar value as v increases.This means that the vibrational excitation of H + 2 does not produce significant increase of the proton-hop mechanism.This can be explained looking at the right panels of Figure 5, where the H hop Figure 5. H-hop and proton-hop cross sections (left panels) and maximum impact parameter, bmax, (right panels) for the H 2 (v=0,j=0) + H + 2 (v , j =0) reaction obtained with QCT calculations for different vibrational states v of H + 2 .
maximum impact parameter, b max , is plotted for each proccess and inital vibrational state.Above 1 eV, b max increases from v = 0 to v =6, in nearly an identical quantity for the two mechanisms.We conclude that the increase of the cross section when varying v is due to the growing of the H + 2 subunit, whose right turning point increases from 1.25 to 2.12 Å, for v =0 and 6 respectively.Therefore, at energies above 1 eV, when long-range interactions are not able to produce important deviations among the reactants, the size of the two reactants approximately determines the value of the maximum impact parameter.At these higher energies a more complex reaction mechanism occurs, in which H 2 may nearly insert in the H + 2 bond, specially at high v excitations.

D 2 + D +
2 collisions in PES1 The reactive cross section for the D 2 (v, j) + D + 2 (v , j ) collisions is shown in Figure 6 for v = j = 0 and v = j = 0.For energies below 1 eV, the results for D 2 + D + 2 closely match those for H 2 + H + 2 .This can be explained in terms of the Langevin model, in which the cross section of Eq. ( 8) does not depend on the mass of the reactants.This explains why the reaction cross sections for D + 4 and H + 4 are nearly the same.However, it should be noted that for reactions involving partially deuterated species, the reaction cross section presents larger differences, as those already reported in the theoretical study of Ref. [61].The shift of the center-of-mass with respect to the geometric center of the two diatomic reagents introduces important differences, specially related to the effect of the rotation.In particular, in the homonuclear neutral H 2 or D 2 case for j=0, only the charge-electric dipole term affects the long-range interaction, determining the reactive cross section below 1 eV.In this case, the chargeelectric quadrupole term vanishes when integrating over the angular coordinate for the isotropic j=0 case (but not for j > 0).

Non-adiabatic effects and fit accuracy
The increasing of the H + 2 (v ) vibrational excitation yields to an increase in the reactive cross section.However, this increasing, even for v =6 does not match the new experimental data by Savic et al. [1].Since the cross section as a function of the v excitation seems to converge to the value of v =6, here after we shall focus on the two limiting cases, v = 0 and 6.
In order to investigate the role of electronic transitions among the lower electronic states, in this work we use the PESTRIM8×8, comparing it to the results on the PESTRIM1×1 model.The dynamical results obtained for these two potentials are shown in Figure 7, and compared with the data obtained with PES1.All the results present a similar behaviour, with small differences in the logarithmic scale used in the figure.All the results converge to the same value at the lowest collision energy considered here, 1 meV, since all the PESs used in this work have similar long range interactions.
At intermediate energies, between 0.01 and 2 eV, the PESTRIM8×8 and PESTRIM1×1 cross sections are slightly different, showing a small effect of nonadiabatic transitions in the dynamics.Curiously, these non-adiabatic effects seem to produce a decrease on the reactive cross section, in contrast to what is needed to match the values at higher energies of the experimental results of Savic et al. [1].Collision energy (eV) . QCT H 2 (v=0,j=0) + H + 2 (v =6,j =0) → H + 3 + H reactive cross sections obtained with the PESTRIM8×8 using surface hoping method (blue), the PESTRIM1×1 (red), the PES1 (black) and the new PES-NN (green) potential energy surfaces.
In figure7 we also compare with the cross section obtained with the adiabatic PES1, which nearly matches the results obtained with the adiabatic PESTRIM1×1 up to 3 eV.Above 3 eV, the results obtained with PES1 are higher than those for PESTRIM1×1, showing that the four body term may be important.However, for the new PES-NN, of higher accuracy than PES1, the reaction cross section is nearly identical to that of PESTRIM1×1 for collision energies below 1 eV and very close to those obtained with PESTRIM8×8 for energies above 4 eV.From this comparison we may conclude that the better accuracy of the four body term of the new PES-NN does not improve significantly the differences between the simulated and measured [1] reactive cross-sections.
In order to analyze whether there are other excited states not well described by the TRIM approximation, we have performed ab initio calculations of the four lower electronic states, considering a larger electronic basis set, with extra orbitals added to describe the 2s and 2p electronic states of atomic hydrogen.With this larger basis set, we have found that the energies obtained (extrapolated to the complete basis set) differ only a few tenths of meV with the previous ab initio calculations, and that no higher electronic state appears below 10 eV.The main mechanism giving rise to the DF channel consists of three steps, as it is shown in the left panels of Figure 9. First, a highly vibrationally excited H + 4 intermediate is formed by insertion of H 2 in the elongated H + 2 , which lives short time.Second, a first H atom is ejected (in the Figure is atom 2 with charge 0), forming a very ex-  16) Savic( 20)

cited (H +
3 ) * , which lives from 80 to 160 fs, approximately.In the third step, one of the atoms of H + 3 (atom 4 in the Figure ) dissociates, carrying the positive charge, thus leading to neutral H 2 .Since this third step occurs much later, it could explain that experimentally the metastable (H + 3 ) * would be detected together with more stable H + 3 products.It is also important to notice that the energy transfer among particles of identical mass is possibly overestimated in a classical treatment, as the one used in this work.If this is the case, it would support the inclussion of the DF cross section as a part of the total reaction cross section.Quantum calculations are needed to solve this problem, but they are rather challengig at the high energy considered.The right panels of Figure 9 show another DF mechanism: in this case H 2 and H + 2 are produced at 70-80 fs.H + 2 is vibrationally very excited and dissociates later, at ≈ 100-110 fs.In this case, there are two degenerate electronic states, each one corresponding to the charge in one of the ejected atoms, at long distances from the H 2 fragment.In the ground electronic state, shown in the lower panels, the charge is exchanged between the two identical atoms, because of this degeneracy, and shows the nature of the surface hopping occurring in the products channel when including several electronic states (PESTRIM8×8).The electronic transition occurs among degenerate electronic states describing H 2 + H + H + and H 2 + H + + H products, what explains the small effect of including electronic transitions in the reaction dynamics.
In addition, the charge transfer in the entrance channel occurs between H 2 and H + 2 , when the two reactants have the same internuclear distance (see bottom panels of Figure 9).In this situation there is a degeneracy between the two lower adiabatic states, as discussed in detail in Ref. [61].
The energy distributions of the H + 3 + H products are shown in Figure 10, and it is nearly quantitatively the same for all the PESs.The energy difference between the two vibrational states of H + 2 (v = 0 and 6) is approximately 1.4 eV, close to the exoergicity.Middle panels shows the energies of the four lower electronic states.Upper panels show some characteristic internuclear distances needed to characterize the trajectory.R ij refers to the distance between atom H i and H j For low collision energies, 0.1-1 meV, the initial vibrational energy of H + 2 ( v = 6) is 1.40 eV higher than that of H + 2 (v = 0), and the vibrational energy of the corresponding H + 3 products is also higher, but only by ≈ 0.9 eV.Therefore the remaining 0.4 eV are nearly equally distributed between the rotational and translational energies of the products.Rotational energy increases with collision energy, except for the v = 6 above 1 eV, where rotational excitation reaches a plateau and seems to start decreasing.The translational energy always increases for v = 0, while for v = 6 it slightly decreases below 0.1 eV, and then increases again.The vibrational energy of products shows two different behaviors: below ≈ 0.7 eV, the vibrational energy slightly decreases (v =0) or remains constant (v =6); above this energy the vibrational energy increases sharply.
Such behavior allows to assign two different reaction mechanisms.Below ≈ 0.7 eV, the impact parameter (in Figure 5) is large, supporting a stripping mechanism, in which the long range interactions attract the reactants to each other, originating orbits enhancing the relative angular momentum between the two reactants.Above ≈ 0.7 eV, however, the impact parameter is rather small, and the reaction occurs by an insertion mechanism, in which the H + 3 is greatly excited vibrationaly, specially as initial vibrational and translational energy of the reactants increases.
The dissociation energy of H + 3 is 4.34 eV, very close to the value reported by [99] of 4.35 eV, and the average vibrational energy distribution of H + 3 reaches values in the 4-5 eV interval, i.e. values above the dissociation energy explaining why H + 3 dissociates, leading to the DF channel.3 + H products, as a function of the collision energy for the H 2 (v = 0) + H + 2 (v ) collisions, for v =0 (left panel) and v =6 (right panel).The origin of energy is in the bottom of the well of the H + 3 products.The initial vibrational energies of the reactants is 0.413 and 1.673 eV for the v = 0 and 6, respectively, with respect to the minimum of each fragment.The potential energy difference between reactants and products is 1.816 eV, and when ZPE are accounted for, the exoergicity of the present PES becomes 1.688 eV.

Plasma modelling
To model the hydrogen plasma, we shall consider a gas-discharge vessel of cylindrical symmetry, so that only z, parallel to the central axis, and R, the distance from the center to the walls, will be considered, with the cylinder of infinite length in this case.The gas is initially in the form of neutral H 2 , and after the discharge ignition new species are formed, H, H + , H + 2 , H + 3 and electrons.H − is neglected under the conditions considered here.To model the abundance of these species in the stationary condition we use a model similar to those already described previously [42,100], which is outlined in the appendix B, including all the processes listed in Table B1, in which the vibrations of molecular species, H 2 , H + 2 and H + 3 are not considered.The plasma model is done for pure hydrogen and pure deuterium gases.For modelling deuterium plasma, the cross sections for electron collisions and radiative transitions of D species are used by those of H species.
The cross section for H 2 + H + 2 and D 2 + D + 2 reactive collisions calculated in this work are included in the models presented below.The plasma modellings are also done using the cross section by Janev et al. (2003) [98], which are compared with the present calculations in Figure 6.
Here, the electron temperature T e = 5.55 eV and electron density n e = 2.07 × 10 12 cm −3 are used, which were measured by a Langmuir probe for D plasma in a long cylindrical vessel [40].The molecular temperature is set to T m = 0.026 eV (300 K) and the atomic temperature is T a = 0.052 eV (600 K).At these realistic temperatures, the relevant energies are below E col = 1 eV, so that the H + 2 vibrational excitation has no significant effect.
Figure 11 shows the resulting population densities of D and H species depending on the cross section used.We can note, the D + 3 and H + 3 population densities vary significantly with the cross section used, while other species are nearly unchanged, as  Table 3. Rate coefficients for D 2 (v=0,j=0) + D + 2 (v =0,j =0) reaction and the modeled density of D + 3 depending on the cross sections shown in Figure 6.a.The values in the [ ] represent for the densities of H + 3 and the rate coefficient for H 2 (v=0,j=0) + H + 2 (v =0,j =0) reaction.

Conclusion
In this work a detailed study on the H 2 + H + 2 (v ) → H + 3 + H reactive cross section has been done using a quasi-classical treatment, for collision energies from 1 meV up to 10 eV and for several vibrational states of the H + 2 reactants and several isotopic variations.To this aim, new potential energy surfaces have been developed, one to include nonadiabatic transitions (PESTRIM8×8) and another to increase the accuracy in the whole energy range up to 17 eV (PES-NN).In all cases, it is found that from 1 meV to ≈ 0.5-1 eV the cross section behaves according to a Langevin law for charge induceddipole long range interactions.For energies above 1 eV, the simulated cross section decreases fast, below the Langevin limit, for all initial vibrational states of H + 2 (v ).However, for E > 1 eV, the reactive cross section exhibits a considerable increase with increasing v , and the results seems to converge at v =6.
It is found that the reactive cross section for v =6, summing the H + 3 + H and H 2 + H + H + channels, match very well the recent experimental measurements by Savic et al. [1], and also previous measurements [55,63], describing a broad energy interval from 0.5 meV to 10 eV.Moreover, the fact that for collision energies above 1 eV the measured reactive cross section show different values in different experiments, can be explained by a different vibrational excitation of H + 2 achieved in each experimental setup.
Experimentally, the H + 3 products are measured [1].The fact that in the QCT simulations, the cross section for the double fragmentation channel, H 2 + H + H + , needs to be considered to get an agreement with the new experimental results can be explained by a classical artifact, since QCT method usually overstimate the energy transfer, specially dealing with systems with equal masses.
The new reactive cross sections obtained for H 2 + H + 2 (v ) and D 2 + D + 2 have been included in a plasma model together with the widely used cross section [98].The resulting population densities of H + 3 and D + 3 are proportional to the rate coefficient, which in turn indicates the sensitivity of the population density to the adopted cross section.The new cross sections for vibrational sates (v ) of H + 2 will be useful for the state-resolved CR modeling in plasma with higher molecular temperature, where higher collision energies (above 0.7 eV) become significant and the differences in the reactivity of the vibrationally excited H + 2 will become important.
a general set of PIP.A graph is a collection of vertices and edges (G = (V, E)) [101].In an undirected graph the edges are non-ordered pairs of vertices.An undirected graph is said to be connected if every pair of vertices are joined by a path.
A relation can be established between the polynomials P n in Eq. ( 5) and an undirected graph, where the vertices represent the particles and the edges the monomials between them.
For instance, the following polynomial can be expressed as the following graph: This corresponds to a disconnected graph since there are various pairs of vertices which are not reachable, for instance from vertex 2 to 3.An example of a polynomial whose graph is connected is: A polynomial P is said to be a N -body polynomial if its corresponding undirected graph is connected.Any other polynomial that arises as ŜP with Ŝ being any operation of a permutation group will be a N -body polynomial if P is.Note that the effect of a permutation operation in the graph only affects the order of the vertices and relabel of the edges: Given a graph, there are simple algorithms as Depth-First Search (DFS) [102] to compute whether it is a connected graph or not, which recursively traverses the graph marking each visited vertice.If at the end of the execution all nodes were visited, the graph is connected.There exist a finite number of connected N -vertices undirected graphs which can be evaluated using the recursive formula [103]: These numbers are tabulated for N up to 16 in [104].One should note that the number of connected undirected graphs increases fast, for instance for a set of six vertices there exist 26704 of those, so in practice an upper limit on the number of edges has to be set.
At this point we have the tools to determine the number of N -body polynomials, as well as, given a polynomial, to determine if it is N -body.If the system presents some kind of permutational symmetry, a minimal set of N -body PIP will be generated by projecting the above polynomials onto the totally symmetric irreducible representation of the permutation group.The dimension of the minimal PIP set is necessarily lower or equal to the minimal one of polynomials, as many of them are related by permutation operations.Note that we never mentioned the exponents l of the monomials, since they play no role in the graph construction.Hence, what we have defined up to here is a generator of N -body polynomials or PIPs.Following the general procedure, we are now free to set the desired maximum polynomial degree and produce the combinations of monomial exponents l to generate a finite set of functions.

Appendix B. Method for the plasma model
The method used to model the plasma is similar to that previously described [42,100], in which each particle's level i density n i can be solved from a continuity equation where D denotes diffusion coefficient and the right hand side is the net particle source and sink term by collisional-radiative (CR) process and particle flux into and out a volume.In steady state of ∂n i /∂t = 0, the density balance equations for a long cylindrical vessel plasma can be expressed as follows.For atomic levels of H (n i , i = 1 − 40) The density of H + , n H + is deduced from the quasi neutrality condition The electron density, n e , is assumed to have a radial distribution close to a Besseltype profile for the ambipolar-diffusion regime considered here and applying the Bohm criterion as boundary conditions between the plasma and the vessel walls with µ = 2.405 for an infinite cylinder of the effective radius of R [35,100].The effective R is set as 40 cm for our plasma device.Diffusion time τ for H atom in the device of radius R d is given by τ = 2R d /v th with the thermal velocity v th = 2 2T a /πM H for atomic temperature T a and the hydrogen mass M H .The wall recombination coefficient γ is given as the empirical expression γ = 0.151 exp(−1.09× 10 3 /T m ) [105], with T m being the temperature of H 2 molecule.
Under the ambipolar diffusion assumption [106] the diffusion coefficient D AH + for H + is given by where p = n H2 T m is in Torr and the electronic temperature T e is in eV.The reduced mobility K 0 1 (cm 2 V −1 s −1 ) for H + and D + are 15.9 and 11.2, respectively [107].The diffusion coefficients for H on the wall are taken from [108].The conversion coefficient γ from H atom to H 2 molecule on the wall is also taken from [108].Q in , Q and V denote the gas flow rate of H 2 in the inlet, the pumping rate in the outlet and the volume of the plasma device, respectively.Q in , Q and V are given by 600 sccm (standard cubic centimeters per minute), 4800 lps (liter per second) and 2.64×10 6 cm 3 , respectively.
The CR processes and the rate coefficient α notations are listed in Table B1 with the references of the collision cross section and the radiative transition rate.The escaping factor η for radiation trapping in optically thick plasma is given as that for the infinite long cylinder [109].[98] Table B1.Collisional-radiative reactions and the rate coefficients considered in our plasma CRM

Reaction
From the cross section the rate coefficient is obtained as ) are solved by the multidimensional secant Broyden's method [113] setting the initial n i as the solution of the linear part of Eqs.B2 and B3 for various T e , n e , T m and T a .

Figure 1 .
Figure 1.The top left panel shows the energy distribution of training data in three regions of the configuration space, as described in the text.Top right and bottom panels show the energy difference between the three PES considered in this work and ab initio calculations.

1 Figure 2 .
Figure 2. H 2 + H + 2 approaches for several θ 1 and θ 2 angles and r 1 = 0.740 Å (for H 2 ) and r 2 = 1.055Å (for H + 2 ) the equilibrium distances of both species.The inset in the right hand side shows the coordinates used.

Figure 3 .
Figure 3. H + 3 + H approaches, as a function of d, the distance of H to the H + 3 center-of-mass.The orientations are preserved in the columns.The rows correspond to different equilateral H + 3 bond distances: top panels 0.85 Å (equilibrium), middle panels to 0.7 Å and bottom panels to 0.6 Å, respectively.In the top panel, corresponding to H + 3 equilibrium configuration, the asymptotic energy is -1.816 eV.

For H + 2 (v = 6 )
the double fragmentation (DF) channel opens at ≈ 2 eV, as it is shown in the left panel of Figure 8.The opening of this channel occurs approximately for values where the extra energy of v =6, 1.67 eV, is added to the D 0 of H + 2 , 2.65 eV If we add the DF channel contribution to the production of H + 3 , as shown in the right panel of Figure 8 the cross section increases considerably, reaching a very good agreement with the new experimental data of Savic et al. [1] above 3 eV.

Figure 9 .
Figure 9. Two typical trajectories leading to double fragmentation (DF), as a function of the collision time in fs.Lower panels show the charge on each atom (Mulliken population) for the ground electronic state.Middle panels shows the energies of the four lower electronic states.Upper panels show some characteristic internuclear distances needed to characterize the trajectory.R ij refers to the distance between atom H i and H j

Figure 10 .
Figure 10.Vibrational, rotational and translational energy distributions of the H + 3 + H products, as a function of the collision energy for the H 2 (v = 0) + H + 2 (v ) collisions, for v =0 (left panel) and v =6 (right panel).The origin of energy is in the bottom of the well of the H + 3 products.The initial vibrational energies of the reactants is 0.413 and 1.673 eV for the v = 0 and 6, respectively, with respect to the minimum of each fragment.The potential energy difference between reactants and products is 1.816 eV, and when ZPE are accounted for, the exoergicity of the present PES becomes 1.688 eV.

2 Figure 11 .
Figure 11.The population densities of D and H species by modeling using the different cross sections shown in Figure6.In the inset, the labels of x-axis correspond to the cross sections used given by the upper number in the legends of the main graph.

+ 2 and H + 3 2 :
are given with the relation D AH + : D AH + conversion coefficients from H ions into H atom (ζ a ) and H 2 molecule (ζ m ) B6)as the averaging over the relative velocity v 12 distribution.When the colliding particles of the masses m 1 and m 2 have the Maxwellian energy distribution with temperatures T 1 and T 2 the rate coefficient α can be expressed as[112] 12 ) exp(−(v 12 /v T12 ) 2 )v3  12 dv 12 , (B7)whereT 12 = (m 2 T 1 + m 1 T 2 )/(m 1 + m 2 ) and v T12 = 2(m 1 + m 2 )T 12 /m 1 m 2 for temperatures in eV.For electron collisions T 12 = T e and v T12 = 2T e /m e , with M e being the electron mass.The Maxwellian rate coefficient α 16 for the heavy particle collision H 2 + H + 2 → H + 3 + H(n = 1) is obtained with T 1 = T 2 = T m and m 1 = m 2 = 2M H .The balance equations of dn i /dt = 0 including the nonlinear terms of η ij (n i ), α 16 n 41 n 42 and D A (n 41

Table 1 .
Root mean squared errors of the ground electronic state for different energy regions, defined for E < Ec.The number of points on which the RMSE is calculated is presented.The zero of energy is set at H 2

Table 2 .
Root mean squared errors of the ground electronic state for different channels, calculated on all the geometries with energy lower than 17 eV.