Canonical Unified Statistical Rate Constants Using a High-Level Composite Coupled Cluster Energy Scheme for the CH 4 + CH Reaction

We report an ab initio study on the kinetics and chemical dynamics of the CH 4 + CH reaction using a coupled cluster based composite scheme that includes unrestricted coupled cluster singles and doubles with perturbative connected triples method, UCCSD(T), energies extrapolated to the complete basis set limit, in addition to core-valence, higher-order correlation, and relativistic corrections. With this protocol, the reaction enthalpies for the two main reaction channels, C 2 H 4 + H and CH 2 + CH 3 , are reproduced with an absolute deviation of 0.04 kcal mol –1 relative to experimental data. At the UCCSD(T)/complete basis set (CBS) level of theory, the formation of the C 2 H 5 intermediate is barrierless, in contrast with the low submerged barrier (–1.87 kcal -mol –1 relative to the reactants) found at the UCCSD(T)/augmented correlation-consistent polarized double-zeta (aug-cc-pVDZ) level of theory. With the correlation-consistent polarized triple-zeta (cc-pVTZ) basis set, we describe structural variations for the reaction bottleneck along the reaction path, finding at least three canonical variational transition states for the temperature range from 20 to 800 K. The thermal rates were obtained via the canonical unified statistical theory (CUS), using the canonical variational transition state theory (VTST) for the inner-transition state and long-range transition state theory (LRTST) for the outer-transition state. Our calculations agree with literature measurements and show inverse Arrhenius behavior, as observed experimentally. At 298 K, the computed rate constant is 2.65 × 10 –10 cm 3 molecule –1 s −1 and reported experimental measurements range from 2.5 × 10 −12 to 3.0 × 10 –10 cm 3 molecule –1 s −1 . Our theoretical study represents an improvement on previous computational investigations and highlights that even relatively simple gas-phase reactions can require high levels of theory to be modeled accurately.


Introduction
Climate change is one of the greatest concerns for humanity's future, and multiple pieces of evidence point towards a worsening of the phenomenon. 1 Particularly relevant to this matter is the increasing concentration of atmospheric greenhouse gases (GHG) as a result of anthropic activities. Methane gas is one of the main GHG which plays an important role in the thermal balance of the planet. 2 It has about 25 times more global warming potential than carbon dioxide and its concentration is increasing by 3.65% every decade. 1 The industrial conversion of methane into value-added organic products arose as an alternative to reduce methane atmospheric concentrations, but, given its high thermodynamic stability, usual chemical processes often demand a high energy consumption. For example, in the so-called dry reforming of methane (DRM), CH 4 + CO 2 → 2 CO + 2 H 2 , the overall endothermicity is 247 kJ mol −1 , hence requiring high temperatures (900-1200 K). 3 In this regard, technologies such as nonthermal plasmas (NTP) offer the advantage of keeping the gases close to room temperature, thus eliminating the problems of degradation of the organic products formed. [4][5][6][7] The chemical activation of a target molecule by electron collisions with NTP technologies provides a good alternative for stable molecule transformations such as methane.
The title reaction also is present in combustion's reactions of methane-air systems, 8 playing a key role in the astrochemistry modeling of planetary atmospheres like Titan. 9 Its experimental rates show an inverse temperature dependence within the range from 23 to 772 K with C 2 H 2 + H as the major products. [10][11][12][13][14][15][16][17][18] The other reaction channel results in CH 2 + CH 3 as products. Anti-Arrhenius behavior can occur for barrierless reactions, or by the presence of one transition state submerged at the potential energy surface (PES), 19,20 due to tunneling, 21 or even by the influence of catalysts. 22 Several computational works aiming to understand the dynamics and kinetics of the CH 4 + CH reaction have been reported in the literature. 10,[23][24][25] Ribeiro and Mebel 25 demonstrated successfully that the reaction occurs via two stages, being, in summary, the long-range reactants approximation forming an intermediate and the proper H-migration from CH 4 to CH radical, that can lead to both products, C 2 H 4 + H and CH 3 + CH 2 . Using the MP5/6-311++G(d,p)//MP4(SDQ)/6-311++G(d,p) level of theory and conventional transition state theory, they obtained rate constants that are between factors of three and six from experimental measurements in the 145-582 K temperature range following the anti-Arrhenius behavior. Furthermore, they described a negative barrier height of 2.6 kcal mol -1 at the CCSD(T) level with explicitly correlated (F12) approaches and complete basis set (CBS) limit extrapolation for single point energies calculated at optimized geometries with density functional theory (DFT) with the B3LYP density functional and the 6-311G(d,p) basis set. Zero-point energies (ZPE) were calculated at the B3LYP/6-311(d,p) level of theory.
In this work, we report calculated thermal rate constants for the reaction in the 20-800 K temperature range using variational transition-state theory (VTST), 26 long-range transition-state theory 27 (LRTST) and canonical unified statistical (CUS) theory. [28][29][30] These calculations are based on an accurate determination of the energetics of the reaction, at a level of theory based on unrestricted coupled cluster singles and doubles with perturbative connected triples method, UCCSD(T), with an estimate of the complete basis set (CBS) limit, scalar relativistic effects, core-valence correlation, full triple and quadruple excitations, using a procedure similar to the one used to determine accurate barrier heights for other gas-phase reactions. [31][32][33][34] We have identified variations on the dynamical bottleneck of reaction, changing the transition state (TS) structure from an early-TS to a late-TS in the investigated temperature range, that allows us to report rate constants that are in much better agreement with experiments than previous computational investigations. 10,24,25 Methodology Characterization of the stationary points Geometry optimizations of stationary points were accomplished using spin-unrestricted coupled-cluster theory with singles, doubles, and perturbative triples clusters 31 (UCCSD(T)) with restricted open-shell Hartree-Fock (ROHF) wave functions. The cc-pVTZ basis set of Dunning 32 was selected for this purpose. To check the nature of the stationary points, harmonic vibrational frequency calculations were performed at the optimized structures using the same level of theory: minima were identified as having no imaginary frequencies, while transitions states with one imaginary frequency. The zero-point energies (ZPE) of all minima and transition states were determined at the UCCSD(T)/cc-pVTZ level of theory.

Electronic energies
At each stationary point on the PES, the total electronic energy, E, was obtained via a composite scheme that loosely follows the Feller-Peterson-Dixon approach and has been shown to provide accurate energy barriers for gas-phase reactions. [33][34][35][36] The present approach includes core-valence correlation, scalar relativistic effects, an estimate to the complete basis set (CBS) limit, as well as corrections due to connected triple triple (CCSDT) and quadruple (CCSDTQ) excitations; see equation 1. The core-valence correlation (∆ CV ) was included at UCCSD(T)/aug-cc-pwCVQZ level of theory. 32,37,38 The scalar relativistic effects (∆ rel ) were included by second-order Douglas-Kroll-Hess 39 (DKH) Hamiltonian with aug-cc-pwCV5Z-DK basis set. 40,41 The reference Hartree-Fock (HF) energy calculations, with the aug-cc-pVnZ (n = T, Q, and 5) basis sets, and the UCCSD(T)/CBS correlation energy, with the aug-cc-pVnZ (n = Q and 5), were extrapolated to the complete basis set (CBS) limit according to equations 2 and 3, respectively. Finally, the triple (∆ T ) and quadruple (∆ Q ) excitations of cluster operator were considered by differences between UCCSDT and UCCSD(T), using cc-pVTZ basis sets for the first case, together with the same operation between UCCSDTQ and UCCSDT, with cc-pVDZ basis sets for the last one. Table S1 in Supplementary Information (SI) section presents the additive energy corrections of all stationary points. All UCCSD(T) calculations were carried out with the Molpro 42,43 software suite, using the MRCC, 44,45 interfaced with Molpro, for UCCSDT and UCCSDTQ calculations. (1) where n is the cardinal number of the correlation consistent basis set (cc-pVDZ = 2, cc-pVTZ = 3, cc-pVQZ = 4, and so on), A and C are extrapolation coefficients of the least squares method.
Spin-orbit coupling for the CH radical The electronic partition functions for the CH radical were obtained using equation 4, considering the spin-orbit (SO) coupling. The internal energy, U el , and entropy, S el , contributions for the same species and degree of freedom were computed using equations 5 and 6 for the same SO effect. Electronic structure programs, as Molpro,42,43 usually only consider the spin-multiplicity contribution to the electronic partition function, q el . The thermochemical quantities computed by Molpro for the CH radical were modified to include the temperature dependent contribution of the spin-orbit effect to the Gibbs free energy (ΔG). 19 (4) where k B is the Boltzmann constant, R is the ideal gas constant, h is the Planck constant, c is the speed of light, and A is the experimental spin-orbit coupling constant, which in our case is equal to 27.9 cm -1 . 46 With these equations, we could obtain Gibbs free energies for the reactants.

Rate constants calculations
Using our final composite electronic energies, we obtained rate constants for the CH 4 + CH reaction. The VTST was applied using the data displayed in Table S3 (SI section). The minimization of recrossing effects was done considering molecular geometries along the reaction path that was approximated by displacing TS1 according to the normal mode with imaginary frequency. The molecular structure with maximum Gibbs free energy at a given temperature along the reaction path defines the canonical variational transition state.
In the low-temperature regime, the long-range interactions between the reactants play an important role in barrierless reactions. 27 We applied equations 7 and 8 using the reactants molecular properties given in Table S4 to obtain the corresponding long-range transition state theory 27 rate constant for the dispersion interaction potential between the reactants.
where µ is the reduced mass, α i is the polarizability, and E i is the ionization energy of the reactants. The dipole-induced dipole forces were neglected in this work due to the small value of the electric dipole moment of the CH radical.
The rates k disp were multiplied by the degeneracy factor, f e , 47 to consider the SO effect of CH (equation 9).
The low and high-temperature regimes were unified through the CUS theory given by equation 10. (10)

Results and Discussion
As a first assessment of our composite method, we have calculated reaction enthalpies (∆ r H) at 0 K for the involved processes. For the CH 4 + CH → C 2 H 4 reaction, ∆ r H is -59.55 kcal mol -1 , while the value of 3.47 kcal mol -1 is obtained for the CH 4 + CH → CH 3 + CH 2 process. These values agree within 0.04 kcal mol -1 to gas-phase experimental data. 48 The availability of experimental geometric parameters of the reactants and products allows us to evaluate the calculated structures at UCCSD(T)/cc-pVTZ level of theory, according to Table 1. In general, the outcomes are in good agreement with the experimental data, with absolute deviations smaller than 0.01 Å for bond distances and smaller than 0.4° for bond angles. Such accurate geometries obtained at the UCCSD(T)/cc-pVTZ level indicate that they are excellent reference structures for the refinement of total electronic energies via the composite method and for the application of TST.
The energy diagram along the reaction path is shown in Figure 1. As seen, the outcome obtained using the composite method at UCCSD(T)/CBS + ∆ rel + ∆ CV + ∆ T + ∆ Q level of theory indicates that, differently from UCCSD(T)/cc-pVTZ calculations, there is an inversion between TS1 and IM1 relative energies. Our high-level calculations indicate that the PES does not have a saddle point in the entrance channel, and the minimum energy reaction path that connects the reactants to the intermediate C 2 H 5 is barrierless.
Analyzing the relative energy inversion between intermediate 1 (IM1) and transition state 1 (TS1) in more detail, we note that the biggest contribution to the lowering of the energy of TS1 comes from the increase of the oneelectron basis sets, with minimal contributions from the higher-order N-particle basis set corrections, i.e., UCCSDT and UCCSDTQ (∆ T and ∆ Q , respectively). Figure S1a shows that the IM1-TS1 barrier height varies from 0.29 kcal mol -1 at UCCSD(T)/cc-pVDZ level to 0.25 kcal mol -1 at UCCSDTQ/cc-pVDZ. Figure S1b evidences that at the UCCSD(T)/aug-cc-pVDZ level TS1 is 0.12 kcal mol -1 higher than IM1; and at the UCCSD(T)/CBS level TS1 lies 0.08 kcal mol -1 below IM1. The improvement of the level of theory leads to a barrierless reaction, i.e., falling directly to the global minimum of the PES. This explains our unsuccessful initial attempts to locate this stationary point with higher levels of theory. Indeed, if we use the B3LYP PES, both transition states are connected to the minimum energy path (MEP) through intrinsic reaction coordinate (IRC) 50,51 as shown in Figure S2 presented in SI section.
The geometry variation of the first TS is shown in Figure 2, which shows three transition states in the temperature range from 20 to 800 K. In fact, from 20 to 230 K the maximum free energy is in the region of negative displacement of the normal coordinate, characterizing an early-TS. 52 Between 235 to 335 K the same maximum is   closer to the proper saddle point TS1. Finally, from 340 to 800 K there is a positive displacement for the same analysis, configuring a late transition state. 52 The best estimate for the computed rate constant, obtained by the CUS method, is shown as an Arrhenius plot in Figure 3, along with the complete set of experimental measurements available in the literature. A more complete version of this plot, with rate constants computed with conventional transition state theory (CTST), VTST, LRTST and CUS, and experimental measurements individualized by reference is available as Figure S3, in SI section. Furthermore, Table 2 allows a quantitative comparison between computed rate constants and selected experimental data available in the literature.
As shown in Figure 3, the computational rate constants display the anti-Arrhenius behavior as observed experimentally. To compare our calculated rate constants with experimental measurements, it is convenient to divide the 20-800 K temperature range into three parts that show different behavior as indicated by our calculations. First, the high-temperature range, from 400 to 800 K, is well described by the VTST. This is consistent with the fact that, at these temperatures, the dynamical bottleneck of the reaction should be associated with a tight inner transition state in the vicinity of the B3LYP saddle point; see Figure S2 (SI section). In fact, the experimental value with the highest temperature available is 4.42 × 10 -11 ± 0.22 × 10 -11 cm 3 molecule -1 s -1 at 772 K 11 while the k VTST is 2.01 × 10 -11 cm 3 molecule -1 s -1 . Ribeiro and Mebel 25 also reported tight-TS and rate constants that are close to the experiment for these temperatures (473 and 581 K).
The low-temperature region, from 20 to 100 K, is very well described by LRTST as the long-range potential determines the rate constant. The experimental value with the lowest temperature available is 3.16 × 10 -10 ± 0.25 × 10 -10 cm 3 molecule -1 s -1 at 23 K, 14 while k disp is 3.76 × 10 -10 cm 3 molecule -1 s -1 . For intermediate temperatures, between 100 and 400 K, both inner-and outer-TS contribute to the determination of the rate constant. Therefore, it is the CUS theory that leads to the closest agreement with the experiment.

Conclusions
To summarize, the electronic energy description at UCCSD(T)/CBS +∆ rel + ∆ CV + ∆ T + ∆ Q level provides rate  The best estimate for the computed rate constant, obtained by the CUS method, is shown as a solid line. The complete set of experimental measurements, including error bars, available in the literature [10][11][12][13][14][15][16][17][18]53,54 is displayed as points. The anti-Arrhenius behavior is noticeable for temperatures higher than 150 K, for computed rate constants, and higher than 100 K, for experimental rate constants, as the rate constants decrease with increasing temperature.
constants within a factor of 0.5-4 from experimental data and points to a barrierless reaction that is the reason for the anti-Arrhenius behavior. Our transition state calculations showed that the dynamical bottleneck for the reaction is very dependent on the temperature. At the low-temperature limit, the bottleneck corresponds to a loose outer long-range transition state that is well described by LRTST based on the dispersion long-range potential between the reactants. At the high-temperature limit, the bottleneck corresponds to a tight inner transition state. A balanced description of the contribution of these transition states over the full 20 to 800 K temperature range is provided by the CUS theory, yielding rate constants that agree with experimental measurements.

Supplementary Information
Supplementary information (all calculated thermal rate constants, additive energies corrections for stationary points, molecular properties information used in our system, experimental and computational enthalpies of reaction, and details of potential energy diagram of IM1 and TS1 stationary points) is available free of charge at http://jbcs.sbq.org.br as PDF file.