# Non-adiabatic Quantum Dynamics of the Dissociative Charge Transfer He^{+}+H_{2} → He+H+H^{+}

^{1}Consiglio Nazionale delle Ricerche, Istituto di Struttura della Materia, Rome, Italy^{2}Departamento de Química Física Aplicada, Facultad de Ciencias, Universidad Autónoma de Madrid, Madrid, Spain^{3}Consiglio Nazionale delle Ricerche, Istituto per i Processi Chimico Fisici, Pisa, Italy

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 H_{2} vibrational states, *v*_{0} = 0 and 1, in the ground rotational state (*j*_{0} = 0) are obtained at collision energies *E*_{coll} ≤ 3 eV. We employ the lowest two excited diabatic electronic states of ${\text{HeH}}_{2}^{+}$ 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 *E*_{coll}, whose amplitudes decrease with the energy, due to an ion-induced dipole interaction in the entrance channel. At higher *E*_{coll}, reaction probabilities and cross sections increase monotonically up to 3 eV, remaining however quite small. When H_{2} is in the *v*_{0} = 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} cm^{3}s^{−1} at *v*_{0} = 0 and 1, respectively. Test calculations for H_{2} at *j*_{0} = 1 show that the probabilities can be enhanced by a factor of ~1/3, that is *ortho-*H_{2} seems ~4 times more reactive than *para-*H_{2}.

## 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}, ${\text{H}}_{2}^{+}$ (Stancil et al., 1993), and HeH^{+} (Zygelman et al., 1998) and then atom+diatom bimolecular collisions as He+${\text{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+${\text{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+${\text{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) $\stackrel{~}{X}$^{2}*A*' of ${\text{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 Ã^{2}*A*' and $\stackrel{~}{B}$^{2}*A*' of ${\text{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*' and $\stackrel{~}{B}$^{2}*A*' adiabatic and (1)^{2}*A*_{1} and (2)^{2}*B*_{2} diabatic PESs *V* of ${\text{HeH}}_{2}^{+}$, where all chemical species are in the ground electronic states, save ${\text{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 ${\text{H}}_{2}^{+}$(*A*^{2}Σ${}_{u}^{+}$) curves at *r* = 2.19 *a*_{0} and *V* = 1.28 eV.

