Theoretical Study of the Thermal Rate Coefficients of the H3+ + C2H4 Reaction: Dynamics Study on a Full-Dimensional Potential Energy Surface

Ion–molecular reactions play a significant role in molecular evolution within the interstellar medium. In this study, the entrance channel reaction, H3+ + C2H4 → H2 + C2H5+, was investigated using classical molecular dynamic (classical MD) and ring polymer molecular dynamic (RPMD) simulation techniques. We developed an analytical potential energy surface function with a permutationally invariant polynomial basis, specifically employing the monomial symmetrized approach. Our dynamic simulations reproduced the rate coefficient of 300 K for H3+ + C2H4 → H2 + C2H5+, aligning reasonably well with the values in the kinetic database commonly utilized in astrochemistry. The thermal rate coefficients obtained using both the classical MD and RPMD techniques exhibited an increase from 100 K to 300 K as the temperature rose. Additionally, we analyzed the excess energy distribution of the C2H5+ fragment with respect to temperature to investigate the indirect reaction pathway of C2H5+ → H2 + C2H3+. This result suggests that the indirect reaction pathway of C2H5+ → H2 + C2H3+ holds minor significance, although the distribution highly depends on the collisional temperature.


Introduction
The H 3 + molecule plays a crucial role in interstellar chemistry due to its efficiency as a proton donor in H 3 + + X → H 2 + HX + ion-molecule reactions, where X represents a neutral molecule in the interstellar medium [1][2][3].Many interstellar molecules detected using spectroscopic techniques fulfill this proton affinity condition.Furthermore, these ion-molecule reactions occur efficiently due to the absence of an entrance reaction barrier, primarily driven by strong attractive forces like charge-dipole or charge-induced-dipole interactions.Consequently, H 3 + had often been referred to as the initiator of interstellar chemistry in previous studies [1,[3][4][5][6].Ethylene (C 2 H 4 ) also plays a role in interstellar chemistry since the molecule has been detected in the circumstellar envelope of the carbon-rich asymptotic giant branch star IRC+10216 [7][8][9][10].Therefore, laboratory studies have been carried out to determine the rate coefficient for H 3 + + C 2 H 4 reactions [11][12][13][14].H 3 + + C 2 H 4 reactions predominantly yield the H 2 + C 2 H 5 + product via a proton-donating process with considerable exothermicity [13,14].Although this barrierless reaction appears to be straightforward, experimental investigations have yet to address the temperature dependence of the reaction rate coefficients [11][12][13][14].While Langevin's theory is commonly utilized for estimating the low-temperature rate coefficients in various ion-molecule reactions, calculating the reaction rate coefficients using first principles may be essential for assessing the validity of Langevin's theory.From an astrochemical perspective, understanding the H + can lead to the formation of three species due to its exothermic nature.This reaction follows an indirect mechanism, where the C 2 H 3 + product is formed from vibrationally excited C 2 H 5 + via a unimolecular reaction, with vibrational energy provided in the initial proton transfer process from H 3 + .Thus, the process can be formally expressed as H + fragments.Experimental measurements using ion cyclotron resonance under room temperature and low-pressure conditions have yielded branching fraction values of 0.40:0.60[13] and 0.31:0.69[14] for [C 2 H 5  + ]:[C 2 H 3  + ], indicating that the C 2 H 3 + ion is the predominant reaction product.It is worth mentioning that previous experimental studies have reported the lifetime of the primary C 2 H 5 + reaction product to be in the range from 10 −8 to 10 −6 s [14].This extended lifetime suggests that the indirect mechanism prevails in the C 2 H 3 + production channel, which is consistent with previous quantum chemistry calculations conducted by Uggerud, who reported barrier and exit energy values of 55.9 and 50.0 kcal/mol (234 and 209 kJ/mol), respectively, at the MP2/6-31G(d,p) level of theory [15].The overall rate coefficient at 300 K published in the KIDA (Kinetic Database for Astrochemistry) database is 2.9 × 10 −9 cm 3 s −1 [16], with a branching fraction of 0.30:0.70 for [C 2 H 5  + ]:[C 2 H 3 + ], as reported in the experiments mentioned earlier and conducted using a drift chamber mass spectrometer [17].Based on the branching fraction results, the C 2 H 5 + * fragment acquires significant rovibrational energies from the H + system were studied using quantum chemical calculations [18][19][20].These studies identified several equilibrium structures of the C 2 H 7 + potential energies, with the proton-bridged (CH 3 -H-CH 3 ) + structure and its associated forms being the most stable [18][19][20].While the dissociation channels of the C 2 H reaction using an ab initio quantum chemical calculation dataset obtained from on-thefly preliminary dynamics simulation results.Monomial symmetrized approach (MSA2) software developed by Bowman and coworkers [21,22] was employed to fit the potential energies and their gradients on a permutational invariant polynomial basis.Details regarding the development of the potential energy surface (PES) are provided in the Methodology section.Subsequently, we present the computational outcomes of the rate coefficients for the H 3 + + C 2 H 4 → H 2 + C 2 H 5 + reaction on the developed PES employing classical molecular dynamic and ring polymer molecular dynamic (RPMD) methods [23][24][25][26], which represent the quantum mechanical behavior of nuclei by considering them as cyclic beads.As mentioned above, there is a substantial barrier in the C 2 H 5 + → H 2 + C 2 H 3 + reaction [15], and the C 2 H 5 + fragment has a relatively long lifetime [14].Therefore, we evaluate the three-body dissociation process using the internal energy of C 2 H 5 + *.S1 in the Supplementary Materials.Although the H 2 bond direction relative to the CC bond differed from the MP2 results, the other geometric parameters were quite similar.In a barrierless ion-molecule reaction involving large exothermicity, the molecule has relatively higher rovibrational energies or a large kinetic energy release (KER).Therefore, the quality of the PES fit in this area is considered to have little effect on these results.Additionally, the MP2 results were compared with high-quality CCSD(T) results, showing consistent findings across all the datasets.While our MP2 results for the C 2 H 5 + → H 2 + C 2 H 3 + reaction showed a slightly higher barrier compared to that of the previous MP2 calculations, they exhibited good correlation with the CCSD(T) results, as shown in Table 1.In a similar reaction, C 2 H 7 + → H 2 + C 2 H 5 + , the barrier height of the CCSD(T) was slightly higher than that of the DFT calculation [18].It should be noted that previous studies indicated that the transition state barrier tends to be higher with MP2 than with CCSD(T) [27][28][29].

Results and Discussion
Molecules 2024, 29, x FOR PEER REVIEW 3 of 13

Potential Energy Profile
Figure 1a,b display the potential energy profiles during the H3 + + C2H4 → H2 + C2H5 + and C2H5 + → H2 + C2H3 + reactions, respectively.In the H3 + + C2H4 reaction, the reactants directly formed the H2 + C2H5 + products through the C2H5 + •••H2 van der Waals well.Conversely, in the generation of H2 + C2H3 + products, there was a potential energy barrier between C2H5 + and the resulting H2 + C2H3 + fragments.Table 1 presents a comparison of the potential energies at the stationary points obtained from the analytical potential energy function and the MP2/cc-pVDZ data for the H3 + + C2H4 → H2 + C2H5 + reaction.As mentioned above, this sampling scheme does not cover the PES region for the C2H5 + → H2 + C2H3 + dissociation process, which will be discussed separately in terms of the vibrational energy distributions of C2H5 + generated from the H3 + + C2H4 → H2 + C2H5 + reaction.However, the C2H5 + •••H2 van del Waals well of the PES is relatively deeper than that in the MP2 results; our constructed PES yielded results that reasonably matched those from the MP2 calculations.The molecular structure of the C2H5 + •••H2 intermediate complex and key geometries for PES and MP2 are compared in Figure S1 in the Supplementary Materials.Although the H2 bond direction relative to the CC bond differed from the MP2 results, the other geometric parameters were quite similar.In a barrierless ion-molecule reaction involving large exothermicity, the molecule has relatively higher rovibrational energies or a large kinetic energy release (KER).Therefore, the quality of the PES fit in this area is considered to have little effect on these results.Additionally, the MP2 results were compared with high-quality CCSD(T) results, showing consistent findings across all the datasets.While our MP2 results for the C2H5 + → H2 + C2H3 + reaction showed a slightly higher barrier compared to that of the previous MP2 calculations, they exhibited good correlation with the CCSD(T) results, as shown in Table 1.In a similar reaction, C2H7 + → H2 + C2H5 + , the barrier height of the CCSD(T) was slightly higher than that of the DFT calculation [18].It should be noted that previous studies indicated that the transition state barrier tends to be higher with MP2 than with CCSD(T) [27][28][29].Table 1.Relative energies for reactions (a) and (b) for MSA2 PES, MP2/cc-pVDZ, DF-CCSD(T)-F12a/cc-pVDZ, and cc-pVTZ, as well as MP2/6-31G(d,p).The zero energies are defined by the potential energy of the reactant.

Rate Coefficients
The collision simulations were conducted using both the classical MD and RPMD methods over a temperature range from 100 to 300 K.The vibrational and rotational energies for the H 3 + + C 2 H 4 reactants, along with their relative translational energy, were thermalized at each temperature in the simulations.In this scenario, the thermal rate coefficient k(T) was obtained using the following equation [30]: where k B , µ, and b max represent the Boltzmann constant, the reduced mass, and the maximum impact parameter, respectively.The standard error ∆k(T) is given as follows: where N r and N t denote the number of reacted and total trajectories.Further details about the collision procedure are provided in the Methodology section below.1b).Consequently, it is anticipated that the H 2 + C 2 H 3 + products will form via an indirect process involving the vibrational excited C 2 H 5 + fragment.
Table 2. Thermal rate coefficients with the standard errors (k(T) ± ∆k(T)) estimated from classical MD and RPMD results, along with the impact parameter (b max ) and the number of reacted (N r ) and total trajectories (N t ), as well as the coefficients in the KIDA database.+ products that were indirectly generated.At 300 K, both sets of our results are somewhat better than the experimental data (green scatter in Figure 2) provided in the KIDA database, but we reasonably reproduced the experimental results.Both the classical MD and RPMD rate coefficients decrease as the temperature decreases.This temperature dependency is similar to the coefficients observed in other ion-neutral reactions, such as H − + C 2 H 2 → H 2 + C 2 H − and H 3 + + CO → H 2 + HCO + /HOC + reactions, as derived from the classical and RPMD simulations [30,31].Notably, this temperature dependency cannot be captured using the temperature-independent Langevin rate equation, which is solely based on the chargeinduced-dipole interaction.The rate coefficients obtained using RPMDs were slightly elevated compared to those obtained using classical MDs, suggesting that the quantum fluctuations in nuclei contribute to their increase.In RPMDs, ion-neutral attraction extends farther than the attraction between classical particles in classical MDs, owing to fluctuations in the nuclear probability densities.This observation aligns closely with findings from the [30].In the following section, we will discuss the energy

