Atom–Diatom Reactive Scattering Collisions in Protonated Rare Gas Systems

The study of the dynamics of atom–diatom reactions involving two rare gas (Rg) atoms and protons is of crucial importance given the astrophysical relevance of these processes. In a series of previous studies, we have been investigating a number of such Rg(1)+ Rg(2)H+→ Rg(2)+ Rg(1)H+ reactions by means of different numerical approaches. These investigations comprised the construction of accurate potential energy surfaces by means of ab initio calculations. In this work, we review the state-of-art of the study of these protonated Rg systems making special emphasis on the most relevant features regarding the dynamical mechanisms which govern these reactive collisions. The aim of this work therefore is to provide an as complete as possible description of the existing information regarding these processes.


Introduction
The presence of rare gas (Rg) atoms and H + in the early stages of the formation of the Universe and their abundance in the interstellar medium (ISM) [1][2][3][4][5][6][7] explains the interest on understanding the relevant features regarding the chemistry involving protonated Rg species. Thus, for example, the helium hydride ion, HeH + , is one of the first bond formed by radiative association in the primordial nucleosynthesis, playing a crucial role in the formation of molecular hydrogen. Its existence was predicted a long time ago, but it has been only in recent years that its detection has been finally confirmed [8]. The presence of ArH + in Crab nebula has been inferred by spectra recorded by the Herschel Space Observatory [9]. Reactive collisions among these species are usually characterized by the absence of potential energy barriers and are therefore possible at the low temperature regime existing in the ISM. The study of the dynamics governing those processes can provide useful information about the evolution of the stellar medium as revealed by the growing interest of these systems in an astrophysical context [10][11][12].
Early accurate ab initio calculation of structure and energies [40] for He 2 H + were followed by the development of PESs [34,35] and the study of bound and quasibound states of the He 2 H + and He 2 D + complexes [34]. QM dynamical investigations of the He + HeH + → HeH + +He reaction by means of time-dependent wave packet (TDWP) calculations on these PESs [29,34,53] revealed noticeable differences between the probabilities obtained with the Coriolis coupled (CC) and centrifugal sudden (CS) approaches.
In this work, we review theoretical state-of-the-art investigations on systems involving two Rg atoms and H + performed by us with a number of methods. The studies involved the development of ab initio PESs describing the existing interactions between the colliding atoms and the application of different numerical techniques to study the dynamics of the process. Both QCT and QM approaches are comparatively applied within their different frameworks: Gaussian and histogram binning for the former and the already mentioned CC and CS alternatives for the latter. Given that these reactions are mediated by the presence of a relatively deep potential well in the intermediate region between reactants and products, a statistical quantum method (SQM) is also employed to test the possible importance of complex-forming mechanisms on the overall dynamics of the process. Our main goal is therefore to provide a complete compilation of recent studies on these reactions, thus giving a detailed description of how to tackle the dynamical analysis of this particular atom-diatom collisions.
The structure of this work is the following. In Section 2, we describe the procedure followed to obtain the PES, including first the ab initio methods (Section 2.1) and second the fitting to an analytical expression (Section 2.2). In Section 3, we discuss the details of the different methods employed to study the dynamics of these reactions. In Sections 3.2-3.6, we show the examples of the He + HeH + → HeH + + He, He + NeH + → HeH + + Ne, Ne + HeH + → NeH + + He, Ne + NeH + → NeH + + Ne and Ar + ArH + → ArH + + Ar, respectively. Finally, in Section 4 concluding remarks are presented.

