Theoretical Study of the Thermal Rate Coefficients of the H₃⁺ + C₂H₄ Reaction: Dynamics Study on a Full-Dimensional Potential Energy Surface

Abstract

Ion–molecular reactions play a significant role in molecular evolution within the interstellar medium. In this study, the entrance channel reaction, H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺, 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 H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺, 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 C₂H₅⁺ fragment with respect to temperature to investigate the indirect reaction pathway of C₂H₅⁺ → H₂ + C₂H₃⁺. This result suggests that the indirect reaction pathway of C₂H₅⁺ → H₂ + C₂H₃⁺ holds minor significance, although the distribution highly depends on the collisional temperature.

1. Introduction

The H₃⁺ molecule plays a crucial role in interstellar chemistry due to its efficiency as a proton donor in H₃⁺ + X → H₂ + 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₃⁺ had often been referred to as the initiator of interstellar chemistry in previous studies [1,3,4,5,6]. Ethylene (C₂H₄) 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₃⁺ + C₂H₄ reactions [11,12,13,14]. H₃⁺ + C₂H₄ reactions predominantly yield the H₂ + C₂H₅⁺ 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₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction rate coefficients and their temperature dependence is quite important.

As mentioned above, this ion–molecule reaction yields H₂ + C₂H₅⁺ predominantly. It is worth noting that the reaction H₃⁺ + C₂H₄ → 2H₂ + C₂H₃⁺ can lead to the formation of three species due to its exothermic nature. This reaction follows an indirect mechanism, where the C₂H₃⁺ product is formed from vibrationally excited C₂H₅⁺ via a unimolecular reaction, with vibrational energy provided in the initial proton transfer process from H₃⁺. Thus, the process can be formally expressed as H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺* → H₂ + H₂ + C₂H₃⁺, where C₂H₅⁺* represents a molecule with sufficient vibrational energy to produce the H₂ + C₂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₂H₅⁺]:[C₂H₃⁺], indicating that the C₂H₃⁺ ion is the predominant reaction product. It is worth mentioning that previous experimental studies have reported the lifetime of the primary C₂H₅⁺ reaction product to be in the range from 10⁻⁸ to 10⁻⁶ s [14]. This extended lifetime suggests that the indirect mechanism prevails in the C₂H₃⁺ 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⁻⁹ cm³ s⁻¹ [16], with a branching fraction of 0.30:0.70 for [C₂H₅⁺]:[C₂H₃⁺], as reported in the experiments mentioned earlier and conducted using a drift chamber mass spectrometer [17]. Based on the branching fraction results, the C₂H₅⁺* fragment acquires significant rovibrational energies from the H₃⁺ + C₂H₄ proton-donating process, leading to its dissociation into H₂ + C₂H₃⁺ fragments. Understanding the internal energy of C₂H₅⁺ in the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction dynamics and its temperature dependence is essential for elucidating the branching mechanism of C₂H₅⁺/H₂ + C₂H₃⁺. Comprehending the branching mechanism as well as the thermal rate coefficient is critical in astrochemistry.

Motivated by this goal, we investigated the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction dynamics. Previously, the potential energy profiles of the entire C₂H₇⁺ system were studied using quantum chemical calculations [18,19,20]. These studies identified several equilibrium structures of the C₂H₇⁺ potential energies, with the proton-bridged (CH₃-H-CH₃)⁺ structure and its associated forms being the most stable [18,19,20]. While the dissociation channels of the C₂H₇⁺ system, such as H₂ + C₂H₅⁺, CH₃⁺ + CH₄ and H + C₂H₆⁺, were considered, the reaction profile for H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ was not examined. Therefore, we initially constructed a full-dimensional potential energy surface for the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction using an ab initio quantum chemical calculation dataset obtained from on-the-fly 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₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ 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₂H₅⁺ → H₂ + C₂H₃⁺ reaction [15], and the C₂H₅⁺ fragment has a relatively long lifetime [14]. Therefore, we evaluate the three-body dissociation process using the internal energy of C₂H₅⁺*.

2. Results and Discussion

2.1. Potential Energy Profile

Figure 1a,b display the potential energy profiles during the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ and C₂H₅⁺ → H₂ + C₂H₃⁺ reactions, respectively. In the H₃⁺ + C₂H₄ reaction, the reactants directly formed the H₂ + C₂H₅⁺ products through the C₂H₅⁺···H₂ van der Waals well. Conversely, in the generation of H₂ + C₂H₃⁺ products, there was a potential energy barrier between C₂H₅⁺ and the resulting H₂ + C₂H₃⁺ fragments.

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.

	PES	MP2	DF-CCSD(T)-F12a		MP2
		cc-pVDZ	cc-pVDZ	cc-pVTZ	6-31G(d,p) a
(a) H3+ + C2H4 → [C2H5+···H2] → H2 + C2H5+
Reactant	0.0	0.0	0.0	0.0	—
Intermediate complex	−72.9	−65.2	−65.3	−63.8	—
Product	−58.5	−63.7	−63.6	−62.0	—
(b) C2H5+ → [C2H3+···H2]‡ → [C2H3+···H2] → H2 + C2H3+
Reactant	—	−63.7 (0.0)	−63.6 (0.0)	−62.0 (0.0)	(0.0)
Transition state	—	2.3 (66.0)	−0.2 (63.4)	−0.6 (61.4)	(55.9)
Intermediate complex	—	−5.7 (58.0)	−4.2 (59.4)	−5.3 (56.7)	—
Product	—	−3.2 (60.5)	−1.5 (62.1)	−2.6 (59.4)	(50.0)

^a Data from Ref. [15].

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 H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction. As mentioned above, this sampling scheme does not cover the PES region for the C₂H₅⁺ → H₂ + C₂H₃⁺ dissociation process, which will be discussed separately in terms of the vibrational energy distributions of C₂H₅⁺ generated from the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction. However, the C₂H₅⁺···H₂ 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₂H₅⁺···H₂ intermediate complex and key geometries for PES and MP2 are compared in Figure S1 in the Supplementary Materials. Although the H₂ 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₂H₅⁺ → H₂ + C₂H₃⁺ 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₂H₇⁺ → H₂ + C₂H₅⁺, 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].

2.2. 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₃⁺ + C₂H₄ reactants, along with their relative translational energy, were thermalized at each temperature in the simulations. In this scenario, the thermal rate coefficient $k\left(T\right)$ was obtained using the following equation [30]:

$$\begin{array}{c}k\left(T\right)=\sqrt{\frac{8k_{B}T}{\pi \mu }}\pi b_{max}^2\frac{N_r}{N_t} ,\end{array}$$

(1)

where $k_B$, $\mu ,$ and $b_{max}$ represent the Boltzmann constant, the reduced mass, and the maximum impact parameter, respectively. The standard error $∆k\left(T\right)$ is given as follows:

$$\begin{array}{c}∆k\left(T\right)=k\left(T\right){\left(\frac{N_t-N_r}{N_t{N}_r}\right)}^{\frac{1}{2}} ,\end{array}$$

(2)

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.

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.

	Classical				RPMD				KIDA
T (K)	(k(T) ± Δk(T)) × 10−9	bmax (Å)	Nr	Nt	(k(T) ± Δk(T)) × 10−9	bmax (Å)	Nr	Nt	k(T) × 10−9
100	2.59 ± 0.03	12.0	3245	4994	2.66 ± 0.09	13.5	453	859	-
150	2.86 ± 0.03	11.5	3183	4995	3.03 ± 0.09	13.0	514	972	-
200	3.14 ± 0.04	11.5	3028	4989	3.21 ± 0.07	12.5	1036	1972	-
250	3.31 ± 0.04	11.5	2856	4981	3.47 ± 0.10	12.0	531	964	-
300	3.43 ± 0.04	11.5	2703	4996	3.59 ± 0.07	11.5	1119	1975	2.9

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₂ + C₂H₅⁺ products. Notably, in our simulations, no trajectories led to 2H₂ + C₂H₃⁺ products, which is consistent with the anticipated high potential energy barrier in the C₂H₅⁺ → H₂ + C₂H₃⁺ reaction (refer to Figure 1b). Consequently, it is anticipated that the H₂ + C₂H₃⁺ products will form via an indirect process involving the vibrational excited C₂H₅⁺ fragment.

Figure 2 illustrates 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₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ and H₃⁺ + C₂H₄ → H₂ + H₂ + C₂H₃⁺ reactions, as the H₂ + C₂H₃⁺ 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₂H₂ → H₂ + C₂H⁻ and H₃⁺ + CO → H₂ + 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 charge-induced–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 H⁻ + C₂H₂ → H₂ + C₂H⁻ reaction [30]. In the following section, we will discuss the energy distribution of the H₂ + C₂H₅⁺ product fragments resulting from proton affinity, aiming to comprehend the indirect mechanism of the C₂H₅⁺ → H₂ + C₂H₃⁺ reaction.

2.3. Internal Energies of Fragments

To estimate the rovibrational energy distributions of the C₂H₅⁺ fragment concerning the indirect reaction C₂H₅⁺ → H₂ + C₂H₃⁺, we examined the system internal energies (ε) of the C₂H₅⁺ fragment derived using the following equation:

$$\begin{array}{c}\epsilon \left({\mathrm{C}}_2{\mathrm{H}}_5^{+}\right)=E_k \left({\mathrm{C}}_2{\mathrm{H}}_5^{+}\right)+\frac{1}{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{\sum }_{j=1}^{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}\left\{V_j\left( {{\mathrm{H}}_2+\mathrm{C}}_2{\mathrm{H}}_5^{+}\right)- V_{{\mathrm{H}}_2}\left(r_j\right)\right\},\end{array}$$

(3)

where $E_k \left({\mathrm{C}}_2{\mathrm{H}}_5^{+}\right)$ represents the vibrational and rotational energies for centroid velocities of the C₂H₅⁺ fragment, and $V_j\left( {{\mathrm{H}}_2+\mathrm{C}}_2{\mathrm{H}}_5^{+}\right)$ denotes the potential energy obtained from the MSA2 PES concerning the nuclear configuration of j-th bead. Additionally, $V_{{\mathrm{H}}_2}$ is a sixth-order polynomial function dependent on the internuclear distance (r) between the H₂ fragments, where the center-of-mass distance (R) between the H₂ and C₂H₅⁺ 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₂ and C₂H₅⁺ 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₂H₇⁺ 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₂H₅⁺ 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₂ + C₂H₃⁺ products, suggesting that the C₂H₅⁺ → H₂ + C₂H₃⁺ dissociation channel was not dominant, although the distribution at T = 300 K indicates that the C₂H₅⁺ fragment with larger internal energies can lead to the H₂ + C₂H₃⁺ 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₃⁺ and C₂H₄ reactants, along with their collision energy, influence the proton abstraction process of C₂H₄ from H₃⁺. Consequently, we infer that the branching fraction for [C₂H₅⁺]:[C₂H₃⁺] is strongly temperature-dependent, potentially leading to the overestimation of the branching fraction of the C₂H₃⁺ 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₂H₅⁺···H₂ 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₂H₅⁺. 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₂ fragment in the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction. The distributions were obtained using the following equations:

$$\begin{array}{c}{\epsilon }_{\mathrm{v}\mathrm{i}\mathrm{b}} \left({\mathrm{H}}_2\right)=E_k^{\mathrm{v}\mathrm{i}\mathrm{b}} \left({\mathrm{H}}_2\right)+\frac{1}{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{\sum }_{j=1}^{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{V}_{{\mathrm{H}}_2}\left(r_j\right) ,\end{array}$$

(4)

and

$$\begin{array}{c}{\epsilon }_{\mathrm{r}\mathrm{o}\mathrm{v}\mathrm{i}\mathrm{b}} \left({\mathrm{H}}_2\right)=E_k^{\mathrm{v}\mathrm{i}\mathrm{b}} \left({\mathrm{H}}_2\right)+E_k^{\mathrm{r}\mathrm{o}\mathrm{t}} \left({\mathrm{H}}_2\right)+\frac{1}{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{\sum }_{j=1}^{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{V}_{{\mathrm{H}}_2}\left(r_j\right) ,\end{array}$$

(5)

where $E_k^{\mathrm{v}\mathrm{i}\mathrm{b}} \left({\mathrm{H}}_2\right)$ and $E_k^{\mathrm{r}\mathrm{o}\mathrm{t}} \left({\mathrm{H}}_2\right)$ represent the vibrational and rotational energies for the centroid velocities of the H₂ molecule. Unlike the distribution of the system’s internal energy for the C₂H₅⁺ fragment, the internal energy of H₂, 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₂, 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:

$$\begin{array}{c}E \left(\mathrm{K}\mathrm{E}\mathrm{R}\right)=\frac{1}{2}{\mu }_{\mathrm{r}\mathrm{e}\mathrm{l}}{\overrightarrow{v}}_{\mathrm{r}\mathrm{e}\mathrm{l}}^2,\end{array}$$

(6)

where ${\mu }_{\mathrm{r}\mathrm{e}\mathrm{l}}$ and ${\overrightarrow{v}}_{\mathrm{r}\mathrm{e}\mathrm{l}}$ denote the reduced mass and relative velocity vector between the H₂ and C₂H₅⁺ fragments. The distributions broadened, reflecting the distribution of the internal energies of the H₂ and C₂H₅⁺ 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₂ + C₂H₅⁺ products of our PES was relatively higher than the MP2 results.

3. Methodology

3.1. Development of a Global Potential Energy Surface

All the quantum chemistry computations to construct the full-dimensional C₂H₇⁺ 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⁻ + C₂H₂ → H₂ + C₂H⁻, qualitatively reproduced the rate coefficient [30]. It should be noted that for the NH₃⁺ + H₂ → NH₃⁺ + 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₃⁺ + C₂H₄ reactants and the H₂ + C₂H₅⁺ 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₂ + C₂H₅⁺ products, without directly producing the 2H₂ + C₂H₃⁺ product. Consequently, data sampling of the H₂ + H₂ + C₂H₃⁺ 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₂H₅⁺··· H₂ 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⁻¹ 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₂H₅⁺··· H₂ 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 coupled-cluster 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.

3.2. 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\left(\mathbf{p},\mathbf{r}\right)$ [30] is defined as follows:

$$\begin{array}{c}H\left(\mathbf{p},\mathbf{r}\right)={\displaystyle \sum _{i=1}^{3n}}\left\{\frac{1}{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}}\frac{{\left(p_i^j\right)}^2}{2m_i}+\frac{m_i{N}_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{2{\beta }^2{\hslash }^2}{\displaystyle \sum _{j=1}^{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}}{\left(r_i^j-r_i^{j-1}\right)}^2\right\}+\frac{1}{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}{\displaystyle \sum _{j=1}^{N_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}}}V\left(r_1^j,\dots ,r_{3n}^j\right),\end{array}$$

(7)

where $\hslash$ and $\beta$ represent the reduced Planck constant and reciprocal temperature, respectively, where $\beta \equiv 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_{\mathrm{b}\mathrm{e}\mathrm{a}\mathrm{d}}$ represents the number of beads, while $r_i^j$ and $p_i^j$ 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₂H₄ and H₃⁺ 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⁴–10⁶ 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].

4. Conclusions

The thermal rate coefficients for the overall reaction, encompassing both the direct reaction H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ and the indirect reaction C₂H₅⁺ → 2H₂ + C₂H₃⁺, were determined across a temperature range from 100 to 300 K. The coefficients exhibited an increase as the temperature increased. The coefficient at 300 K was slightly higher than the value in the KIDA database, but we still reasonably reproduced it. Comparatively, the rate coefficients obtained from the RPMD simulations showed a slight increase compared to those from the classical MD simulations, which were attributed to fluctuations in the nuclear probability densities [30]. Moreover, in order to explore the indirect mechanism of the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ → 2H₂ + C₂H₃⁺ reaction, we analyzed the distributions of internal energies for the H₂ + C₂H₅⁺ fragments and the relative translational energy (E (KER)). Despite the influence of the proton affinity of H₃⁺ on generating the C₂H₅⁺ fragment, its internal energies were insufficient to reach the dissociation limit of the H₂ + C₂H₃⁺ products. While the ε of C₂H₅⁺ fragment decreased with decreasing temperature, the internal energy of H₂ and the E (KER) exhibited a weak correlation with the temperature. Based on these findings and the estimation of the internal energies of the C₂H₅⁺ fragments, we suggest that the H₃⁺ + C₂H₄ → H₂ + C₂H₅⁺ reaction prevailed, while the C₂H₅⁺ → 2H₂ + C₂H₃⁺ reaction was of lesser significance. Here, it should be mentioned that the comparison of theoretical rate coefficients at 300 K with the experimental results and the temperature dependency of the ε of C₂H₅⁺ fragment was investigated qualitatively using the developed PES. In the future, we recommend an experimental investigation into the temperature dependency of the branching fraction for [C₂H₅⁺]:[C₂H₃⁺], complemented by theoretical results based on a high-quality PES constructed using sophisticated techniques such as the Δ-machine learning algorithm [32,33,34].