Internal Energies of Fragments
To estimate the rovibrational energy distributions of the C2H5 + fragment concerning the indirect reaction C2H5 + → H2 + C2H3 + , we examined the system internal energies (ε) of the C2H5 + fragment derived using the following equation: where  C H represents the vibrational and rotational energies for centroid velocities of the C2H5 + fragment, and  H + C H denotes the potential energy obtained from the MSA2 PES concerning the nuclear configuration of j-th bead.Additionally,  is a sixthorder polynomial function dependent on the internuclear distance (r) between the H2 fragments, where the center-of-mass distance (R) between the H2 and C2H5 + fragments is fixed at 20 Å. Figure 3a,b illustrate the relaxed potential energy curves as a function of R and r and RPMD outcomes (red line).The green scatter plot denotes the rate coefficients extracted from the KIDA database for the overall reactions.

Internal Energies of Fragments
To estimate the rovibrational energy distributions of the C 2 H 5 + fragment concerning the indirect reaction C 2 H 5 + → H 2 + C 2 H 3 + , we examined the system internal energies (ε) of the C 2 H 5 + fragment derived using the following equation: where E k C 2 H + 5 represents the vibrational and rotational energies for centroid velocities of the C 2 H 5 + fragment, and V j H 2 + C 2 H + 5 denotes the potential energy obtained from the MSA2 PES concerning the nuclear configuration of j-th bead.Additionally, V H 2 is a sixth-order polynomial function dependent on the internuclear distance (r) between the H 2 fragments, where the center-of-mass distance (R) between the H 2 and C 2 H 5  + fragments is fixed at 20 Å. Figure 3a,b illustrate the relaxed potential energy curves as a function of R and r, respectively.Notably, no interaction was observed between the H 2 and C 2 H 5 + fragments when R exceeded 13 Å, as depicted in Figure 3a.Additionally, Figure 3b confirms that the sixth-order polynomial function effectively reproduced our PES for the C 2 H 7 + system.For the final step, the trajectory data for R exceeding 13 Å are represented in the right-hand side of Equation (3).
Figure 4 displays the ε of the C 2 H 5 + fragment at temperatures of T = 100, 200, and 300 K. Notably, the peaks of ε using the classical MD and RPMD techniques at all the temperatures were below the asymptotic region (60.5 kcal/mol at the MP2 level) for the H 2 + C 2 H 3 + products, suggesting that the C 2 H 5 + → H 2 + C 2 H 3 + dissociation channel was not dominant, although the distribution at T = 300 K indicates that the C 2 H 5 + fragment with larger internal energies can lead to the H 2 + C 2 H 3 + dissociation channel.Both ε values decrease using the classical MD and RPMD techniques as the temperature declines, suggesting that the temperature-dependent behaviors of the internal energies of the H 3 + and C 2 H 4 reactants, along with their collision energy, influence the proton abstraction process of C 2 H 4 from H 3 + .Consequently, we infer that the branching fraction for [C 2 H 5  + ]:[C 2 H 3  + ] is strongly temperature-dependent, potentially leading to the overestimation of the branching fraction of the C 2 H 3 + product in previous experiments [13,14,17].It is worth noting that the MSA2 PES was developed using the MP2 results for the dynamics simulations, and the C 2 H 5 + •••H 2 van del Waals well of the PES was relatively deeper than that in the MP2 results, as mentioned above.Further elaboration, along with quantitative results, can be expected by conducting theoretical dynamic simulations using a high-quality PES constructed with sophisticated techniques, such as the ∆-machine learning algorithm [32][33][34].Moreover, as the temperature decreases, the RPMD distribution becomes narrower, indicating that at low temperatures, the quantum nuclei exhibit characteristic quantum behavior in the internal energy of C 2 H 5 + .It should be noted that the energy distributions were derived from a snapshot at the final step rather than from the Fourier transformation of the auto-correlation function.
FOR PEER REVIEW 6 of 13 respectively.Notably, no interaction was observed between the H2 and C2H5 + fragments when R exceeded 13 Å, as depicted in Figure 3a.Additionally, Figure 3b confirms that the sixth-order polynomial function effectively reproduced our PES for the C2H7 + system.For the final step, the trajectory data for R exceeding 13 Å are represented in the right-hand side of Equation ( 3). Figure 4 displays the ε of the C2H5 + fragment at temperatures of T = 100, 200, and 300 K. Notably, the peaks of ε using the classical MD and RPMD techniques at all the temperatures were below the asymptotic region (60.5 kcal/mol at the MP2 level) for the H2 + C2H3 + products, suggesting that the C2H5 + → H2 + C2H3 + dissociation channel was not dominant, although the distribution at T = 300 K indicates that the C2H5 + fragment with larger internal energies can lead to the H2 + C2H3 + dissociation channel.Both ε values decrease using the classical MD and RPMD techniques as the temperature declines, suggesting that the temperature-dependent behaviors of the internal energies of the H3 + and C2H4 reactants, along with their collision energy, influence the proton abstraction process of C2H4 from H3 + .Consequently, we infer that the branching fraction for [C2H5 + ]:[C2H3 + ] is strongly temperaturedependent, potentially leading to the overestimation of the branching fraction of the C2H3 + product in previous experiments [13,14,17].It is worth noting that the MSA2 PES was developed using the MP2 results for the dynamics simulations, and the C2H5 + •••H2 van del Waals well of the PES was relatively deeper than that in the MP2 results, as mentioned above.Further elaboration, along with quantitative results, can be expected by conducting theoretical dynamic simulations using a high-quality PES constructed with sophisticated techniques, such as the Δ-machine learning algorithm [32][33][34].Moreover, as the temperature decreases, the RPMD distribution becomes narrower, indicating that at low temperatures, the quantum nuclei exhibit characteristic quantum behavior in the internal energy of C2H5 + .It should be noted that the energy distributions were derived from a snapshot at the final step rather than from the Fourier transformation of the auto-correlation function.Figure 5 illustrates the distributions of vibrational (ε vib ) and rovibrational states (ε rovib ) for the H 2 fragment in the H 3 + + C 2 H 4 → H 2 + C 2 H 5 + reaction.The distributions were obtained using the following equations: and where E vib k (H 2 ) and E rot k (H 2 ) represent the vibrational and rotational energies for the centroid velocities of the H 2 molecule.Unlike the distribution of the system's internal energy for the C 2 H 5 + fragment, the internal energy of H 2 , having fulfilled its role as a proton donor, is minimally affected by the temperature.As the temperature rises, the distributions for both ε vib and ε rovib expand due to the incorporation of thermal energies.The internal energies of classical MDs are notably low, primarily because the classical scheme does not account for zero-point energy.The peaks of the vibrational state distributions for RPMDs on the potential energy curve, accounting for anharmonicity, are approximately 4 kcal/mol.Notice that the zero-point energy of H 2 , which was calculated using the harmonic analysis of the MP2/cc-pVDZ level, is 6.4 kcal/mol.Figure 5 illustrates the distributions of vibrational (εvib) and rovibrational states (εr for the H2 fragment in the H3 + + C2H4 → H2 + C2H5 + reaction.The distributions were tained using the following equations: and where  H and  H represent the vibrational and rotational energies for centroid velocities of the H2 molecule.Unlike the distribution of the system's internal ergy for the C2H5 + fragment, the internal energy of H2, having fulfilled its role as a pro donor, is minimally affected by the temperature.As the temperature rises, the distri tions for both εvib and εrovib expand due to the incorporation of thermal energies.The ternal energies of classical MDs are notably low, primarily because the classical sche The kinetic energy release (KER), which was not distributed along the internal modes, is depicted in Figure 6.These distributions were obtained using the following equations: where µ rel and → v rel denote the reduced mass and relative velocity vector between the H 2 and C 2 H 5 + fragments.The distributions broadened, reflecting the distribution of the internal energies of the H 2 and C 2 H 5 + fragments.However, no significant temperature dependence was observed in the relative energy for both the classical MD and RPMD results.It should be noted that the KER might have been underestimated because the exit channel for the H 2 + C 2 H 5 + products of our PES was relatively higher than the MP2 results.
Molecules 2024, 29, x FOR PEER REVIEW 8 of 13 4 kcal/mol.Notice that the zero-point energy of H2, which was calculated using the harmonic analysis of the MP2/cc-pVDZ level, is 6.4 kcal/mol.The kinetic energy release (KER), which was not distributed along the internal modes, is depicted in Figure 6.These distributions were obtained using the following equations: where  and  ⃗ denote the reduced mass and relative velocity vector between the H2 and C2H5 + fragments.The distributions broadened, reflecting the distribution of the internal energies of the H2 and C2H5 + fragments.However, no significant temperature dependence was observed in the relative energy for both the classical MD and RPMD results.It should be noted that the KER might have been underestimated because the exit channel for the H2 + C2H5 + products of our PES was relatively higher than the MP2 results.