Methods
In the literature, protonated Rg systems have been investigated at different levels of theory. However, the "gold standard" coupled-cluster singles, doubles, and perturbative triples (CCSD(T)) method with augmented cc-pVnZ (n = T, Q, 5) basis sets provide the best results compared to other levels. For He 2 H + , [HeHNe] + , Ne 2 H + and Ar 2 H + systems, analytical PESs were constructed using ab initio energies obtained using CCSD(T) electronic structure calculations. While the d-aug-cc-pVTZ basis set was used for He 2 H + , the augcc-pVQZ basis set was chosen for all the other systems while carrying out the ab initio calculations. The ab initio calculations included the computation of energies for both the corresponding triatomic and three diatomic fragments.
For the construction of PESs for the triatomic system ABC, ab initio energies were computed along the grids defined in internal coordinates (r AB , r BC and θ), where r AB and r BC are the two internuclear distances between hydrogen and Rg atoms and θ is the angle between r AB and r BC as shown in Figure 1. Adiabatic PESs for ABC type systems can be represented as a many-body expansion function [66]: ABC (r AB , r BC , r AC ). (1) Here, V ABC (r AB , r BC , r AC ) is the potential energy of the total system, V (1) i s are the energies of the free atoms in their corresponding states, V (2) i s are two-body (2B) interaction energies, and V Coordinate considered for ab initio electronic structure calculations and dynamical simulations. r AB , r BC , and r AC are the distances between atoms A and B, B and C, and A and C, respectively.
θ is the angle between r AB and r BC . While the set r AB , r BC and θ is used for the ab initio points, the Jacobi coordinates set of R, r and γ is used for the propagation, where R is the distance from the center-of-mass of BC to A and γ is the angle between R and r.
To represent the 2B interaction energies, analytical PESs for all the possible diatomic fragments were obtained using either cubic spline interpolation or nonlinear curve fitting methods. For the nonlinear curve fitting, the diatomic terms were expressed using a polynomial form proposed by Aguado and Paniagua [67]: where M is the order of the polynomial, ρ AB = r AB e −β AB r AB , c i s are linear parameters, and α and β are nonlinear parameters. The 3B interaction energy terms are expressed as [67] V where d ijk are the linear parameters and ρ i = r i e −β i r i . The constraints i + j + k = i = j = k and i + j + k ≤ M are used in Equation (3) to make the 3B interaction energies zero at asymptotic limits. The Levenberg-Marquardt nonlinear optimization algorithm [68] was used to determine the linear and nonlinear parameters in Equation (3).

