Reactions Of OH Radicals With Benzene And Toluene In The Gas-Phase

Aromatic hydrocarbons such as benzene and toluene are ubiquitous both in natural and anthropogenic emissions, and they play a central role in atmospheric chemistry, combustion processes, and environmental pollution.^1,2 Among the reactive species in the atmosphere, the hydroxyl radical (^•OH) is widely recognized as the primary “detergent” of volatile organic compounds (VOCs), initiating their degradation via either hydrogen abstraction or addition to the aromatic ring.^3,4

The reactions of ^•OH with benzene and toluene are particularly important, because these compounds are not only common VOCs, but also precursors to secondary pollutants such as ozone and secondary organic aerosols (SOA), especially via radical-initiated oxidation pathways. Moreover, fundamental understanding of their kinetics and mechanistic details under a wide range of conditions (temperature, pressure) remains critical for atmospheric modeling, combustion mechanisms, and toxicological risk assessment.^5,6 Historically, experimental studies have measured absolute rate constants for OH + benzene and OH + toluene over a range of temperatures and pressures, revealing a mechanism that changes with conditions: at low temperature (near ambient), addition of ^•OH to the aromatic ring tends to dominate; but at elevated temperatures (T > 500 K), decomposition of the adduct becomes rapid and the net reaction shifts toward hydrogen-abstraction channels (ring H-abstraction for benzene, side-chain H-abstraction for toluene). This temperature- (and pressure-) dependence has substantial implications for combustion chemistry as well as for high-temperature atmospheric environments (e.g., urban pollution, thermal degradation).^5-8 However, despite decades of research, accurate theoretical descriptions remain challenging. The branching between addition and abstraction pathways, and the relative importance of different abstraction channels (e.g., ring vs. methyl-side chain) for toluene, have been subject to significant uncertainties due to sensitivity on computational method, treatment of anharmonicity, and the potential for pressure- and temperature-dependent kinetics. Recent high-level theoretical studies combining sophisticated quantum-chemical calculations with master-equation (RRKM/ME) or variational transition-state theories reveal that treatment of anharmonic vibrational modes, tunneling, and multistructural effects can dramatically influence the predicted rate coefficients and branching ratios.^9,10 In this context, the present work employs high-level ab initio methods (CCSD(T)) and density functional theory (DFT), together with canonical and statistical kinetic treatments (TST and RRKM/ME), to compute rate constants for OH + benzene and OH + toluene over a broad temperature range (300–2000 K). The aim is to elucidate, in a unified framework, which reactive channels dominate under which conditions, and to quantify the influence of anharmonicity (particularly low-frequency bending modes in the transition states) on the predicted kinetics. Such a study not only refines fundamental mechanistic understanding, but also provides benchmark data useful for atmospheric modeling, combustion mechanisms, and for comparison with experimental measurements. Specifically, we show that for OH + benzene, hydrogen-atom abstraction becomes the dominant pathway under the studied conditions, while for OH + toluene, the ipso-H abstraction channel prevails. Moreover, the inclusion of anharmonic treatment for the RC–H–OH bending vibrations in transition states is found to be crucial for reliable calculation of rate constants. Our findings highlight the importance of proper theoretical treatment (method level, anharmonicity, statistical kinetics) when modeling radical–aromatic reactions, and suggest that similar care should be taken for other H-abstraction systems.

Finally, by providing temperature-dependent Arrhenius expressions for both title reactions, this work offers data that can be directly implemented in kinetic models for atmospheric chemistry and combustion, thereby contributing to more accurate predictions of VOC fate, pollutant formation, and radical-driven reaction networks.

COMPUTATIONAL METHODS

Electronic structure calculations

All quantum-chemical calculations were carried out using the Gaussian 16 package.¹¹ First, the geometries of all species involved in the reaction systems (reactants, transition states, intermediates, and products) were optimized at the DFT level using the B3LYP functional in conjunction with the 6-311++G(3df,2p) basis set.^12,13 Harmonic vibrational frequency analyses and moments of inertia were also computed at the same level to verify the nature of the stationary points (real frequencies for minima; one imaginary frequency for transition states) and to obtain thermochemical data (zero-point energies, vibrational partition functions, rotational constants, etc.). To improve the accuracy of energetic parameters (e.g., relative energies and reaction barriers), single-point energy calculations were then performed at the CCSD(T) level using the same 6-311++G(3df,2p) basis set, based on the DFT-optimized geometries. The CCSD(T) method,¹⁴  often considered the “gold standard” for accurate ab initio energetics for closed-shell molecules,  accounts for electron correlation beyond mean-field, including a perturbative treatment of triple excitations, and thus provides high-accuracy relative energies. We checked the reliability of the coupled-cluster wavefunctions by evaluating T₁ diagnostics and spin-contamination (< S²>) for all relevant species; only results satisfying standard thresholds (e.g., T₁ below an accepted limit, negligible spin contamination) were considered for subsequent kinetic modeling.¹⁵

Kinetic modeling: TST and RRKM / Master Equation treatment

Rate constants for all reaction routes (e.g., H-abstraction, possible addition or other channels) were computed using both conventional Transition State Theory (TST) and statistical RRKM/Master-Equation (ME) approaches^16,17 to account for pressure and energy-transfer effects. For the TST calculations, we used the ChemRate code; for RRKM/ME modeling we employed the Mesmer code. In the partition-function evaluation, besides the standard rigid-rotor and harmonic-oscillator approximations, we treated low-frequency modes (in particular torsional modes such as C–CH₃ bond rotations) using a hindered-rotor (HIR) model rather than the pure harmonic oscillator approximation. This choice helps to avoid overestimate of vibrational partition functions (and thus to improve the accuracy of calculated rate constants) when the harmonic model becomes invalid for floppy or highly anharmonic motions. For the Master Equation modeling of pressure- and collision-induced energy transfer, we adopted a downward energy transfer parameter <?E_down> = 400 cm^-1. The bath gas was assumed to be N₂; the Lennard–Jones (L–J) parameters used for N₂ were ?/k = 97.530 K and σ = 3.621 Å. For intermediate species (e.g., adducts or radical intermediates), we employed L–J parameters based on those of a proxy molecule: Phenol (C₆H₅OH), namely ?/k = 464.8 K and σ = 5.29 Å.¹⁸ These parameters define collisional frequencies and energy transfer probabilities in the kinetic simulations. Sum of states (SoS) and density of states (DoS) for each species were computed internally by the Mesmer code, using the standard algorithms implemented therein.

RESULTS AND DISCUSSION

Mechanism and kinetics of the C₆H₆ + OH reaction

The potential energy surface (PES) for the C₆H₆ + OH system has been constructed to account for both possible reaction pathways: the H-abstraction channel and the OH-addition channel. The energy profile (in kcal/mol, relative to the reactants as zero) is shown in Figure 1.

Figure 1. PES of the C₆H₆ + OH reaction at the CCSD(T)//B3LYP level of theory. Energy values are in units of kcal/mol.

According to the H-abstraction channel, the reaction proceeds via a relatively “tight” transition state (TS1) with an activation barrier of 2.8 kcal/mol. Once this barrier is overcome, the system evolves to yield the bimolecular products C₆H₅ + H₂O. Importantly, these products lie at -5.5 kcal/mol relative to the reactants, indicating that the overall H-abstraction process is exothermic and thermodynamically favorable under the conditions considered. In contrast, the OH-addition route initially forms a pre-reaction complex (COM1), lying at 2.6 kcal/mol below the reactant asymptote, rather than crossing a classical barrier at the first step. From this complex, the system can proceed through a shallow transition state (TS2, 0.1 kcal/mol) to form a stable adduct (IS1) at -17 kcal/mol. Thus, the addition of OH to the aromatic ring is energetically favored at this intermediate-formation stage. However, to evolve from the adduct (IS1) toward the final bimolecular products (C₆H₅OH + H), the reaction must surmount a relatively high barrier via TS3, located at 8.3 kcal/mol. This barrier makes the addition channel significantly less favorable for product formation compared to the abstraction pathway, at least on the static PES. Overall, from the energetic standpoint, the H-abstraction channel is both kinetically and thermodynamically more favorable for producing stable products (C₆H₅ + H₂O) than the OH-addition channel is for yielding C₆H₅OH + H. The favorable energetic profile of the abstraction pathway (low barrier and exothermic product formation) suggests that H-abstraction will dominate under many conditions, especially when sufficient energy is available to overcome TS1. Nevertheless, the presence of a deep adduct well (IS1 at -17 kcal/mol) in the addition route indicates that under conditions where collisional stabilization is efficient (e.g., higher pressures, rapid energy transfer), the adduct may be long-lived; under such conditions, addition might compete, or lead to alternative reaction fates. This observation underlines the importance of kinetic modeling (e.g., RRKM / Master-Equation analysis) to evaluate the real contribution of each channel under realistic environmental conditions, beyond the static PES picture. Moreover, the contrast between a “tight” barrier for abstraction and a “barrierless/low-barrier” formation of the addition adduct (but a high barrier for its conversion to products) reflects a common pattern in radical–aromatic reaction mechanisms. This also connects with classical ideas regarding saddle-points and minima on a PES: while minima correspond to stable (or metastable) species, transition states correspond to first-order saddle points along the reaction coordinate. To calculate rate constants for the reaction system, the TST model was used for the H-abstraction channel, while the RRKM/ME approach was utilized for the OH-addition channel. The resulting rate constants were fitted to a modified Arrhenius expression; for example, for the OH + benzene reaction we obtain: k(OH + benzene) = 7.8.0×10^-11 exp(-6.36 kcal.mol^-1/RT) cm³/molecule/s.

Mechanism and kinetics of the C₆H₅CH₃ + OH reaction