Development of a Global Potential Energy Surface
All the quantum chemistry computations to construct the full-dimensional C2H7 + PES were conducted at the MP2/cc-pVDZ level using Gaussian09 [35].This level of calculation was chosen due to its ability to provide a substantial amount of energy and gradient data across various structures within a reasonable computational timeframe.This study uses a barrierless proton transfer reaction in the electronic singlet ground state.This proton transfer reaction, not involving electron transfer, has a relatively small electron correlation energy value.Table 1 shows that the MP2 energy is quite similar to the CCSD(T) results.Our previous dynamics study for the ion-molecule reaction, such as H − + C2H2 → H2 + C2H − , qualitatively reproduced the rate coefficient [30].It should be noted that for the NH3 + + H2 → NH3 + + H reaction, which involves a substantial barrier to hydride transfer, the CCSD(T) level was employed [36].
In the barrierless ion-molecule reactions, strong attractive forces like charge-dipole or charge-induced-dipole interactions play a crucial role.Therefore, quasi-classical trajectory calculations were computed for both the H3 + + C2H4 reactants and the H2 + C2H5 + products to sample the structural data points.A total of 220 trajectories from the reactant side  PES were conducted at the MP2/cc-pVDZ level using Gaussian09 [35].This level of calculation was chosen due to its ability to provide a substantial amount of energy and gradient data across various structures within a reasonable computational timeframe.This study uses a barrierless proton transfer reaction in the electronic singlet ground state.This proton transfer reaction, not involving electron transfer, has a relatively small electron correlation energy value.Table 1 shows that the MP2 energy is quite similar to the CCSD(T) results.Our previous dynamics study for the ion-molecule reaction, such as , qualitatively reproduced the rate coefficient [30].It should be noted that for the NH 3 + + H 2 → NH 3 + + H reaction, which involves a substantial barrier to hydride transfer, the CCSD(T) level was employed [36].
In the barrierless ion-molecule reactions, strong attractive forces like charge-dipole or charge-induced-dipole interactions play a crucial role.Therefore, quasi-classical trajectory calculations were computed for both the H 3 + + C 2 H 4 reactants and the H 2 + C 2 H 5 + products to sample the structural data points.A total of 220 trajectories from the reactant side and 60 trajectories from the product side with collision energies of 5, 10, and 20 kcal/mol were run, yielding 465,000 data points.All the trajectories from the reactants reached the H 2 + C 2 H 5 + products, without directly producing the 2H 2 + C 2 H 3 + product.Consequently, data sampling of the H 2 + H 2 + C 2 H 3 + fragments was not performed.These data points were fitted to an analytical function comprising permutationally invariant polynomial basis sets using the MSA2 code developed by Bowman and coworkers [21,22].The preliminary fit helped identify the unphysical leaky holes, which frequently occur in PES regions with insufficient data points [37].An additional 2790 data points were added to fill these holes.The structures with higher energies, specifically those up to 0.2 hartrees from the C 2 H 5 + ••• H 2 complex, were subsequently excluded, reducing the number of data points to approximately 330,000.A fourth-order polynomial was employed, resulting in a final root-mean-square error of 959 cm −1 for the fit.Figure S2 in the Supplementary Materials displays the plotted fitted potential energies relative to the MP2/cc-pVDZ energies.The fitted potential energy surface (PES) reasonably replicated the reaction energies for both the reactants and products, including the C 2 H 5 + ••• H 2 intermediate complex, along with the vibrational frequencies for this intermediate complex (see Figure 1 and Table S1).For comparison with the MP2/cc-pVDZ results, the density fitting explicitly correlated coupledcluster singles and doubles plus perturbative triples (DF-CCSD(T)-F12) method [38,39] with cc-pVDZ and cc-VTZ basis sets implemented using Molpro 2023.2 was utilized [40].Global Reaction Route Mapping (GRRM) calculations [41][42][43] were performed to obtain the equilibrium and transition state structures.