Analytical PESs
Global analytical PESs for RgH + , Rg 2 H + , or [RgHRg'] + systems constructed from ab initio energies are reported in many instances in the literature. Details of electronic structure calculation methods and basis sets, number of ab initio energies used, and root mean square errors (RMSE) for the analytical PESs discussed in this review are tabulated in Table 1 as reported in the original references. A large number of ab initio energies are used in the fitting procedures except for He 2 H + . The fitting errors are less than 0.1 kcal/mol for all the cases, which guarantees high quality of those analytical PESs. For He 2 H + , apart from the PES reported in Table 1, another analytical PES has also been constructed using 15682 energies calculated at multi-reference configuration interaction and d-aug-cc-pV5Z level of theory [35] with an RMSE of 0.048 kcal/mol. Similarly, a global PES for Ar 2 H + has been reported using 7040 QCISD/6-311++G(3df,3pd) energies with fitting RMSE of 0.143 kcal/mol [39]. Color maps of the analytical PESs reported in Table 1 are shown in Figure 2 for different possible reaction channels. Figure 2 shows the existence of deep potential well regions in all the systems for collinear and near collinear geometries with the hydrogen atom positioned between two Rg atoms. Geometries for the global minima for all the systems correspond to linear configurations, which are symmetric for Rg 2 H + systems. Energies of the global minima and the possible reactant/product channels for the PESs are given in Table 2. The depth of the global minima from the nearest reactant asymptotes lie between ∼12-16.5 kcal/mol, with the lowest depth (12.014 kcal/mol) for [HeHNe] + and highest depth for (16.605 kcal/mol) for Ne 2 H + . Hence, these triatomic cation systems are quite stable unlike noncovalent weak van der Waals complexes. The RgH + bond distances for the equilibrium configurations of the triatomic systems and for the diatomic molecular ions are also given in Table 2. As it is observed, the internuclear separations between H and Rg atoms in the triatomic complexes get elongated compared to the free diatomic molecular ions. It is also seen that the He-H bond distance is longer in [HeHNe] + than in He 2 H + . On the other hand, the Ne-H bond distance is longer in Ne 2 H + compared to the bond in [HeHNe] + . This shows that size of the Ne atom plays an important role in this case.
In cases of Ne + HeH + and Ar + ArH + reactive systems, the existence of other interesting stable species is also observed in Figure 2. Shallow potentials can be seen in Ne+HeH + and Ar+ArH + color maps when the free Rg atom approaches the diatomic molecular ion from the other Rg atom side. This suggests formation of weak Rg· · · RgH + complexes. The bond distance of one RgH is very close to the free RgH + . For [HeHNe] + , the second minimum (Ne· · · HeH + ) is 1.36 kcal/mol lower than the Ne+HeH + asymptote, whereas for Ar 2 H + , the other minimum (Ar· · · ArH + ) is 2.74 kcal/mol lower in energy compared to the Ar+ArH + asymptote. Although existence of a second minimum was reported for Ne 2 H + in Reference [40], no such minimum was found for Ne 2 H + in [70]. 6    Minimum energy paths for possible atom+molecular ion collisions for these triatomic systems are shown for different internal angles in Figure 3. It can be seen that all the reactions are barrierless in their entrance channels. Collinear and near-collinear approaches are most favorable for these reactions and a potential well with depth ∼12-16.5 kcal/mol can be seen in the strong interaction region. However, for small attacking angles, barriers emerge in the same regions due to strong repulsion between the Rg atoms. The barrier height increases with the size of rare gas atoms, and with decreasing attacking angles.

Methods
In this section, we summarize the most relevant aspects of the theoretical methods employed in the study of the title reactions.

Time-Dependent Quantum Mechanical Method
Details of time-dependent (TD) wave packet (WP) methodology is well documented, and here we provide a brief outline only. In this method, the TD Schrödinger equation is solved on a grid using WPs. A set of body-fixed reactant Jacobi coordinates, as shown in Figure 1, is used in the calculations. An initial WP (either complex or real) prepared in the reactant asymptotic region is propagated using either the split operator [72] or Chebyshev real WP propagation [73,74] methods. The initial WP is constructed on equidistant grids along the radial (R, r) coordinates and Gauss-Legendre quadrature points along the angular degree of freedom. In this grid representation, the full Hamiltonian is written as a tridiagonal matrix as [75,76] Here, V(R, r, γ) is the PES, j is the rotational quantum number of the reactant diatom, J is the total angular momentum of the reactants, and K is the projection of J on the bodyfixed z-axis. The reduced masses along R and r coordinates are µ R and µ r , respectively. λ in Equation (4) is defined as λ ± JK = J(J + 1) − K(K ± 1). Reactive scattering calculations are performed using either the full Hamiltonian as given in Equation (4) including the CC or within CS approximation where the off-diagonal terms in the Equation (4) are neglected, and only K = 0 is considered.
In the WP propagation, the action of radial kinetic energy operators are evaluated using the fast Fourier transformation technique and the action of the angular momentum operator is evaluated in associated Legendre polynomial basis. After sufficient propagation, the energy dependent total reaction probabilities are calculated by summing the total flux passing through a fixed surface located at a large distance in the product channel. Spherical Hankel functions are used to compute the energy weightage of the initial translational WP.

Time-Independent Quantum Mechanical Method
In order to study the reactive scattering for the proton exchange reactions within the time-independent approach, the ABC program [77] was used. The methodology for the time-independent quantum mechanical (TIQM) approach followed in the ABC package is documented in [77] and references therein. The time-independent Schrödinger equation is solved using the CC approach in hyperspherical coordinates. The diatomic ro-vibrational wavefunctions of all the available reactant/product channels for the given maximum energy and j max (maximum number of rotational states allowed in each channel) are used to construct the CC basis functions. The CC hyper-radial equations are then solved by using a constant reference potential log derivative method between ρ max to ρ min (ρ is the hyperradius) in n sec sectors. Finally, the S-matrix elements for a particular J with diatomic parity (p) and triatomic parity (P) for all the channels are calculated from the final log derivative matrix by applying scattering boundary conditions.

Quasiclassical Trajectory Calculation
The standard QCT methodology for atom-diatom collisions, discussed in details in [78][79][80][81][82], is followed to study the H + transfer processes between Rg atoms. Initial conditions for the trajectories are sampled using standard Monte Carlo sampling in reactant Jacobi coordinates. The impact parameters in cross sections and rate constant calculations are sampled following either the normal procedures or using stratified sampling scheme [79,81,82]. Twelve coupled Hamilton's equations of motion are then integrated numerically in reactant Jacobi coordinates using the fourth-order Runge-Kutta method. The total energy and total angular momentum were conserved up to sixth and tenth decimal places for all the trajectories. The ro-vibrational states for the reactant and product diatoms were determined following either the semiclassical theory of bound states or using the discrete variable representation based Colbert-Miller method [83]. The final quantum numbers of the product were assigned using both histogram binning (HB) and Gaussian binning (GB) methods [81,84,85].
The reaction probability for a selected initial ro-vibrational state and a given total angular momentum can be computed as where N r is the number of reactive trajectories and N tot is the total number of trajectories. The reaction cross section for a given initial ro-vibrational state is then calculated as where b max is the maximum impact parameter for which reactive collision can occur. The initial state selected differential cross sections (DCS) (dσ r /dΩ) is computed as where k = 2µ R E/h 2 and, N r (J, θ) and N tot (J) are the number of reactive trajectories scattered at an angle θ and the total number of trajectories run, respectively, for a given J. N r (J, θ) can be determined by using histograms along θ.

Statistical Quantum Mechanical Method
For those reactions which proceed via a complex-forming mechanism it is possible to apply statistical techniques. One of such methods is the statistical quantum mechanical (SQM) developed by Manolopoulos and coworkers [86,87]. A detailed description of the numerical details of this technique can be found in previous references [88][89][90][91][92][93][94], where the SQM approach has been employed to study atom-diatom reactions. Under the assumption of a complex-forming mechanism governing the overall dynamics of the process, the state-to-state reaction probability can be approximated as where v, j, and l refer to the diatomic vibrational, rotational, and orbital angular momentum quantum numbers, respectively, and l is the orbital angular momentum In the above expression, Equation (8), p J vjl (E) is the capture probability or probability of forming the collision complex from the rovibrational state vjl at the reactant arrangement at the total angular momentum J and the energy E, while p J v j l (E) is the probability for the collision complex to decay to the final v j l state of the product channel respectively. Indexes in the denominator run for all energetically accessible states. The above expression for the reaction probability in Equation (8) reveals that the SQM approach does not provide any information regarding the amplitude of scattering matrix and therefore it can reproduce any possible existing resonance structure, yielding exclusively to an average value.
These capture probabilities are obtained by solving a set of CC equations for each arrangement with the form where the interaction matrix W(R) is expressed as being µ the 3B collision mass [95] and E vj is the rovibrational energy of the diatom (v, j) state. The interaction potential matrix V J v j l ,vjl can be finally expressed in terms of the diatomic vibrational wavefunctions, vector-coupling coefficients and spherical harmonics as explained in Reference [86].
A computationally cheaper alternative to the CC scheme described in Equation (10) is the CS approximation where the coupled-channel equations are expressed in smaller sets for each value of K, the projection of the angular momentum on the atom-diatom axis, is obtained, where The solution of these equations is performed by means of a TIQM using the corresponding full ab initio PES within the region defined between asymptotic distances and a capture radius at which the collision complex is supposed to be formed. Therefore, the intermediate region where the PES for complex-forming processes usually display a relatively deep well is neglected. The above state-to-state probabilities are employed for the calculation of the corresponding integral cross sections (ICS) according to the following expression: with k 2 vj = 2µ(E − E vj )/h 2 , E vj being the energy of the initial rovibrational state vj of the reactant diatom and g the electronic degeneracy. Finally, the ICS of Equation (13) evaluated in the collision energy E c = E − E vj , is employed for the calculation of the rate constants as where β = (k B T) −1 .
The calculation of the differential cross section (DCS) by means of the SQM method requires an extra approximation due to the lack of information regarding the amplitude of the scattering matrix S vj,v j mentioned above in Equation (8): As a consequence of this expression, the statistical angular distributions are therefore, symmetric, predicting equal peaks both at the forward and backward scattering directions.

He + HeH + → HeH + + He
The He + HeH + (v = j = 0) → HeH + + He process was explored for the first time using the TDQM method [36] on the ab initio energy-based analytical PES developed in [34]. Reaction probabilities were computed for this reaction within the CS approximation which are shown in Figure 4 for some selected J values. It is worth mentioning that exact quantum dynamical simulations for this reaction has also been carried out by Xu and Zhang [29] using the PES constructed in [35]. To be consistent with the PESs mentioned in this review, exact quantum dynamics as well as QCT calculations have been performed, and total reaction probabilities obtained from TDQM-CC and QCT calculations are shown in Figure 4. The QM probabilities for J = 10 and 20 oscillate mostly around 0.2, but remain slightly larger at low energies. Few sharp peaks are obvious for J = 10 at low energies which suggests formation of a metastable He 2 H + complex in the potential well (see Figures 2 and 3). Probability curves for the high Js are relatively smoother with broad peaks. It was observed in [29] that the TDQM CS probabilities for this reaction computed using References [34,35] PESs are very similar which is also obvious in this work for the TDQM-CC probabilities when compared to the results presented in [29]. This suggests that although different level of theories were used to calculate the two PESs, their global topology does not differ significantly. As can be seen in Figure 4, the QCT method nicely reproduces the overall behavior of the exact QM results. Results obtained from both QCT binning schemes are quite similar. However, due to inherent zero point energy leakage in QCT, finite probabilities for this reaction are predicted for energies below threshold for higher Js.
Total ICSs as a function of collision energy and rate constants as a function of temperature for the He + HeH + (v = j = 0) → HeH + + He reaction are shown in Figure 5. TDQM-CS results in this figure are taken from Reference [36]. The ICSs for this reaction have high values at low collision energies, but their magnitudes decrease with the increase of collision energies. The TDQM-CS method underestimates the ICSs in the entire energy range while results seem to be overestimated by both QCT approaches up to energies about 0.3 eV. Beyond that energy, the two QCT alternatives agree quite well in reproducing the ICSs. QM ICSs, specially those obtained with the TDQM-CC approach, exhibit many oscillations as a function of the energy. Rate constants obtained using different dynamical methods for He + HeH + (v = j = 0) → HeH + + He reaction are found to be independent of temperature beyond 250 K, following the classical Langevin capture model [96,97] for a barrierless ion-molecule reaction. While the TDQM-CS method underestimates the reaction rate constants, predictions obtained with the QCT approaches remain above the QM results. It is also observed that both QCT binning schemes produce very similar results except at low temperatures.

He + NeH + → HeH + + Ne
One the characteristic features of reactions mediated by the presence of relatively deep potential wells is the existence of numerous resonance peaks in the probabilities as a function of the energy. In this particular case, the PES for the He + NeH + reaction exhibits a well of~16.14 kcal mol −1 for the [HeHNe] + species between reactants and products [53,69]. As shown in Figure 6 for the case of the He + NeH + → HeH + + Ne reaction, even for relatively large values of the total angular momentum, such as J = 30 and 50, the WP results from in [53] exhibit noticeable maxima attributed to the formation of an intermediate complex supporting a number of quasi bound and resonance states. The cases shown in Figure 6 correspond to two different initial reactant states, in particular, NeH + (v = 0, j = 0) and NeH + (v = 0, j = 1). Peaks, which are still narrow at J = 30, become however progressively broader as the value of the angular momentum increases, as revealed from the comparison with the case of J = 50. Differences between the CC and CS schemes for the WP calculation are seen for J = 30, but it is clear from Figure 6 that this effect is more important as J gets larger. The statistical predictions obtained by means of the SQM approach, also included in the figure for comparison, reveals certain independence with respect to the specific framework considered, as the CS and CC results do not seem too different. The description provided by the SQM values of the WP reaction probabilities remain as an acceptable average for J = 30 with no information regarding the resonance peaks, but for J = 50 the CC (CS) prediction overestimate (underestimate) the corresponding WP results.
The threshold for reaction is ≈ 0.29 eV for the case of the ground rovibrational state NeH + (v = 0, j = 0), a value which is increased for high values of the total angular momentum due to the increase of the centrifugal barriers. For this vibrationless case, thresholds exhibited by the QM-CC and QM-CS probabilities are different, and Figure 6 shows how the CS results have larger threshold energies. This feature is however not observed in the statistical case, with almost identical reaction probabilities for the SQM-CC and SQM-CS cases.
A similar comparison between QM WP and statistical predictions estimated by means of the SQM approach can be established through the opacity functions, that is, the reaction probability as a function of the total angular momentum J at a specific value of the collision energy. More precisely, Figure 7 shows the partial wave contribution (2J + 1)P J (E c ) as a function of J for the initial-state selected cases (v = 0, j = 0) and (v = 1, j = 0) at E c = 350 meV and 450 meV. As seen in Figure 6 for the case of the energy, the WP probabilities exhibit also an oscillatory behaviour with respect to J. Whereas significant differences are observed between WP-CC and WP-CS results, the SQM predictions show little dependence on the specific scheme, either CC or CS, employed in the calculations. The statistical estimations for the reaction initiated from the ground rovibrational state NeH + (v = 0, j = 0) provide a good average description of the corresponding QM WP probabilities. Furthermore, good agreement is also observed for the case of v = 1, with the SQM probabilities providing a correct reproduction of the existing decrease of the opacity functions at J ∼ 60 displayed by the QM-CS results. Partial waves coming from the larger J seem to contribute more noticeably in the QM-CC calculation rather than in the QM-CS approach where contributions coming from lower values of the total angular momentum are significant. The present theoretical analysis of the dynamics of the He + NeH + (v, j) → Ne + HeH + reaction includes the calculation of the ICSs for the (v = 0, j = 0) and (v = 1, j = 0) cases. In Figure 8, cross sections obtained with the WP-CS and WP-CC approaches are compared with the corresponding statistical predictions for a collision energy up to 0.5 eV. As suggested by the different comparison seen for the reaction probabilities (see Figure 7) for both initial rovibrational states between the WP and the SQM methods, the agreement for the (v = 0, j = 0) case contrasts with the differences observed with (v = 1, j = 0). Thus, whereas the energy threshold for the cross section when the reaction proceeds from the ground rovibrational state is well reproduced, the behaviour at the low energy regime exhibiting the characteristic trend of a barrierless process in the (v = 1, j = 0) case is clearly underestimated by the statistical prediction.

Ne + HeH + → NeH + + He
The existence of a relatively deep well (∼18.45 kcal mol −1 ) in the PES for this reaction also leads to the presence of narrow resonances in the corresponding probabilities. Examples of such reaction probabilities for three different values of J and three possible initial states HeH + (v, j): (v = 0, j = 0), (v = 0, j = 1), and (v = 1, j = 0) obtained with a WP method [63] are shown in Figure 9. Results from a QCT calculation, also included in the figure, reproduce in all cases the average trend of the WP probabilities although no information regarding the resonance peaks is recovered. The SQM predictions, in turn, despite to provide the correct threshold for reaction, overestimate noticeably the WP results as the energy increases for J = 10 and 40. For the highest value of the total angular momentum shown in Figure 9, the SQM probabilities seem to remain below both the WP and QCT results. These discrepancies between QM results and the corresponding statistical predictions are usually interpreted as deviations of the overall dynamics of the reaction from a purely complex-forming process [90,91,98,99].
The oscillatory behaviour of the time-dependent QM (TDQM) probabilities of Figure 9 disappears when we obtain the ICSs, washed out as a result of the partial waves averaging effect. The corresponding cross section for the (v = 0, j = 0) case between 10 −3 eV and 0.5 eV is shown in Figure 10. The comparison with the QCT and SQM results, also shown in the figure, reveals significant discrepancies with these two approaches at the lower energies [63]. The trend followed by the TDQM ICSs as the energy decreases is a consequence of the difficulties of the WP techniques at this regime, but as shown in [64], it can be solved by means of a TIQM calculation. Computationally much cheaper, it is possible to show by means of the SQM approach that the proper trend exhibited by the ICSs at such low-energy regime is directly related with a correct description of the asymptotic region employing sufficiently large distances in the calculation [63,88,100,101]. The ICS obtained with the TIQM method is in a nice agreement with both the QCT and SQM predictions.  The rate constants between T = 10 K and 1000 K obtained from the ICSs shown in Figure 9 are presented in Figure 11 for the three initial rovibrational states of the reactant HeH + under consideration. The comparison of the WP values reveals that the ground state (v = 0, j = 0) yields the larger k(T), followed by the one for the first rotationally excited state (v = 0, j = 1). The lowest rate constant in the comparison shown in Figure 11 is obtained for the first vibrationally excited states (v = 1, j = 0). In particular, at T = 100 K, k v=0,j=0 = 6.67 × 10 −10 cm 3 s −1 , whereas k v=0,j=1 = 5.81 × 10 −10 cm 3 s −1 and k v=1,j=0 = 4.49 × 10 −10 cm 3 s −1 . The same sequence is observed among the corresponding SQM rates, but, the QCT calculation, on the contrary predicts certain crossing between the k v=0,j=1 and k v=1,j=0 rates around T ∼ 130 K. WP SQM QCT v=0, j=0 v=0, j=1 v=1, j=0 Ne + HeH + Figure 11. Rate constants (in 10 −1 cm 3 s −1 ) as a function of the temperature (in K) for the Ne + HeH + (v, j) → He + NeH + reaction for v = 0, j = 0 (black), v = 0, j = 1 (red) and v = 1, j = 0 (blue) obtained with the WP (solid lines), SQM (dashed lines) and QCT (dotted lines) calculations of Reference [63]. Logarithmic scales have been used in both axes.
State-to-state rovibrational distributions of the Ne + HeH + (v, j) → He + NeH + (v , j ) reaction for specific values of the collision energy have been also obtained by means of the theoretical approaches discussed in this review. In particular, final state resolved cross sections at two energies, 100 meV and 500 meV, for the reaction initiated from the ground rovibrational HeH + (v = 0, j = 0) state obtained by means of TIQM, SQM, and QCT calculations are compared in Figure 12. Significant discrepancies with the QM result for the production of NeH + in its vibrationaless state v = 0 are observed at E c = 100 meV: On the one hand, the statistical prediction displays a maximum peak for the central rotational states NeH + (v = 0, j ∼ 6), and on the other hand, the QCT distributions peaks at too low final rotational state j = 2 in comparison with the TIQM values.
As the energy increases, on the contrary (see the case of E c = 500 meV at the bottom panel of Figure 12), both SQM and QCT calculations provide a fairly good counterpart of the TIQM rotational distributions for all possible product vibrational states NeH + (v = 0 − 2, j ). The oscillations of the TIQM and QCT rotational distributions as a function of the rotational states for the production of vibrationless NeH + (v = 0) contrast however with the smooth trend followed by the statistical predictions. The comparison of the different theoretical calculations for the DCS also reveals some differences with the energy variations. In Figure 13, results of TIQM, SQM, and two alternatives (GB an HB) for QCT calculations are compared for the same two collision energies considered for the rotational distributions (see Figure 12). The TIQM DCSs exhibit some asymmetry between the forward (θ ∼ 0 degrees) and backward (θ ∼ 180 degrees) scattering directions, a feature which is only described by the QCT calculation at 500 meV (see bottom panel of Figure 13). The imposed-by-construction forward-backward symmetry of the statistical approach fails therefore to give a proper description of the QM result, although both the forward peak at E c = 100 meV and the backward peak at E c = 500 meV.
The two options within the QCT calculation, the GB and HB, lead to similar DCS at both collision energies. The comparison with the TIQM result reveals a nice agreement with the TIQM distribution for the higher energy (500 mev) and noticeable discrepancies for the lower energy (E c = 100 meV).

Ne + NeH + → NeH + + Ne
Also mediated by the presence of a potential well in the intermediate region, [NeHNe] + , of~16.60 kcal mol −1 between reactants and products, the probabilities for the Ne + NeH + → NeH + + Ne reaction display a rich structure of narrow resonances, especially at the low energy regime (see to panels of Figure 14). The SQM predictions remain on a constant value of 0.5, as a consequence of having equal channels for both reactants and products. Despite thresholds for reactions seem to be well described with this approach, the statistical value constitutes a poor average description of the WP reaction probabilities, with the only exception, perhaps of the largest value of the total angular momentum considered in the figure, J = 120.
The comparison with the QM calculations becomes much more favourable for the statistical predictions in the case of the ICS. Figure 15 shows the cross sections for the SQM, TIQM and TDWP approaches in the collision energy between 1 meV and 0.4 eV when the reaction is initiated from the ground rovibrational state NeH + (v = 0, j = 0) The already mentioned difficulties of the WP methods at relatively low energies here manifested in an spurious behaviour of the reaction probability as the energy decreases sufficiently (below ∼0.02 eV). The statistical cross sections, despite an overall qualitatively fair agreement with the TIQM result, remain below up to E c = 0.2 eV.  Ne+NeH + (v=0,j=0) TIQM Figure 15. Cross sections as a function of the collision energy (in eV) for the Ne + NeH + → Ne + NeH + calculated with the TIQM (black line), TDWP (red line), and SQM (blue line) methods in [70].
These quantitative discrepancies seen in the ICSs have also a consequence in the corresponding rate constants. In Figure 16, values of k(T) obtained with the TIQM and SQM cross sections of Equation (15) are compared between 50 K and 900 K. The TIQM rate constants are always larger than the statistical counterparts and, for example, at T = 900 K, k TIQM (T) = 4.46 × 10 −10 cm −3 s −1 , whereas the SQM yields a value of 3.82 × 10 −10 cm −3 s −1 . The constant Langevin value of 4.63 × 10 −10 cm −3 s −1 , also included in Figure 16 for comparison, remain above. Ne + NeH + Figure 16. Rate constants (in 10 −10 cm 3 s −1 ) for the Ne + NeH + reaction calculated with TDWP (black line) and SQM (red line) approaches of Reference [70]. Langevin prediction is also included (dotted line).

Ar + ArH + → ArH + + Ar
This reaction is investigated using quantum and classical dynamical simulations on the analytical PES in [71]. Reaction probabilities for Ar + ArH + (v = 0, j = 0) → ArH + + Ar obtained from TDQM, TIQM and QCT methods are shown in Figure 17 for few selected Js. Numerous resonances are noticeable in the QM probability curves at low collision energies for J = 0. As can be seen in Figures 2 and 3, the potential for Ar 2 H + contains a 15.55 kcal/mol deep well for linear and near collinear Ar-H-Ar geometries which supports the numerous meta-stable collision complexes during the proton exchange reaction. It has been reported in [71] that these intermediate species have average lifetime of ∼0.7 ps determined from the width of the resonances. The probability curves are comparatively smoother for high J values which suggests stripping mechanisms for the reaction.
As can be seen in Figure 17, the agreement between the TIQM and TDQM results is excellent except for the lowest collision energies. These differences originate from the finite time propagation of the WP in the TDQM method. Reaction probabilities calculated using the QCT simulations are also shown in Figure 17 and the results are in a fair agreement with the QM results. For J = 0 and 100, QCT-HB probabilities lie between 0.4 and 0.6, and the curves are rather flat while the QCT-GB probabilities are larger than QCT-HB probabilities in the low-energy region. However, both GB and HB schemes produce similar results for collision energies higher than 0.35 eV for low values of Js. The QCT method largely underestimates the QM probabilities at collision energies ∼0.15-0.4 eV, whereas those are in good agreement beyond 0.4 eV except for J = 200. For J = 300, the QCT-HB method successfully reproduces the overall behavior of the QM probabilities while QCT-GB method slightly overestimates the QM results. The effect of vibrational excitation on the Ar + ArH + → ArH + + Ar reaction was also investigated in [71] using TDQM method. It was found that the vibrational excitation of ArH + greatly reduces the reactivity for this reaction which is similar to the results for Ne + HeH + reaction.
Total ICSs and rate constants for the Ar + ArH + (v = 0, j = 0) → ArH + + Ar reaction computed using TIQM and QCT methods in Reference [71] are given in the left panel of Figure 18. The cross sections are very large at low collision energies and decrease rapidly as E c increases reaching a plateau region. The overall behavior of the cross sections is quite similar to the results obtained for Ne + NeH + . Resonances present in the probability curves are mostly washed out due to the J averaging, and only a few broad peaks can be seen at very low energies. QCT-GB scheme reproduces the QM ICSs quite accurately up to 0.1 eV of collision energy but underestimates the QM results for E > 0.1 eV. The QCT-HB cross sections have a constant slope in the log-log plot which suggests a simple capture type cross sections typical for ion-molecule reactions. Initial state selected rate constants for the proton exchange reaction between two Ar atoms with the reactant ArH + in its ground ro-vibrational state are shown in the right panel of Figure 18. The QM rate constant increases with increase in temperature and becomes almost constant with values ∼6×10 −10 cm 3 s −1 at T > 600 K. The QCT-GB results follow a similar trend whereas QCT-HB rate constants remain almost constant at ∼6×10 −10 cm 3 s −1 in the entire temperature range. The rate constant calculated following the Langevin's capture model [96,97] for ion-molecule reaction is ∼6.88×10 −10 cm 3 s −1 which is slightly higher than those obtained from dynamical simulations. However, it was found from QCT trajectories that only a fraction of trajectories which form collision complexes leads to products via proton exchange reaction and others produce reactants again.

Conclusions
Reactions in protonated rare gas (Rg) systems have been intensely investigated in recent years due to their astrophysical interest. In this review, we present results of previous studies on the reactions He + HeH + → HeH + + He; Ne + NeH + → NeH + + Ne; Ar + ArH + → ArH + + Ar; He + NeH + → HeH + + Ne and Ne + HeH + → NeH + + He. Different observables such as reaction probabilities, cross sections, and rate constants for those reactions have been obtained by means of exact quantum and quasi-classical trajectory (QCT) calculations on high-level ab initio energy-based analytical potential energy surfaces (PESs). The existence of a relatively deep well in the intermediate region between reactants and products enable as well the use of a statistical quantum method (SQM) which provides a fairly good average description of the overall dynamics. Different approximations introduced to both the time-dependent and time-independent quantum and statistical techniques such as the centrifugal sudden (CS) option have been successfully tested and alternatives for the binning procedure in the QCT such as the histogram and Gaussian approaches were comparatively employed. In general, all methods are found to be adequate for the study of the title processes although the CS approximation can introduce significant discrepancies with respect to the more rigorous coupled channel scheme. Moreover, the description provided by the time dependent techniques can be affected by some deficiencies of the wave packet propagation at the low energy regime. We think that one of the challenges in this area is the comparative analysis of these reactions when heavier Rg atoms are considered. It would be interesting, for instance, to see whether or not in these cases, the complex-forming mechanisms play a major role on the overall dynamics of the process, thus justifying the use of statistical techniques such as the SQM discussed here. Finally, we conclude that this complete sort of investigations which include as a fundamental first step the development of full ab initio PES represent a convenient strategy to understand the dynamics of the title reactions.

Abbreviations
The following abbreviations are used in this manuscript: