Nonthermal rate constants for $CH4*$ + X $→$ CH₃ + HX, X = H, O, OH, and O₂

The chemical properties of energetic gas phase environments, such as homogeneous catalytic reactors, shockwaves, plasmas, engines, and flames, result from the competition of reactive collisions and nonreactive thermalizing collisions. In many contexts, inert collisions greatly outnumber reactive ones, and it would seem to be a safe assumption to treat species as thermalized in chemical modeling. Recent work, however, has indicated the importance of reactions involving nonthermalized species.^1–8 Notably, these modeling studies demonstrated that nonthermal reactions can influence macroscopic phenomena and complicate experimental interpretations. For example, relatively small changes in the radical pool due to the rapid nonthermal dissociation of hot HCO⁷ and other weakly bound radicals⁹ can have a measurable impact on predicted flame speeds, an important indicator of global reactivity in combustion systems.

In radical rich environments, nonthermal reactivity may be promoted via radical–radical association reactions, as the initial molecular adducts are highly energized. For example, the reaction

produces a CH₄^* adduct with an initial excess of ∼100 kcal/mol of internal energy. The assumption that, before it can react, the energized adduct thermalizes via so-called “third body” collisions with an inert species M,

is so commonplace in chemical modeling that Eqs. (1) and (2) are almost always written together as

Not all third bodies are inert, however, and here we distinguish reactive third bodies X from inert ones M. In hydrocarbon flames,^10,11 peak concentrations of the reactive radicals X = O, OH, and H can be as large as 0.5%–5%, depending on the physical conditions (pressure, temperature, and fuel stoichiometry) and reactive-flow processes (including both diffusion and chemistry) involved. Thermal rate constants k for radical abstraction reactions involving these radicals, i.e.,

increase with increasing temperature and range from 10⁻¹³ to 10⁻¹¹ cm³ molecule⁻¹ s⁻¹ at the temperatures relevant to flames and combustion chemistry (∼500–2000 K). While the analogous nonthermal rate constants k^* for

are not generally known, they may be anticipated to be significantly larger. For example, simply equating the vibrational energy of CH₄^* to a temperature (104 kcal/mol/9k_B = 5800 K, where we have partitioned the energy among methane’s nine vibrational modes, and k_B is Boltzmann’s constant) suggests that k^*/k may be quite large. Indeed, one might expect the rate constant at such a high temperature to be near the collision rate constant (10⁻¹⁰–10⁻⁹ cm³ molecule⁻¹ s⁻¹), with k^*/k then in the range from about 10 to 1000.

From a chemical modeling perspective, reactions (1) and (5) can be represented together as the termolecular reaction

The branching between thermalizing collisions in reaction (3) and reactive ones in (6) is then of considerable interest, as this branching alters the available radical pool and in turn can affect the overall reactivity of the system.

Here, quasiclassical trajectories are used to predict nonthermal rate constants k^* and to characterize their dependence on the initial internal energy of CH₄^* as well as the temperature of the remaining degrees of freedom. As abstractors, we consider the doublet and triplet radicals X = H, OH, and ³O, each of which can become relatively abundant in combustion environments. We also consider X = ³O₂, which has a much larger abstraction barrier (∼55 kcal/mol) and significantly lower thermal rate constants (<10⁻¹⁵ cm³ molecule⁻¹ s⁻¹ below 2000 K) than the other abstractors considered here but is, of course, present in relatively large concentrations in some environments (21% in air).

While we are not aware of any previous studies of bimolecular reactions involving a polyatomic reactant with state-averaged excitations close to its dissociation threshold, there have been considerable influential experimental and theoretical efforts in closely related areas. For example, the state-selected reactivity of excited methane^12–18 and other polyatomic species^19–24 has been extensively studied, and mechanistic details fundamental to chemical reactivity have been explored. Of particular relevance here is the relative importance of translational and internal energy in promoting reactivity that has been addressed for methane abstraction reactions,^25–27 with large abstraction rates observed experimentally with sufficient state-selected internal excitation of methane.²⁸ As another set of examples, detailed quasiclassical trajectory studies of atom–diatom reactions have been carried out to characterize reactivity at the nonthermal (hyperthermal) translational energies of practical relevance to the modeling of scramjets, shockwaves, orbiting bodies, and other nonequilibrium flows.^29–37 

Abstraction reactions involving methane are often used as prototypical chemical systems, and a great number of potential energy surfaces (PESs) describing these reactions have been developed. As an incomplete set of examples, we note Ref. 38, which compares some of the earliest PESs for CH₄ + H, Ref. 39, which reviews more recent work, and Refs. 13–16, and 18 as representative of fits from a few groups. These reactions are often used as benchmarks in the development and validation of new dynamics and kinetics theories,^40–42 including some of the most ambitious full quantum dynamics calculations to date.^39,43–45

Here, we use the analytic potential energy surfaces (PESs) of Corchado, Espinosa-Garcia, and co-workers for X = H,⁴⁶ OH,^47,48 and O.⁴⁹ The fitting strategy employed by these authors⁵⁰ is particularly useful in the present context, as the PESs are constructed from physically motivated analytic functional forms and can therefore be used with CH₄ excited well above thermal energies without promoting artificial dissociation (e.g., there are no turnovers at high-energies in the functional forms used to describe methane’s repulsive interactions). These PESs were designed to reproduce a variety of experimental and calculated properties of the reactants, products, and saddle points, and the original authors extensively validated them against experimental thermal kinetics, measured kinetic isotope effects, and collision-energy-dependent reactive cross sections. Some of these validation studies are summarized next. The PESs used here are all freely available via POTLIB-online.⁵¹ 

The CH₄ + H PES was taken from Ref. 46 where it was parameterized against CCSD(T)/cc-pVTZ energies. The fitted zero-point-exclusive (classical) barrier height (15.0 kcal/mol) was shown to differ from the computed one by just 0.1 kcal/mol, and subsequent⁵² variational transition state theory (TST) calculations including multidimensional tunneling reproduced experimental rate constants over a broad temperature range. Specifically, from 1000 K (the lowest temperature of relevance in the present work) to 2000 K (the highest temperature considered in Ref. 52), TST rate constants were shown to be within the error bars of two experimental recommendations,^53,54 differing from them by less than 15%.

Similarly, for CH₄ + OH, the fitted⁴⁸ classical barrier height of 6.40 kcal/mol is in close agreement with a subsequent CCSD(T)-F12/aug-cc-pVTZ calculation¹⁵ (6.29 kcal/mol), and, TST rate constants⁴⁷ were again found to be in close agreement with the curated experimental values of Baulch et al.,⁵³ differing by less than 20% from 1000 to 2000 K. For CH₄ + O,⁴⁹ the authors noted the small energy difference (∼0.2 kcal/mol) in the saddle point energies for the two reactive surfaces, and they stated that their surface was intended to represent reactivity on either surface. Their computed A′ (lower-energy) saddle point at the CCSD(T)/aug-cc-pVQZ//CCSD(T)/cc-pVTZ level of theory (14.0 kcal/mol) was reproduced well by their fitted surface (14.1 kcal/mol),⁴⁹ which was parameterized against global ab initio data computed at the same level of theory. We will further validate these PESs below via comparisons of our thermal quasiclassical trajectory results with experimental rate constants.

Direct dynamics was used for CH₄ + O₂. Several low-level quantum chemistry approaches with computational efficiencies suitable for direct dynamics were tested. The B3LYP/6-31G method was chosen, which predicts a 0 K reaction enthalpy of 55.8 kcal/mol, which compares favorably with the ANL0 computed value of 55.4 ± 0.3 kcal/mol.⁵⁵ For this reaction, Ref. 56 reported a CCSD(T)/aug-cc-pVTZ//CCSD(T)/aug-cc-pVDZ computed saddle point at −1.9 kcal/mol and a product van der Waals well at −2.6 kcal/mol relative to products and again including zero point corrections. The B3LYP/6-31G method does not predict a saddle point. Nevertheless, with the high energies and temperatures of interest here, and since the van der Waals well is only slightly bound relative to the saddle point (0.6 kcal/mol), the absence of this saddle point in our trajectory calculations is not expected to be a significant source of additional uncertainty.

Thermal rate constants were computed using the usual quasiclassical prescriptions⁵⁷ and the trajectory code DiNT.⁵⁸ Briefly, the rotational and vibrational states of the reactants were sampled from thermal distributions at T_vib and T_rot, respectively, and using the harmonic oscillator and rigid rotor approximations. The vibrational phase of each mode and the overall orientation of each reactant were randomized. The relative energy was sampled from a thermal distribution at T_rel, and the square of the impact parameter was sampled uniformly from 0 to b_max², where b_max is a parameter along with the initial and final center of mass distances R_i and R_f. Results were confirmed to be converged within our reported statistical uncertainties with respect to the choices of b_max = 6 Å, R_i = 10 Å, and R_f = 10 Å. Trajectories were propagated using a variable step size integrator⁵⁹ and integrator tolerances of 10⁻⁵ and 10⁻³ E_h (1 E_h = 627.5 kcal/mol) for the analytic potential energy surfaces and direct dynamics calculations, respectively.

When computing thermal rate constants, k(T), all three temperatures were set to the nominal temperature, T = T_vib = T_rot = T_rel, which defines the trajectory ensemble. Nonthermal ensembles are defined by two parameters: E_vib, the initial vibrational energy of the energized reactant, and T, the temperature of the remaining degrees of freedom. Nonthermal ensembles were generated as above but with T_vib for the energized reactant temporarily set to E_vib/(3n − 6)k_B, where n is the number of atoms in the molecule. Then, the potential energy of the molecule was computed at the sampled geometry and the vibrational component of the momentum was scaled such that the energized reactant contained an internal vibrational energy of exactly E_vib. In the rare cases where the potential energy exceeded the desired energy, the sampled geometry was abandoned and the procedure was restarted. The rotational state of CH₄ was sampled from an isolated thermal distribution at temperature T_rot and using a rotational constant for methane of 5.24 cm⁻¹. The instantaneous projection K along the symmetry axis is not conserved during vibrational motion, and the rotational energy associated with K evaluated at the quasiclassically sampled geometry of CH₄ was included in the scaling of E_vib. The remaining rotational energy, E_rot, is associated with the quantum number J, which is conserved during vibrational motion, and E_rot was not included when scaling E_vib. The initial total internal energy of CH₄ is therefore E_vib + E_rot. Nonthermal rate constants are labeled k^*(E_vib, T).

Reactive events were monitored based on C–H bond distances, and rate constants were computed using the standard formula⁵⁷ 

where μ is the reduced mass of the reactants, g_e is the ratio of the electronic partition functions for the reactive system and the reactants, and N_R/N is the fraction of reactive trajectories. The one-sigma statistical error in k was computed, again using standard approaches, as

and two-sigma error bars are reported throughout this work.

In ongoing work, comprehensive modeling studies of the effect of nonthermal kinetics on combustion observables are underway. In the meantime, here we modeled the competition between thermalization [reaction (3)] and nonthermal reaction [reaction (6)] at a single temperature, 1000 K, as follows: we assumed that the reactive environment consists only of radicals X and inert species M, with concentrations of n_X and n_M, respectively. The rate of undergoing a nonthermal reaction relative to the rate of collisions with M is

where Z is the CH₄ + M collision rate constant, and we have suppressed noting the T-dependence in p(E). The collisional cooling of CH₄ was modeled by the successive lowering of the energy E, and no attempt was made to distinguish vibrational and rotational energy or to accurately describe the distribution of E that occurs in real systems. Instead, we used

where z counts collisions and $ζ1/2=ln2ζ$ is the number of collisions required to cool CH₄ by half its energy. Equation (10) is equivalent to assuming first-order decay of the energy E with a rate constant given by Zn_M/ζ.

Nonreactive CH₄^* + N₂ trajectory calculations were carried out to predict the collisional cooling rate at T = 1000 K using methods described in detail elsewhere,^60–63 and the range parameter ζ of Eq. (10) was determined in two ways. In the first way, the average vibrational energy of an ensemble of 12 800 CH₄ + N₂ trajectories was monitored over 256 successive collisions, with collisions defined using the Lennard–Jones collision rate Z. The initial vibrational energy was set to D₀ = 104 kcal/mol, and the temperatures of the bath and of the rotational degrees of freedom of CH₄ were set to 1000 K. The average vibrational energy was found to decrease exponentially with the number of collisions z, from which the range parameter $ζ1/2$ = 165 was extracted. This value of $ζ1/2$ corresponds to a decrease in the vibrational energy of ∼150 cm⁻¹ per collision for early collisions where E ≈ D₀.

Because at 1000 K the vibrational temperature of CH₄ is effectively much hotter than its rotational temperature (and the bath), early collisions tend to heat up rotations (and the bath) such that per-collision decreases in the total rovibrational energy of CH₄ are smaller than per-collision decreases in the vibrational energy quantified above. For early collisions, the trajectory studies showed that this difference is roughly a factor of four. The simple collisional cooling model used here does not differentiate between vibrational and total energy, and so there is ambiguity in which cooling rate is more appropriate. Here, we simply show the sensitivity to this parameter by first generating results using the vibrational energy cooling rate, where $ζ1/2=165$⁠, and then again with a value of $ζ1/2$ that is four times larger corresponding to the slower cooling rate for the total energy.

The fraction of CH₄^* surviving after z collisions is

where F₀ is initially 1. Equation (11) converges with respect to z because p(E_z) → 0 as z → ∞, and therefore,

is this model’s estimate of the fraction of hot CH₄^* that reacts to form CH₃ + HX instead of cooling to CH₄. This model is intended to be qualitative, and more detailed modeling studies are underway.

We considered the accuracy of the following simple method for estimating nonthermal rate constants, k^*, from thermal ones, k. In this approach, an effective temperature, T_eff, is first estimated for the nonthermal system via

where

E_vib and m_vib are the vibrational energy and number of vibrations of the hot reactant, respectively, T is the temperature of the m = 3n − 3 − m_vib remaining degrees of freedom, and n is the total number of atoms of the reactive system. Equation (13) is simply the weighted average of the temperatures of all the degrees of freedom of the system (with vibrations contributing twice as much as rotations in accordance with the equipartition theorem and ignoring overall translation), with the energy of the nonthermal reactant converted to units of temperature via k_B as in Eq. (14).

The relationship in Eq. (14) does not reflect the distribution of energies that temperature is intended to represent. We therefore consider a second simple model where temperature is equated with the median energy of the exponential Boltzmann distribution instead of its mean energy as in Eq. (14). This has the effect of increasing T_vib by a factor of 1/ln 2 = 44%.

The effective temperature T_eff, whether estimated using the mean or median energy-temperature relationship, provides a means of approximating nonthermal rate constants from thermal ones via

so long as thermal rate constants are known to sufficiently high T_eff. The accuracy of this approximation will be quantified below by comparing k(T_eff) with our quasiclassical trajectory results for k^*.

Figure 1 compares quasiclassical and experimental thermal rate constants, k, for CH₄ + H as a function of T, as well as quasiclassical nonthermal rate constants, k^*. We first consider k. The present thermal results are in good agreement with the experimental recommendations of Baulch et al.,⁵³ especially at the higher temperatures included in their evaluation. At 2000 K, for example, the quasiclassical trajectory result (1.88 ± 0.15 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹, where a two-sigma statistical uncertainty is given) is just 15% larger than the experimental recommendation (1.63 × 10⁻¹¹ cm³ molecule⁻¹ s⁻¹), and this difference is in fact smaller than the assigned experimental uncertainty of Δlog₁₀ k = 0.2 or ∼50%.⁵³ 

At lower temperatures, errors in the quasiclassical rate constant increase and are 40% and 180% at 1500 K and 1000 K, respectively. Corchado et al.⁴⁶ validated the PES used here against experimental rate constants using variational transition state theory calculations including multidimensional tunneling, and they demonstrated quantitative agreement with experiment over a broad temperature range, including 1000 K. We therefore attribute the significant over-prediction of our quasiclassical trajectory results at T ≤ 1000 K to our use of classical mechanics (including errors due both to the neglect of tunneling as well as differences in the classical and quantum mechanical state counts and thresholds) and not to errors in the potential energy surface. One expects errors due to our neglect of these quantum mechanical effects to decrease rapidly with increasing temperature, as observed here.

We are primarily interested here in nonthermal rate constants, and, as discussed next, these are typically much larger than the 2000 K thermal rate constant. We therefore chose to not complicate the present calculations with semiclassical corrections for tunneling or zero point energy maintenance. Instead, we anticipate errors in our nonthermal rates to be similar to or smaller than the errors estimated here in the thermal rate constant at 2000 K, at least for rate constants larger than ∼10⁻¹¹ cm³ molecule⁻¹ s⁻¹. Larger relative errors may be present in the smaller nonthermal rate constants reported below, but these rates are sufficiently small that knowing them imprecisely should not affect the present analyses and conclusions.

Figure 1 shows the experimental recommendation of Baulch et al.,⁵³ including an extrapolation of this expression above 2500 K and outside of the temperature range of the experiments used to inform it. The extrapolated experimental rate constants and the present calculated rate constants continue to agree within a factor of 2.5, even at the highest temperature considered (8000 K). At high temperatures where no experimental information is available, we expect the quasiclassical result to be more accurate than the extrapolated experimental expression. Quantifying this discrepancy is useful in the context of the effective temperature model described by Eqs. (13)–(15) and validated below, as T_eff may exceed temperatures that have been accessed experimentally. This comparison shows that the experimental expression for k can continue to be reliably used (within a factor of 2–3) for large T_eff.

Nonthermal rate constants, k^*(E_vib, T), are shown in Fig. 1 for E_vib = D₀, D₀/2, D₀/4, and 0, where the threshold energy D₀ = 104 kcal/mol. The results are generally shown to be strong functions of both E_vib and T. For E_vib = D₀, which may be thought of as representing the reaction prior to any thermalizing collisions, the nonthermal rate constants are very large, even at 1000 K, with k^*(D₀, T) > 1.5 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹. Thermal rate constants this large are typically associated with barrierless radical–radical reactions (for CH₃ + H, for example, k ≈ 4 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹).⁶⁴ In both cases, the rate constants are an appreciable fraction of the collision rate.

For lower E_vib (i.e., as CH₄^* cools), k^*(E_vib, T) decreases, and the effect is larger at lower T, where D₀ is initially a larger fraction of the total energy available to the reaction. We see that k^*(E_vib, T) ≈ k(T) around T = 2500, 4000, and 6750 K for E_vib = D₀/4, D₀/2, and D₀, respectively. These temperatures may be rationalized via the effective temperature model described in Sec. II D. Specifically, Eq. (14) predicts crossover temperatures of 1800, 3270, and 6130 K for E_vib = D₀/4, D₀/2, and D₀, respectively (or 2230, 4130, and 7840 K using the modified median energy relationship) in fair agreement with the observed crossover points in Fig. 1. This model further predicts that these crossover temperatures should not depend strongly on the nature of the radical X, as will be shown in Sec. III B.

A more detailed test of the effective temperature model is given next. Figure 2 compares k^* computed via quasiclassical trajectories with k(T_eff) estimated using Eqs. (13)–(15). For E_vib = D₀ and D₀/2, the error in the effective temperature model is less than a factor of three, and the error is less than a factor of four for E_vib = D₀/4. Larger differences are found for E_vib = 0. Although these discrepancies are larger than the errors assigned here to the results of the quasiclassical trajectory calculations, errors of this size are likely small enough to provide semiquantitative guidance in modeling studies of nonthermal effects when no other information about k^* is available. The accuracy of the effective temperature model is empirically found to be better using the median energy relationship than the mean energy relationship.

Physically, deviations in the trajectory-based values of k^* from those obtained using the effective temperature approach are one measure of nonthermal and nonstatistical dynamics in CH₄^* + H. When E_vib and/or T are large, such deviations are small (less than a factor of three), presumably due both to the rapid sharing of energy among the modes during reaction and to the smallness of the saddle point energy relative to the total energy. At low E_vib and low T, the rate constants depend more sensitively on the initial distribution of energy among the modes and how strongly these modes are coupled to the reaction coordinate, larger nonstatistical effects are expected, and the effective temperature method fails qualitatively.

Figure 3 shows k^*(E_vib, T) as a function of E_vib for T = 1000 K. Notably, the nonthermal rate constants are significantly larger than the thermal rate constants, even after substantial cooling. The nonthermal enhancement is as large as k^*(E_vib, T)/k(T) = 707 for E_vib = D₀ (using the experimental thermal rate constant as the reference) and remains two orders of magnitude larger down to E_vib ≈ D₀/2. In Sec. II, we estimated the number of collisions required to reduce the vibrational and total rovibrational energies by half to be 165 and 660, respectively. Although this analysis is imprecise (in particular, by neglecting the distribution of energies central to the concept of temperature), it does suggest that nonthermal reaction rates will remain fast for on the order of hundreds of collisions.

The above analysis was repeated for X = O and OH, and similar results were obtained. Two of the three triplet states of O are reactive, and, as noted by its authors, the potential energy surface used here was designed to represent reaction occurring on either surface.⁴⁹ We estimated the total rate constant for CH₄ + O by multiplying our computed rate constants by (1 + exp(−ΔV^‡/k_BT)), where ΔV^‡ = 0.2 kcal/mol is the energy difference in the saddle point energies of the ground and first-excited states.⁴⁹ This correction was applied to both k and k^* and, due to the smallness of ΔV^‡, is nearly equal to simply doubling the rate constant.

The CH₄ + O reaction has been extensively studied experimentally and theoretically, and here we again restrict our validation study to comparisons with the recommended thermal rate constants of Baulch et al.⁵³ The calculated values of k are within 30% of the recommended experimental values up to 2500 K, the highest temperature included in the Baulch evaluation, and these differences are smaller than the factor of two uncertainty that they assigned. Transition state theory calculations by the developers of the PES used here (using an earlier but similar fit⁶⁵) predicted rate constants that are ∼30%–40% lower than the Baulch et al. expression (again reported up to 2500 K) and are in close agreement with the present trajectory calculation above 1500 K. For CH₄ + OH, the present quasiclassical trajectory results agree with Baulch et al.’s recommended values for k to better than 40% at 1000 K and 20% above 2000 K. These results are consistent with the error analysis given above for CH₄ + H, as one expects somewhat smaller errors in the quasiclassical approach for the heavier O and OH abstractors than for H, and this trend is observed here. The overall error in the dynamics and in the fitted potential energy surfaces for these systems is again likely small for rate constants larger than 10⁻¹¹ cm³ molecule⁻¹ s⁻¹.

Figures 4 and 5 show thermal and nonthermal rate constants for CH₄ + OH and CH₄ + O, respectively. The computed nonthermal rate constants are very large for highly-energized CH₄, with rate constants exceeding 3 × 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹ for E_vib = D₀. These values are even larger than rate constants for the analogous barrierless CH₃ + X reactions (both were previously calculated^66,67 to be ∼10⁻¹⁰ cm³ molecule⁻¹ s⁻¹).

Supporting the effective temperature model, the crossover temperatures at which k^*(E_vib, T) = k(T) are very similar for X = H, O, and OH, as anticipated above in our discussion of CH₄ + H. Although not included here, plots analogous to Fig. 2 for X = O and OH were made to test the effective temperature model and similar errors of a factor of 2–3 were found for E_vib ≥ D₀/4, with somewhat larger errors when both E_vib and T are low.

Thermal rate constants for X = O₂ are significantly lower than those for the other abstractors considered here, with k(1000 K) = 3 × 10⁻²² cm³ molecule⁻¹ s⁻¹.⁵⁶ Furthermore, our use of direct dynamics limits the number of trajectories that can be computed, and so we do not attempt a detailed study of the thermal and nonthermal rate constants for this abstractor. Instead, nonthermal rate constants were computed for E_vib = 104, 96, and 88 kcal/mol and T = 1000 K using ensemble sizes of 1000–2000 trajectories. The computed values of k^* at these three energies are 7.4 ± 2.6, 5.7 ± 2.6, and 2.8 ± 2.1, respectively, where the units are 10⁻¹² cm³ molecule⁻¹ s⁻¹ and two-sigma statistical uncertainties are given. Despite these relatively large uncertainties, these values again validate the median-energy effective temperature model, which here predicts nonthermal rate constants within a factor of ∼2 of the quasiclassical trajectory values of k^*. In the modeling studies below, we therefore used the effective temperature model for X = O₂ and divided the resulting rate constant by 1.8 to match our computed values at high E_vib.

Figure 6 compares the dependence of k^* on E_vib for the abstractors at T = 1000 K. In every case, k^* ≫ k until CH₄ is cooled nearly completely. Prior to any cooling, the nonthermal enhancement at T = 1000 K is at least two orders of magnitude, with k^*(D₀, T)/k(T) = 707, 143, and 360 for X = H, OH, and O, respectively, and the nonthermal rate constants are enhanced relative to the thermal ones by more than a factor of 10 down to E_vib = D₀/4. For low E_vib, k^* is shown to increase nearly linearly with energy up to ∼60 kcal/mol for X = OH and O and up to ∼D₀ for X = H. At high energies, the nonthermal rates are relatively independent of E_vib, are roughly equal to one another, and have magnitudes typically associated with collision rate constants. The curves suggest a transition in the dynamical bottleneck from the more conventional picture of energy exchange among the modes until the reaction barrier is surmounted at lower energies to collision rate–controlled reactivity at high energies.

We have much less data for X = O₂. The nonthermal rate constants k^* for X = O₂ are 30–100 times smaller than those for the other abstractors studied here (H, O, and OH), even for energies near threshold. These values of k^* for X = O₂ are ∼10¹⁰ larger than the thermal rate constants for this reaction at 1000 K, however.

The data in Fig. 6 were used along with the estimates for the collisional cooling rate discussed in Sec. II to evaluate P, the fraction of hot CH₄^* that reacts with X before cooling [Eqs. (9)–(12)]. In Fig. 7, P is shown at 1000 K as a function of the mole fraction of X. Typical peak mole fractions of O, OH, and H in methane-air flames can approach 0.5%, 1%, and 2% respectively,^10,11 and these are marked on the curves in Fig. 7. At these values, we see that our model predicts that nonthermal abstractions will divert between 35% and 50% of CH₄ formed via H + CH₃ to CH₃ + HX before CH₄^* can be stabilized to CH₄. For a given mole fraction, the effect of OH and O atoms is roughly twice that of H atoms due to the nonthermal rate constants for H being half as large as those for O and OH at high energies. H atoms can build up in higher concentrations than O and OH, however, and so nonthermal effects from all three abstractors may be anticipated generally to have similar impacts in combustion modeling.

Although the nonthermal rates considered here have been computed as functions of the vibrational energy of CH₄, E_vib, a more complete model would consider both the vibrational and rotational energies of CH₄. In particular, we noted in Sec. II that the total rovibrational energy of CH₄ initially cools approximately four times more slowly than its vibrational energy. If we apply our same model but with a cooling rate that is four times smaller, we find that P is increased by roughly a factor of two, leading to an even more complete conversion of CH₄^* to CH₃ + HX via the nonthermal reaction channels. This model likely overestimates the true effect of including rotations into the cooling model, but it is suggestive that the present predictions may in fact underestimate the importance of the nonthermal channels considered here.

Figure 8 shows P for X = O₂ as a function of the mole fraction of O₂. Although the nonthermal rate constants for this reaction are much smaller than those for the other abstractors, the increased relative abundance of O₂ (21% in air relative to 0.5%–2% for the other radicals) leads to the conversion of roughly 10% of CH₄^* via the nonthermal channel, which is ∼$14$ of the effect noted for the other systems.

Overall, the branching fractions in Figs. 7 and 8 suggest that nonthermal reactions may have a non-negligible effect on combustion simulations, as flame speeds are often sensitive to the available radical pool. Preliminary master equation simulations were performed to calculate rate constants for CH₄^* + X, X = H, O, OH, O₂, as a function of T, pressure, and n_X. These results were represented as the termolecular processes

for use in detailed chemical kinetic modeling, and preliminary modeling studies of laminar flame speeds in adiabatic freely-propagating CH₄-air flames at standard conditions (1 atm and 300 K) were performed. Noticeable increases in laminar flame speeds were observed when these termolecular reactions were added along with the corresponding lowering of the stabilization rates k₃. The latter are known to be sensitive chain termination reactions in CH₄-air flames, and we can attribute the majority of the changes in the predicted laminar flame speeds to the reduction in k₃ by 20%–50%. In fact, if we artificially include the termolecular channels without reducing k₃, a small decrease in the flame speed is observed. Detailed simulations are in progress to quantify the effects of nonthermal reactions more broadly and to better represent these processes in detailed kinetics models. Similar termolecular processes may be important for the adducts of other small radical–small radical recombination reactions involving common combustion radicals, such as H + OH + M/X, OH + OH + M/X, H + HCO + M/X, CH₃ + OH + M/X, and CH₃ + CH₃ + M/X.

Quasiclassical trajectory calculations, along with the previously published potential energy surfaces of Corchado, Espinosa-Garcia, and co-workers, were used to compute nonthermal rate constants for CH₄ + X (X = H, O, and OH) as a function of temperature and of the vibrational energy of methane, E_vib. For large values of E_vib, CH₄^* was found to react with X as rapidly as a comparably sized radical (CH₃), with rate constants close to those of barrierless reactions and approaching the collision limit (>10⁻¹⁰ cm³ molecule⁻¹ s⁻¹). These nonthermal rate constants represent an enhancement for near-threshold values of E_vib of more than two orders of magnitude at 1000 K.

This nonthermal rate enhancement was rationalized by computing an effective temperature T_eff from E_vib and the temperatures of the other degrees of freedom and noting that k^*(E_vib, T) ≈ k(T_eff), at least within a factor of 2–3. This model suggests that significant nonthermal effects are less likely for larger systems and for radical adducts, where T_eff is reduced by the increased number of degrees of freedom and the typically smaller bond energies, respectively. These considerations may limit the potential importance of nonthermal effects of the type considered here in combustion to reactions involving just a few small molecules, e.g., CH₄^*, CH₂O^*, C₂H₆^*, H₂O^*, H₂O₂^*, and CH₃OH^*.

By varying the internal energy of CH₄^*, E_vib, we found that large nonthermal rate constants persisted even after ∼100 s of thermalizing collisions and significant cooling of CH₄. Each abstractor showed a somewhat different dependence on E_vib, but in general k^* was shown to increase linearly with E_vib at low energies until k^* ≈ 10⁻¹⁰ cm³ molecule⁻¹ s⁻¹. At higher energies, where k^* is presumably collision-limited, k^* is only weakly dependent on E_vib.

The competition between nonthermal reaction and collisional cooling was studied at 1000 K using a simple model. We showed that due to the large nonthermal rates for H, O, and OH, as well as their weak dependence on collisional cooling at large E_vib, nonthermal reactions can account for from one-third to one-half of the fate of CH₄^* formed via H + CH₃ when the concentrations of the abstracting radicals are at their peak concentrations of ∼0.5%–2%.

A more limited study of X = O₂ was carried out. Although the nonthermal rate constants k^* are smaller for this system than for the other abstractors, the relative enhancement k^*/k is many orders of magnitude and the abundance of O₂ in flames is higher. Altogether, we found that around 10% of CH₄^* was diverted to CH₃ + HO₂ instead of stabilized CH₄ at 1000 K, which could nonetheless be an importance source of radicals under some conditions.

More detailed modeling studies are underway to quantify the importance of these nonthermal reactions in combustion simulations of practical relevance. Preliminary results suggest that these nonthermal channels lead to measurable changes in flame speeds under conditions where H + CH₃ reactions are important.

An important limitation on the accuracy of the branching fractions predicted here is the model used to describe collisional cooling. This potential source of error was indicated by the fairly large discrepancy in the solid and dashed curves in Figs. 7 and 8, which were generated using the two different cooling rates representing either vibrational or total rovibrational cooling. Detailed models of collisional cooling parameterized using classical trajectories have been shown to provide quantitative descriptions of collisional energy transfer for a priori predictions of thermal rate constants, but only when rotational and vibrational cooling are both explicitly treated.⁶³ Similarly, a more accurate model for nonthermal reactivity might track both energy and angular momentum while cooling. This could be achieved, for example, by monitoring the fate of an ensemble of independent CH₄^* molecules as they stochastically encounter collisions with X and M and where the outcome of each collision is determined via explicit bimolecular trajectories.

Finally, we note that in studies of similar systems^68–70 large collision energies were shown to promote the high-energy process X + CH₄ → H + XCH₃, which were not considered here. Although it seems less likely that internal excitations would be as effective as high translational energies at promoting these processes, it may be worthwhile to consider these channels in future studies.

This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences through Argonne National Laboratory. Argonne is a U.S. Department of Energy laboratory managed by UChicago Argonne, LLC, under Contract No. DE-AC02-06CH11357. We gratefully acknowledge computing resources provided by Bebop, a high-performance computing cluster operated by the Laboratory Computing Resource Center at Argonne National Laboratory. Support for R.S. and S.J.K. was provided as part of the Argonne–Sandia Consortium on High-Pressure Combustion Chemistry (Grant No. ANL FWP 59044).