Procedure for Molecular Dynamics
All the nuclear dynamics simulations conducted in this study were based on the fitted PES.Initially, path integral molecular dynamic (PIMD) simulations were executed to derive the initial coordinates and momenta for collisional simulations spanning temperatures from T = 100 K to 300 K.The ring polymer Hamiltonian H(p, r) [30] is defined as follows: where ℏ and β represent the reduced Planck constant and reciprocal temperature, respectively, where β ≡ 1/k B T. m i stands for the atomic mass of the i-th nucleus, and V signifies the potential energy of the system.N bead represents the number of beads, while r j i and p j i denote the Cartesian coordinates and their conjugated momentum vectors of the j-th bead for the i-th atom, respectively.Utilizing the ring polymer Hamiltonian [30] defined by the PIMD method, phase information was obtained in accordance with the quantum Boltzmann distribution, effectively capturing the nuclear quantum effects, including discretized vibrational energies with an appropriate number of beads (N bead ).The convergence of N bead is demonstrated in Figure S3 in the Supplementary Materials, illustrating its association with the system's internal energy.In this study, RPMD simulations were conducted with varying numbers of beads, N bead = 96, 64, 48, 48, and 24 for temperatures T = 100, 150, 200, 250, and 300 K, respectively.To maintain the configurations where the C 2 H 4 and H 3 + reactants do not interact, a harmonic bias potential was employed alongside a Nosé-Hoover thermostat for canonical ensembles.The integration of the equations of the motion of the ring polymer Hamiltonian was carried out using the velocity-Verlet method with a time step of ∆t = 0.10 fs, totaling 10 4 -10 6 simulation steps.Subsequently, RPMD simulations were performed, extending the PIMD method to enable real-time dynamics simulations, which are particularly adept at capturing nuclear quantum effects, such as zero-point energy and tunneling [30,36,[44][45][46][47][48][49].The impact parameter (b) for collisional simulations was set below b = b max ζ 1/2 , where b max represents the maximum impact parameter, and ζ is a random number within the range of [0, 1].Following this, the reactant coordinate was adjusted based on the specified impact parameter.Further specifics regarding the collisional simulation can be found in our previous publication [30].Eventually, the RPMD trajectories were propagated by solving the equations of the motion of the ring polymer Hamiltonian without a thermostat, employing a time step of ∆t = 0.10 fs.The total simulation time spanned 1.5-22 ps.Additionally, classical MD simulations were conducted for comparison with the RPMD results, following a similar procedure except for the number of beads, which remained constant at one for all the temperatures in the classical MD simulations.All the calculations for PIMDs, RPMDs, and classical MDs were performed using the open-source code PIMD.ver.2.6.0 [50].