**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) and $\stackrel{~}{B}$^{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 ${\text{H}}_{2}^{+}$(*A*^{2}Σ${}_{u}^{+}$, unbound).

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 1*s* 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.

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 infinite-order-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.

## Theory and Calculations

### Potential Energy Surfaces and Non-adiabatic Coupling

As we said in the Introduction, the ${\text{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*' and $\stackrel{~}{B}$^{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 *V*_{A} and *V*_{B} of a 2 × 2 matrix $\text{}\left(\begin{array}{cc}{V}_{11}& {V}_{12}\\ {V}_{12}& {V}_{22}\end{array}\right)\text{}$, 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*' and $\stackrel{~}{B}$^{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, 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,

fulfills the anti-symmetric requirement.

The DCT is produced through the CI, as shown in Figure 2 using nuclear Jacobi coordinates *R, r* = *R*${}_{\text{H}{\text{H}}^{\text{\u2032}}}$, and γ , at *R* = 4 *a*_{0} and γ = 0. The top panel shows the adiabatic Ã^{2}*A*' and $\stackrel{~}{B}$^{2}*A*' PESs of ${\text{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*)+${\text{H}}_{2}^{+}$(${A}^{2}{\Sigma}_{u}^{+}$, unbound) for *r* > 2 *a*_{0}. To analyze the accuracy of the fit, we compare the analytical non-adiabatic 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 ${\text{H}}_{4}^{+}$ and ${\text{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, 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.

**Figure 2**. Top panel: Ã^{2}*A*' and $\stackrel{~}{B}$^{2}*A*' adiabatic PESs at reactant Jacobi coordinates *R* = 4 *a*_{0}, *r*, and γ = 0. The vibrational *v*_{0} = 0 and 1 probability density is plotted at the corresponding vibrational energy. Bottom panel: Comparison of the absolute values of the fitted and *ab initio* Multi-Reference Configuration-Interaction (MRCI) NACME.

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*${}_{\text{H}{\text{H}}^{\text{\u2032}}}$ 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).

**Figure 3**. Diabatic PESs: *V*_{12}/(*V*_{22}-*V*_{11}) as function of the nuclear distances *R*_{HeH} and *R*_{HH}' = *r*, and of the included angle θ . Solid/dashed curves correspond to positive/negative values.

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 ${\text{HeH}}_{2}^{+}$ spinless rovibronic Hamiltonian Ĥ, which contains the electronic Ĥ

^{el}and the total angular momentum $\widehat{J}$ 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, ℏ

*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 2

*J*+1 values of

*K*are thus factorized in two non-interacting 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 ${\text{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} > (Gray and Balint-Kurti, 1998). 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

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}_{abs}^{R}$(*R*–*R*_{abs})^{2}] and $exp[{C}_{abs}^{r}{(r-{r}_{\text{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}_{{v}_{0}}^{J}$(*E*_{coll}), with the initial WP on the ${\text{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 coupled-channel 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.

## Collision Results

### Reaction Probabilities

We plot in Figure 4 the opacities functions (2*J*+1)${P}_{{v}_{0}}^{J}$(*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 ${\text{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 *K*s 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 nano-second. 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 one-two 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.

**Figure 4**. *v*_{0} = 0, 1; *J* = 20, 40; *kstep* = 5. Convergence of the opacities (2*J*+1)${P}_{{v}_{0}}^{J}$ with respect to *K*_{max}.

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.

**Figure 5**. *v*_{0} = 0; *J* = 10, 20; *kstep* = 120. Convergence of the opacities (2*J*+1)${P}_{{v}_{0}}^{J}$ with respect to *K*_{max}.

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.

**Figure 6**. *v*_{0} = 1; *J* = 5, 10, 15, 20; *K*_{max} = 3. Convergence of the opacities (2*J*+1)${P}_{{v}_{0}}^{J}$ with respect to *kstep*.

### 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 0.05 eV. The results in Figure 7 confirm the general scenario of the opacities: sharp resonances below 0.3 eV and a monotonic increase at larger *E*_{coll}. The *J* sum now merges the resonances 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 (1994, 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.

**Table 3**. Present quantum ICSs/Å^{2} vs. those semiclassical (Aguillon, 1998).

Finally, we report in Table 4 and Figure 8 the initial-state-resolved rate constants *k*_{v0}(*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 rate-constant maximum was also found in the adiabatic ground PES dynamics of H+HeH^{+} (Esposito et al., 2015) using the PES of 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 *para*- and *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 ${\text{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 Ã^{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).

**Figure 9**. *v*_{0} = 0, *J* = 0, *kstep* = 120. Adiabatic Ã^{2}*A*' (black) and non-adiabatic (1)*A*_{1}/(2)*B*_{2} (red) probabilities.

## 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 ${\text{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 charge-transfer formation (Mrugala et al., 2013) of stable ${\text{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, 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.

## Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

## 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-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. FIS2014-52172-C2-2-P and FIS2017-83473-C2-2-P. CINECA is also acknowledged for computer time awarded via the ISCRA programme.

## References

Aguado, A., and Paniagua, M. (1992). A new functional form to obtain analytical potentials of triatomic molecules. *J. Chem. Phys.* 96, 1265–1275. doi: 10.1063/1.462163

Aguado, A., Suárez, C., and Paniagua, M. (1993). Accurate fit of the two lowest excited-state potential-energy surfaces for doublet ${\text{HeH}}_{2}^{+}$. *J. Chem. Phys.* 98, 308–315.

Aguado, A., Suárez, C., and Paniagua, M. (1994). Accurate global fit of the H_{4} potential energy surface. *J. Chem. Phys.* 101, 4004–4010. doi: 10.1063/1.467518

Aguado, A., Tablero, C., and Paniagua, M. (1998). Global fit of *ab initio* potential energy surfaces I. Triatomic systems. *Comp. Phys. Comm*. 108, 259–266. doi: 10.1016/S0010-4655(97)00135-5

Aguado, A., Tablero, C., and Paniagua, M. (2001). Global fit of *ab initio* potential energy surfaces: II.1. Tetraatomic systems ABCD. *Comp. Phys. Comm*. 134, 97–109. doi: 10.1016/S0010-4655(00)00181-8

Aguillon, F. (1994). Semi-classical coupled wavepacket study of the dissociative charge exchange He^{+}+H_{2}→ He+H+H^{+}. *Chem. Phys. Lett.* 222, 69–74. doi: 10.1016/0009-2614(94)00330-0

Aguillon, F. (1998). A new treatment of nonadiabatic dynamics: application to the determination of the He^{+}+H_{2} → He+H+H^{+} differential cross section. *J. Chem. Phys.* 109, 560–571. doi: 10.1063/1.476592

Aquilanti, V., Cavalli, S., De Fazio, D., Simoni, A., and Tscherbul, T. V. (2005). Direct evaluation of the lifetime matrix by the hyperquantization algorithm: narrow resonances in the F+H_{2} reaction dynamics and their splitting for nonzero angular momentum. *J. Chem. Phys*. 123:054314. doi: 10.1063/1.1988311

Aquilanti, V., Cavalli, S., Simoni, A., Aguilar, A., Lucas, J. M., and De Fazio, D. (2004). Lifetime of reactive scattering resonances: Q-matrix analysis and angular momentum dependence for the F+H_{2} reaction by the hyperquantization algorithm. *J. Chem. Phys.* 121, 11675–11690. doi: 10.1063/1.1814096

Bovino, S., Tacconi, M., Gianturco, F. A., and Galli, D. (2011). Ion chemistry in the early unverse: revisiting the role of HeH^{+} with new quantum calculations. *Astronomy Astrophys.* 529:A150. doi: 10.1051/0004-6361/201116740

Cederbaum, L. S., and Friedman, R. S. (2003). Conical intersections and bound molecular states embedded in the continuum. *Phys. Rev. Lett*. 90:013001. doi: 10.1103/PhysRevLett.90.013001

Coppola, M., D'Introno, R., Galli, D., Tennyson, J., and Longo, S. (2012). Non-equilibrium H_{2} formation in the early universe: energy exchanges, rate coefficients, and spectral distortions. *Astrophys. J. Suppl. Ser.* 199:16. doi: 10.1088/0067-0049/199/1/16

Coppola, M., Longo, S., Capitelli, M., Palla, F., and Galli, D. (2011). Vibrational level population of H_{2} and ${\text{H}}_{2}^{+}$ in the early universe. *Astrophys. J. Suppl. Ser.* 193:7. doi: 10.1088/0067-0049/193/1/7

De Fazio, D. (2014). The H+HeH^{+} → He+${\text{H}}_{2}^{+}$ reaction from the ultra-cold regime to the three-body breakup: exact quantum mechanical integral cross sections and rate constants. *Phys. Chem. Chem. Phys.* 16, 11662–11672. doi: 10.1039/C4CP00502C

De Fazio, D. M., de Castro-Vitores, M., Aquilanti, V., and Cavalli, S. (2012). The He+${\text{H}}_{2}^{+}$ → 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.* 137:244306. doi: 10.1063/1.4772651

Dhuicq, D., Jugi, B., Benoit, C., and Sidis, V. (1998). The He^{+}+H_{2} → HeH^{+}+H*(2s,2p)+H_{2} reaction: study of the *s* and *p* contributions in a coincidence-correlation experiment. *J. Chem. Phys.* 109, 512–524. doi: 10.1063/1.476588

Dhuicq, D., Lehner, O., Linder, F., and Sidis, V. (1996). Angular and energy analysis of HeH^{+} products from the charge exchange excitation reaction He^{+}+H_{2} → HeH^{+}+H*. *Chem. Phys.* 206, 139–159. doi: 10.1016/0301-0104(96)00014-6

Dunning, T. H. Jr. (1989). Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen. *J. Chem. Phys.* 90, 1007–1023.

Esposito, F., Coppola, C. M., and De Fazio, D. (2015). Complementarity between quantum and classical mechanics in chemical modeling. The H+HeH^{+} → ${\text{H}}_{2}^{+}$+He reaction: a rigourous test for reaction dynamics methods. *J. Phys. Chem. A* 119, 12615–12626. doi: 10.1021/acs.jpca.5b09660

Galli, D., and Palla, F. (2013). The dawn of chemistry. *Annu. Rev. Astron. Astrophys*. 51:163. doi: 10.1146/annurev-astro-082812-141029

Gamallo, P., Akpinar, S., Defazio, P., and Petrongolo, C. (2013). Conical-intersection quantum dynamics of OH(*A*^{2}Σ^{+})+H(^{2}*S*) collisions. *J. Chem. Phys.* 139:94303. doi: 10.1063/1.4819355

Gamallo, P., Akpinar, S., Defazio, P., and Petrongolo, C. (2015). Born–Oppenheimer and Renner–Teller quantum dynamics of CH(X^{2}Π)+D(^{2}S) reactions on three CHD potential surfaces. *J. Phys. Chem. A* 119, 11254–11264. doi: 10.1021/acs.jpca.5b08891

Giese, C. F., and Maier, W. M. I. I. (1961). Ion-molecule reactions studied with mass analysis of primary ion beam. *J. Chem. Phys.* 35, 1913–1914. doi: 10.1063/1.1732184

Gioumousis, G., and Stevenson, D. P. (1958). Reactions of gaseous molecule ions with gaseous molecules. V. Theory. *J. Chem. Phys.* 29, 294–299. doi: 10.1063/1.1744477

Gray, S. K., and Balint-Kurti, G. G. (1998). Quantum dynamics with real wave packets, including application to three-dimensional (*J* = 0)D+H2 → HD+H(*J* = 0) reactive scattering. *J. Chem. Phys*. 108, 950–962. doi: 10.1063/1.475495

Guo, H. (2012). Quantum dynamics of complex-forming bimolecular reactions. *Int. Rev*. *Phys. Chem.* 31, 1–68. doi: 10.1080/0144235X.2011.649999

Herbst, E., and Klemperer, W. (1973). The formation and depletion of molecules in dense interstellar clouds. *Ap. J.* 185:505. doi: 10.1086/152436

Johnsen, R., Chen, A., and Biondi, M. A. (1980). Dissociative charge transfer of He^{+} ions with H_{2} and D_{2} molecules from 78 to 330 K. *J. Chem. Phys.* 72, 3085–3088. doi: 10.1063/1.439512

Lepp, S., Stancil, P. C., and Dalgarno, A. (2002). Atomic and molecular processes in the early universe. *J. Phys. B At. Mol. Opt. Phys.* 35, R57–R80. doi: 10.1088/0953-4075/35/10/201

Martínez, E., Grimbert, D., Aguillon, F., and Sidis, V. (2000). Calculations of quantum cross-sections for dissociative charge transfer in the He^{+}+H_{2} collision in the 10–50 eV center of mass energy range. *Chem. Phys. Lett.* 322, 103–110. doi: 10.1016/S0009-2614(00)00363-8

McGuire, P., and Kouri, D. J. (1974). Quantum mechanical close coupling approach to molecular collisions. jz-conserving coupled states approximation. *J. Chem. Phys*. 60, 2488–2499.

McLaughlin, D. R., and Thompson, D. L. (1979). Ground- and lower excited-state discrete *ab initio* electronic potential energy surfaces for doublet ${\text{HeH}}_{2}^{+}$. *J. Chem. Phys.* 70, 2748–2769.

Meijer, J. H. M., Goldfield, E. M., Gray, S. K., and Balint-Kurti, G. G. (1998). Flux analysis for calculating reaction probabilities with real wave packets. *Chem. Phys. Lett.* 293, 270–276. doi: 10.1016/S0009-2614(98)00743-X

Mrugala, F., and Kraemer, W. P (2013). Radiative charge transfer in He^{+} +H_{2} collisions in the milli- to nano-electron-volt range: a theoretical study within state-to-state and optical potential approaches. *J. Chem. Phys.* 138:104315. doi: 10.1063/1.4793986

Pack, R., Walker, R., and Kendrick, B. (1998). Three-body collision contributions to recombination and collision-induced dissociation. I. Cross sections. *J. Chem. Phys.* 109, 6701–6713. doi: 10.1063/1.477349

Pack, R. T. (1974). Space-fixed vs body-fixed axes in atom-diatomic molecule scattering. Sudden approximations. *J. Chem. Phys*. 60, 633–639. doi: 10.1063/1.1681085

Petrongolo, C. (1988). Nonadiabatic theory of triatomics: general formalism and application to Renner-Teller and conical-intersection effects. *J. Chem. Phys.* 89, 1297–1308. doi: 10.1063/1.455181

Preston, R. K., Thompson, D. L., and McLaughlin, D. R. (1978). A theoretical prediction of vibrational enhancement for dissociative charge transfer in the ${\text{HeH}}_{2}^{+}$ system. *J. Chem. Phys.* 68:13.

Ramachandran, N., De Fazio, D., Cavalli, S., Tarantelli, F., and Aquilanti, V. (2009). Revisiting the potential energy surface for the He+${\text{H}}_{2}^{+}$ → HeH^{+}+H reaction at the full configuration interaction level. *Chem. Phys. Lett.* 469, 26–30. doi: 10.1016/j.cplett.2008.12.035

Rydberg, R. (1931). Rydberg-Klein-Rees method for determining potential energy curves of diatomic molecules from band spectra. *Ann. Phys*. 73:376.

Sanz-Sanz, C., Aguado, A., Roncero, O., and Naumkin, F. (2015). Non-adiabatic couplings and dynamics in proton transfer reactions of ${\text{H}}_{n}^{+}$systems: application to H_{2}+${\text{H}}_{2}^{+}$→*H*+*H*${}_{3}^{+}$ collisions. *J. Chem. Phys.* 143:234303. doi: 10.1063/1.4937138

Schauer, M., Jefferts, S. R., Barlow, S. E., and Dunn, G. H. (1989). Reactions of H_{2} with He^{+} at temperatures below 40 K. *J. Chem. Phys.* 91, 4593–4596. doi: 10.1063/1.456748

Sidis, V. (1996). Diabatic excited states of the (HeH2)^{+} molecular ion for the charge exchange-excitation reaction: He^{+}+H_{2} → HeH^{+}+H*. *Chem. Phys.* 209, 313–326. doi: 10.1016/0301-0104(96)00060-2

Stancil, P. C., Babb, J. F., and Dalgarno, A. (1993). The radiative association of He^{+} and He and H^{+} and H. *Ap. J.* 414:672. doi: 10.1086/173113

Stedeford, J. B. H., and Hasted, J. B. (1955). Further investigations of charge exchange and electron detachment—I. Ion energies 3 to 40 keV. *Proc. Royal Soc. A* 227, 466–486.

Szidarovszky, T., and Yamanouchi, K. (2016). Photodissociation dynamics of weakly bound ${\text{HeH}}_{2}^{+}$ in intense light fields. *Phys. Rev. A* 94:063405. doi: 10.1103/physreva.94.063405

Tennyson, J., and Miller, S. (1987). Predicted vibrational-rotation levels of H_{2}He^{+} and its isotopomers. *J. Chem. Phys.* 87, 6648–6652. doi: 10.1063/1.453399

Werner, H.-J., Knowles, P. J., Knizia, G., Manby, F. R., Schütz, M., Celani, P., et al. (2018). *MOLPRO, Version 2015.1, a Package of ab initio Programs*. Available online at: http://www.molpro.net

Woon, E., and Dunning, T. H. Jr. (1994). Gaussian basis sets for use in correlated molecular calculations. IV. Calculation of static electrical response properties. *J. Chem. Phys.* 100, 2975–2988.

Keywords: He^{+}+H_{2}, wavepacket, conical intersection, non-adiabatic, quantum, dynamics

Citation: De Fazio D, Aguado A and Petrongolo C (2019) Non-adiabatic Quantum Dynamics of the Dissociative Charge Transfer He^{+}+H_{2} → He+H+H^{+}. *Front. Chem.* 7:249. doi: 10.3389/fchem.2019.00249

Received: 08 February 2019; Accepted: 27 March 2019;

Published: 16 April 2019.

Edited by:

Antonio Aguilar, University of Barcelona, SpainReviewed by:

Tomás González-Lezana, Spanish National Research Council (CSIC), SpainRicardo Gargano, Universidade de Brasília, Brazil

Copyright © 2019 De Fazio, Aguado and Petrongolo. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Dario De Fazio, defazio.dario@yahoo.it