Figure 1 .
Figure 1.Diagram of reaction coordinates, illustrating the geometries at stationary points for (a) the H3 + + C2H4 → H2 + C2H5 + reaction obtained using MSA2 PES and (b) the C2H5 + → H2 + C2H3 + reaction calculated at the MP2/cc-pVDZ level.Zero energy is defined as the H3⁺ + C2H4 reactant energy level in (a), while it is defined as the C2H5⁺ energy level in (b).

Figure 1 .
Figure 1.Diagram of reaction coordinates, illustrating the geometries at stationary points for (a) the H 3 + + C 2 H 4 → H 2 + C 2 H 5 + reaction obtained using MSA2 PES and (b) the C 2 H 5 + → H 2 + C 2 H 3 + reaction calculated at the MP2/cc-pVDZ level.Zero energy is defined as the H 3 + + C 2 H 4 reactant energy level in (a), while it is defined as the C 2 H 5 + energy level in (b).

Figure 2 .
Figure 2. The thermal rate coefficients for the overall reactions, encompassing both H3 + + C2H4 → H + C2H5 + and H3 + + C2H4 → 2H2 + C2H3 + , derived from the classical MD results (blue line) and RPMD outcomes (red line).The green scatter plot denotes the rate coefficients extracted from the KIDA database for the overall reactions.

Figure 3 .
Figure 3.The potential energy curves (red lines) of the C2H7 + system as a function of (a) the centerof-mass distance (R) between the H2 and C2H5 + fragments and (b) the internuclear distance (r) between the H2 fragments, respectively.The molecular structures of R and r are depicted as insets in Figure3a.The blue lines in Figure3brepresent the 6th order polynomial function dependent on r.The zero energies are defined at the potential energy of the asymptote region between H2 and C2H5 + .

Figure 3 .
Figure 3.The potential energy curves (red lines) of the C 2 H 7 + system as a function of (a) the center-ofmass distance (R) between the H 2 and C 2 H 5 + fragments and (b) the internuclear distance (r) between the H 2 fragments, respectively.The molecular structures of R and r are depicted as insets in Figure 3a.The blue lines in Figure 3b represent the 6th order polynomial function dependent on r.The zero energies are defined at the potential energy of the asymptote region between H 2 and C 2 H 5 + .

Figure 5 .
Figure 5.The vibrational (εvib) and the rovibrational states (εrovib) of the H2 fragment for (a) classica MDs at T = 100 K, (b) RPMDs at T = 100 K, (c) classical MDs at T = 200 K, (d) RPMDs at T = 200 K (e) classical MDs at T = 300 K, and (f) RPMDs at T = 300 K.The purple and blue plots depict εvib and εrovib for classical MDs, whereas the magenta and red illustrate εvib and εrovib for RPMDs.

Figure 6 .
Figure 6.Kinetic energy release E (KER) corresponding to the relative translational energy between H 2 and C 2 H 5 + fragments for (a) classical MDs at T = 100 K, (b) RPMDs at T = 100 K, (c) classical MDs at T = 200 K, (d) RPMDs at T = 200 K, (e) classical MDs at T = 300 K, and (f) RPMDs at T = 300 K.

3 . Methodology 3 . 1 .
Development of a Global Potential Energy SurfaceAll the quantum chemistry computations to construct the full-dimensional C 2 H 7 + 3 + + C 2 H 4 → H 2 + C 2 H 5 + reaction rate coefficients and their temperature dependence is quite important.As mentioned above, this ion-molecule reaction yields H 2 + C 2 H 5 + predominantly.It is worth noting that the reaction H 3 + + C 2 H 4 → 2H 2 + C 2 H 3 •••H 2 van del Waals well of the PES is relatively deeper than that in the MP2 results; our constructed PES yielded results that reasonably matched those from the MP2 calculations.The molecular structure of the C 2 H 5 + •••H 2 intermediate complex and key geometries for PES and MP2 are compared in Figure 2.1.Potential Energy Profile Figure 1a,b display the potential energy profiles during the H 3 + + C 2 H 4 → H 2 + C 2 H 5 Table 2 provides the essential details required for computing the rate coefficients.The parameter N r presented in Table 2 denotes the count of the trajectories resulting in H 2 + C 2 H 5 + products.Notably, in our simulations, no trajectories led to 2H 2 + C 2 H 3 + products, which is consistent with the anticipated high potential energy barrier in the C 2 H 5 + → H 2 + C 2 H 3 + reaction (refer to Figure Figure2illustrates the thermal rate coefficients across temperatures ranging from 100 to 300 K. Our results from both the classical MD and RPMD simulations can be interpreted as the rate coefficients for the overall reaction, encompassing both the H 3 + + C 2 H 4 → H 2 + C 2 H 5