Temperature and Pressure-Dependent Rate ... - ACS Publications

May 15, 2015 - Temperature and Pressure-Dependent Rate Coefficients for the. Reaction of Vinyl Radical with Molecular Oxygen. C. Franklin Goldsmith,*...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCA

Temperature and Pressure-Dependent Rate Coefficients for the Reaction of Vinyl Radical with Molecular Oxygen C. Franklin Goldsmith,*,†,‡ Lawrence B. Harding,† Yuri Georgievskii,† James A. Miller,† and Stephen J. Klippenstein† †

Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, Illinois 60439, United States School of Engineering, Brown University, Providence, Rhode Island 02912, United States



S Supporting Information *

ABSTRACT: State-of-the-art calculations of the C2H3O2 potential energy surface are presented. A new method is described for computing the interaction potential for R + O2 reactions. The method, which combines accurate determination of the quartet potential along the doublet minimum energy path with multireference calculations of the doublet/quartet splitting, decreases the uncertainty in the doublet potential and thence the rate constants by more than a factor of 2. The temperature- and pressuredependent rate coefficients are computed using variable reaction coordinate transition-state theory, variational transition-state theory, and conventional transition-state theory, as implemented in a new RRKM/ME code. The main bimolecular product channels are CH2O + HCO at lower temperatures and CH2CHO + O at higher temperatures. Above 10 atm, the collisional stabilization of CH2CHOO directly competes with these two product channels. CH2CHOO decomposes primarily to CH2O + HCO. The next two most significant bimolecular products are OCHCHO + H and 3CHCHO + OH, and not C2H2 + HO2. C2H3 + O2 will be predominantly chain branching above 1700 K. Uncertainty analysis is presented for the two most important transition states. The uncertainties in these two barrier heights result in a significant uncertainty in the temperature at which CH2CHO + O overtakes all other product channels.

1. INTRODUCTION The vinyl radical (C2H3) is a key intermediate in hightemperature oxidation and pyrolysis. At high temperatures, alkyl radicals rapidly fragment to yield smaller radicals and a corresponding alkene. This process terminates with the formation of ethylene (C2H4). The cleavage of this final carbon−carbon bond to produce CO or CO2 occurs through the oxidation of either vinyl or the ketenyl radical (HCCO): C2H3 + O2 → products

(R1)

C2H3 → C2H 2 + H

(R2)

C2H 2 + O → HCCO + H

(R3)

© XXXX American Chemical Society

HCCO + O2 → CO + CO2 + H

(R4a)

HCCO + O2 → CO + CO + OH

(R4b)

The C2H3 + O2 reaction has numerous product channels, some of which are chain branching, and some of which are chain propagating. Understanding the branching among the Special Issue: 100 Years of Combustion Kinetics at Argonne: A Festschrift for Lawrence B. Harding, Joe V. Michael, and Albert F. Wagner Received: February 2, 2015 Revised: May 13, 2015

A

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

of CO:CO2 to be 17:1. Oguchi et al. used laser-induced fluorescence at room temperature and reported a branching fraction for CH2CHO + O of 0.14 to 0.21 between 10 and 100 Torr.27 Chishima et al. used cavity ring-down spectroscopy (CRDS) to measure the pressure dependence of the formation of CH2CHOO at room temperature and pressures between 10 and 100 Torr.28 The most recent experimental data to date are cavity-ring down spectroscopy measurements for C2H3 + O2 and HCO + O2 by Matsugi and Miyoshi.29

various product channels, as well as the competition between vinyl oxidation and decomposition is key in describing the growth of the radical pool in high-temperature oxidation. Vinyl can also add to the double bond in alkenes, which leads to the molecular growth of dienes and ultimately to soot formation.1,2 In this regard, the branching between fuel oxidation versus molecular growth depends upon the fate of vinyl, and accurate rate coefficients for the competing reactions are essential in predictive modeling. Notably, we calculate that the largest sensitivity coefficient for benzene formation in a C2H4 burnerstabilized flame with an equivalence ratio between 1.0 and 1.5 is for C2H3 + O2. The temperature- and pressure-dependent rate coefficients for vinyl decomposition,3−5 vinyl addition,6−8 and ketenyl oxidation9 are well characterized. The rate coefficients for vinyl oxidation, in contrast, remain poorly understood. A thorough, quantitative analysis is complicated by the complexity of the potential energy surface (PES), which (in addition to the initial adduct, vinylperoxy (CH2CHOO)) has several bimolecular product channels: C2H3 + O2 → CH 2CHOO

C2H3 + O2 → CH 2O + HCO → CH 2O + H + CO

(R1b)

→ CH 2CHO + O

(R1c)

→ C2H 2 + HO2

(R1d)

→ OCHCHO + H

(R1e)

→ CH 2CO + OH

(R1f)

→ CH 2O + HCO

(R1g)

→ CO + CH3O

(R1h)

→ CO2 + CH3

(R1i)

(R1g‐ii)

The C−H bond strength in HCO is a mere 14.7 kcal/mol, so CH2O + H + CO is still 72.6 kcal/mol below C2H3 + O2 in energy. Thus, a large fraction of the HCO will dissociate promptly to H + CO. In an earlier work by some of the coauthors, direct dynamics calculations were initiated for two transition states on the C2H3O2 potential energy surface.30 They reported that approximately 70% of the HCO will directly dissociate. Thus, although the reaction is typically written as C2H3 + O2 → CH2O + HCO, it is arguably more accurate to represent it as C2H3 + O2 → CH2O + H + CO. Similarly, though less significantly, the vibrationally excited methoxy in CH3O + CO also could dissociate prior to collisional stabilization, thereby yielding the same trimolecular products, CH2O + H + CO. For the purposes of this manuscript, however, we will continue to use the standard notation of incipient bimolecular products. A quantitative determination of the temperature dependence of the branching between HCO and H + CO is beyond the scope of this paper, but it will be the subject of a future direct dynamics investigation. The original measurements of Gutman and co-workers, as well as the more recent measurements by Eskola and Timonen, report ion counts for m/z of 29 and 30 for HCO and CH2O, respectively. However, the interpretation of the HCO signal was for a relative HCO yield, because accurate radical cross sections were not available. Thus, this technique can confirm that HCO and CH2O were formed with the same characteristic rise time, but it cannot determine if they are formed with the same yield, and therefore it cannot measure the branching fraction between CH2O + HCO and CH2O + H + CO. To the best of our knowledge, the only experimental data that are capable of distinguishing these two channels are from Matsugi and Miyoshi, who report that a significant percent of the incipient HCO will dissociate prior to collisional stabilization. According to their CRDS measurements, the branching fraction of HCO, reaction R1g-i, at room temperature and 20 Torr was 0.22 ± 0.03. Matsugi and Miyoshi also included a statistical theoretical prediction of the HCO yield, based upon semiclassical approximations to the HCO density of states, and obtained a prediction of 35%. One consequence of the theoretical analyses was that CH2CHO + O could be the primary product channel at higher temperatures. RRKM/ME calculations in ref 21 predicted that the rate constant for CH2CHO + O would exceed that for CH2O + HCO at 800 K in 760 Torr of N2. The location of this crossover temperature has profound consequences for kinetic modeling, because CH2CHO + O is chain branching but nearly thermoneutral (e.g., −7 kcal/mol), whereas CH2O + HCO is chain propagating but highly exothermic (−87 kcal/mol). Indeed, attempts to model ethylene oxidation in a high-

(R1a)

→ 3CHCHO + OH

(R1g‐i)

The net rate of C2H3 + O2 → products was measured at room temperature by Krueger and Weitz10 and by Fahr and Laufer,11 and later by Knyazev and Slagle over the temperature range of 300 to 1000 K.12 The dominant product for C2H3 + O2 originally was assumed to be C2H2 + HO2, on the basis of analogy with C2H5 + O2.13 However, Baldwin and Walker were unable to detect acetylene in flow cell experiments of C2H4 in reactive mixtures of H2 + O2.14 Instead, the stable products at 753 K were CO and CH2O, which lead them to suggest that a direct pathway to CH2O + HCO was the likely product of this reaction. This hypothesis was confirmed by Gutman and coworkers for temperatures between 297 and 602 K;15,16 they suggested the possibility of a four-member ring intermediate that would lead to the formation of CH2O + HCO. Several theoretical analyses of the potential energy surface attempted to explain the prompt formation of CH2O + HCO.17−23 In ref 21, Mebel et al. concluded that the most likely intermediate following vinylperoxy was a three-member ring, rather than a four-membered ring. Furthermore, the theoretical analyses suggested that additional, highly exothermic bimolecular products are possible. Eskola and Timonen provided further confirmation that CH2O + HCO is the dominant product channel between 200 and 362 K.24 Time-resolved FTIR emission spectroscopy revealed that CO2 + CH3, and possibly C2H2 + HO2, CH2CHO + O, and OCHCHO + H, are potential product channels;25,26 no rate coefficients were reported, but they estimated the ratio B

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 1. C2H3O2 potential energy surface. The pathways to the main product channel at lower temperatures, CH2O + HCO, are shown in black; the pathway to the main product channel at higher temperatures, CH2CHO + O, is shown in blue.

pressure flow reactor using the rate coefficients of ref 21 overpredicted the rate of fuel consumption by a factor of 7.31 The previously published product-specific rate coefficients for this system are Mebel et al.,21 Carriere et al.,31 Lopez et al.,23 and Matsugi and Miyoshi.29 As will be seen later, these studies differ substantially in their prediction of the individual rate constants, in some instances by orders of magnitude. The potential energy surface (PES) is shown in Figure 1. The branching between CH2O + HCO at lower temperatures and CH2CHO + O at higher temperatures is sensitive to two key transition states. The C2H3 + O2 → CH2O + HCO pathway proceeds through at least four intermediates (shown in black), and the kinetic bottleneck is the isomerization from cCH(CH2)OO (dioxiranylmethyl) to c-CH(O)CH2O (oxiranyloxy). The C2H3 + O2 → CH2CHO + O pathway (shown in blue), in contrast, proceeds through a single intermediate, and the kinetic bottleneck is the O−O bond fission in CH2CHOO. These two transition states have proven to be difficult to analyze using standard computational quantum chemistry methods. Mebel et al. used the G2M(RCC,MP2) method to compute the PES in ref 21. Carpenter demonstrated that the transition state for dioxiranylmethyl isomerization to oxiranyloxy has a strong multireference character.32 Using multireference perturbation theory (CASPT2), Carpenter obtained a transition-state barrier height that was 11.5 kcal/mol lower in energy than that originally obtained in ref 21. By lowering the barrier height for this transition state, we find that CH2O + HCO persists as the dominant product channel for significantly higher temperatures. In fact, if the c-CH(CH2)OO → cCH(O)CH2O barrier is lowered by this amount, then the CH2CHO + O channel becomes kinetically irrelevant, because at temperatures sufficiently high to activate this channel, vinyl decomposition will dominate oxygen addition. In a response to Carpenter, Mebel and Kislov used both coupled cluster [CCSD(T)] and multireference configuration interaction (MRCI) with larger basis sets to study this transition state in greater detail.33 They obtained transition-barrier heights that were 7−8 kcal/mol lower in energy than in ref 21. The energies obtained by Mebel and Kislov were computed using the same geometries as by Carpenter in ref 32. Although Carpenter presents CASPT2(23,15) energies, the geometry

optimization was done using CASSCF(23,15). Despite the large active space, CASSCF is nonetheless a poor method for geometry optimization, and the frequencies obtained from CASSCF are not reliable for kinetics. Consequently, neither Carpenter’s nor Mebel and Kislov’s barrier height can be assumed to be quantitatively accurate for this critical transition state. The transition state for CH2CHOO → CH2CHO + O has received less attention, but it too is problematic, owing to a low-lying excited state and the nature of the bond fission. A more rigorous treatment of these two key transition states is critical to understanding the branching fraction between CH2CHO + O versus CH2O + HCO. The present work presents the most thorough analysis to date of the C2H3O2 potential energy surface. These highaccuracy electronic structure results are combined with state-ofthe-art RRKM/ME methods to predict the temperature- and pressure-dependent rate coefficients for all elementary reactions on the C2H3O2 PES. Uncertainty analysis is performed to propagate the uncertainties in key transition states into the final rate constants and product branching fractions.

2. COMPUTATIONAL METHODS 2.1. Electronic Structure. For each stationary point on the PES shown in Figure 1, geometry optimization and normalmode analysis were performed using UCCSD(T)/cc-pVTZ with an RHF reference wave function. For each stationary point, all possible conformers were considered. Five postCCSD(T) corrections were applied to the electronic energy, similar to those included in HEAT34 and W4,35 to reduce the uncertainty at two standard deviations to within 0.3 kcal/mol for minima. Estimation of the uncertainties for saddle points is considerably more difficult, and they are assumed to be larger. The total zero-point corrected electronic energy relative to C2H3 + O2, ΔECC,total, was computed as follows. First the CCSD(T) energy at the complete basis-set limit, ΔECBS, was extrapolated from aug-cc-pVQZ and aug-cc-pV5Z basis set calculations, assuming an inverse power law; 36,37 the augmented basis sets were modified to include the diffuse functions for only the s and p orbitals in C and O and the s function in H. Next, corrections for higher level excitations, C

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A ΔET/Q, were computed from the difference between CCSDT(Q)/cc-pVDZ and CCSD(T)/cc-pVDZ energies. Corrections for core−valence correlations, ΔECV, were determined from the difference in the CCSD(T)/cc-pcV∞Z (based upon the extrapolations from TZ and QZ calculations) energies where the core electrons are frozen (uncorrelated) or not. Relativistic effects, ΔErel, are approximated from the difference in CI energy with and without using the Douglas−Kroll one-electron integrals.38 Lastly, corrections for anharmonicity, ΔEanh, were approximated from the difference in the B3LYP/6-311+ +G(d,p) harmonic zero-point energy and that obtained using a second-order perturbative anharmonic analysis.39 All wave function calculations were performed in MOLPRO;40 the density functional theory calculations were performed in Gaussian09.41

are required, CASPT2 is chosen over MRCI for its computational efficiency. In each energy calculation, the relative orientation of the two fragments frozen at their equilibrium geometries is adjusted, and the 1D corrections are applied to the resulting energy. For CH2CHOO → CH2CHO + O, variational transitionstate theory was used to optimize the location of the dynamic bottleneck for both the ground state and the first electronic excited state. For all other transition states, conventional transition-state theory was used.44,45 For these transition states, tunneling was approximated using an asymmetric Eckart model. For both the sum and density of states, the nontorsional vibrational modes were assumed to be independent rigid-rotor harmonic oscillators. Torsional modes were represented as 1D, and in some cases 2D, internal rotors, with rotational potentials calculated in 10° increments using M06-2X/cc-pVTZ. 2.3. RRKM/ME. The temperature- and pressure-dependent phenomenological rate coefficients were computed using the newly developed RRKM/ME code, PAPER.46,47 A key feature of PAPER is the ability to automatically lump two or more wells into a single species when the corresponding eigenvalue for isomerization crosses the threshold from chemically significant eigenvalues into the fast internal relaxation eigenspace. It also has procedures in place for treating the numerical instabilities that arise when various parts of the potential have vastly different thermochemistries, such as when some stationary points are highly exothermic relative to the entrance channel. These features are critical when rate constants are computed for multiwell systems in which isomerization and/or decomposition barriers are greatly submerged relative to the inlet channel. Clearly that is the case for the present C2H3O2 system, where most of the right-hand side of Figure 1 is more than 60 kcal/mol below the C2H3 + O2 reactants. Collisional energy transfer was represented by a single exponential down model, given by 200(T/300)0.85 cm−1, based upon prior work for similarly sized systems.48,49 A LennardJones (LJ) collision model was used to determine the collision frequency. The LJ parameters for the bath gas, σHe = 2.55 Å and εHe = 6.96 cm−1, were taken from literature,50 and the LJ parameters for the C2H3O2 intermediates, σC2H3O2 = 5.18 Å and εC2H3O2 = 285.2 cm−1, were taken from ref 48.

ΔECC,total = ΔECBS + ΔE T/Q + ΔECV + ΔErel + ΔEanh + ΔZPECCSD(T)/cc ‐ pVTZ,harmonic ΔECBS = ΔECCSD(T)/aV5Z + (ΔECCSD(T)/aV5Z 54 − ΔECCSD(T)/aVQZ) 4 6 − 54 ΔE T/Q = ΔECCSDT(Q)/cc ‐ pVDZ − ΔECCSD(T)/cc ‐ pVDZ ΔECV = ΔECCSD(T)/cc ‐ pV ∞ Z,all − ΔECCSD(T)/cc ‐ pV ∞ Z,frozen core ΔErel = ΔECI/aug ‐ cc ‐ pcVTZ,DK − ΔECI/aug ‐ cc ‐ pcVTZ,non ‐ relativistic ΔEanh = ΔZPE B3LYP/6 ‐ 311 ++G(d,p),anharmonic − ΔZPE B3LYP/6 ‐ 311 ++G(d,p),harmonic (E1)

Multireference effects were significant for six transition states. The selection of the appropriate active space for a transition state is highly problem dependent, so a description of the approach used for each transition state is provided in the Results section below. For the multireference calculations, the complete basis set limit was extrapolated from cc-pVTZ and ccpVQZ single-point calculations: ΔEcc ‐ pV ∞ z = ΔEcc ‐ pVQZ − (ΔEcc ‐ pVTZ − ΔEcc ‐ pVQZ)

44 5 − 44 4

(E2)

3. RESULTS 3.1. Electronic Structure. The zero-point corrected electronic energy, including the individual corrections, for the intermediates and products of the C2H3O2 potential energy are listed in Table 1; the corresponding transition-state barrier heights are listed in Table 2. Also listed is the CCSD(T)/ccpVTZ T1 diagnostic,51 which provides some indication of multireference effects (i.e., a T1 diagnostic of >0.04 implies strong multireference effects, whereas a value >0.03 implies some cause for concern). The CCSD(T)/cc-pVTZ optimized geometries are provided in the Supporting Information. Of the transition states listed in Table 2, six of them required multireference theory to quantify the barrier height: TS 1, TS 2, TS 9, TS 10, TS 11, and TS 20. These transition states are described individually below. The T1 diagnostics are also greater than 0.03 for TS 3, TS 5, TS 6, TS 7, TS 8, TS 12, TS 15, and TS 17; however, these transition states are of comparatively minor importance to the overall kinetics, and the deviations from 0.03 are more modest. Thus, multireference aspects were not explored for those transition states.

2.2. Microcanonical Rate Constants. Microcanonical rate constants were computed for each transition state in Figure 1. For two transition states, C2H3 + O2 → CH2CHOO and CHCHOOH → 3CHCHO + OH, variable reaction coordinate transition-state theory (VRC-TST) was used.42,43 In VRC-TST, the internal degrees of freedom for the evaluation of the partition function are separated into conserved and transitional modes. The conserved modes are the normal modes of the two fragments at infinite separation, which are treated as standard rigid-rotor harmonic oscillators. The transitional modes, in contrast, correspond to the coupled, anharmonic motion of the two fragments, which are evaluated using Monte Carlo integration in classical phase space. The minimum energy path is computed first for the frozen geometries and again allowing for relaxation. The difference in energy for these calculations is used as a one-dimensional (1D) relaxation correction for the reference energies. Similar corrections are obtained for basis set extrapolation and active space size. The interaction potential is calculated on-the-fly using multireference theory. Because thousands of single-point calculations D

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

considered for the capture to form CH2CHOO is 9 electrons in 7 orbitals; the selected orbitals for this active space include the 4 π, π* orbitals in O2, the 2 π, π* orbitals in C2H3, and the carbon radical. The minimum energy path was calculated as a function of C−O distance, from r = 1.39 Å (the A″ CH2CHOO equilibrium geometry) to r = 20 Å. For each fixed distance r, the remaining degrees of freedom were optimized using both CASPT2(9,7)/cc-pVTZ and MRCI+Q(9,7)/cc-pVTZ, where +Q refers to the Davidson correction. The CASPT2 calculations included a level shift of 0.2 and an IPEA shift of 0.25.52 Additional single-point calculations were performed on the (9,7) geometries to test the effects of larger active spaces. Specifically, CASPT2(11,9) and MRCI+Q(11,9) calculations included the O 2 σ, σ* orbitals, and CASPT2(13,11) calculations also included the C−C σ, σ* orbitals (the latter active space was computationally prohibitive for MRCI+Q). The results of these calculations are shown as the broken lines in Figure 2. For R + O2 reactions, the dynamic bottleneck typically is in the region of 2−5 Å. Although the CASPT2 and MRCI+Q calculations each appear to be individually converged with respect to the active space, the two methods consistently differ by roughly 1 kcal/mol in the dynamic bottleneck region, with the CASPT2 potentials being more attractive. In a separate work, some of the authors demonstrated that the error introduced by the internal contractions used in the MRCI implementation in MOLPRO is greatest in the bottleneck region, and the present work appears to confirm this finding.53 To reduce the uncertainty in the interaction potential for R + O2 reactions, a new method was developed, on the basis of the doublet/quartet splitting. For each fixed distance r, two sets of single-point calculations were performed on the quartet surface using the optimized doublet geometries. In the first set of calculations, the doublet/quartet splitting was computed using the same multireference method as for the doublet geometry

Table 1. Local Minima on the C2H3O2 Potential Energy Surface

3.1.1. TS 1: CH2CH + O2 → CH2CHOO. The addition of molecular oxygen to the vinyl radical is barrierless, and multireference methods are required to describe the interaction potential far from the equilibrium geometry (e.g., a C−O distance of 2.0 Å or greater). The minimum active space

Table 2. Transition-State Barrier Heights on the C2H3O2 Potential Energy Surface energy relative to vinyl + O2 (kcal/mol) no. TS TS TS TS TS TS TS TS TS TS TS TS TS TS TS TS TS TS TS TS

1a 2 3 4 5 6 7 8 9a 10b 11b 12 13 14 15 16 17 18 19 20b

reaction R to W1 W1 to P2 W1 to P3 W1 to P5 W1 to W2 W1 to W3 W1 to W4 W2 to P3 W2 to P1 W3 to W5 W4 to W6 W5 to W7 W5 to W6 W7 to W8 W7 to P6 W8 to P7 W8 to P8 W6 to P6 W6 to P4 R to P3

T1

total

CBS

T/Q

CV

rel

anh

ZPE

0.029 0.037 0.023 0.031 0.038 0.039 0.034

−7.5 −1.2 6.5 −5.7 −20.1 −7.8 2.2

−8.1 −1.2 6.5 −6.1 −23.1 −9.7 1.5

−1.6 −0.3 −0.5 −0.6 −0.7 −1.9 −0.8

0.2 −0.1 0.1 0.1 0.2 0.2 0.0

0.0 0.1 0.1 0.2 0.1 0.1 0.2

0.0 0.2 0.1 0.0 0.0 −0.1 0.0

2.1 0.1 0.2 0.8 3.3 3.4 1.3

−84.6 −83.3 −71.1 −78.7 −90.0 −97.1 −83.3 −64.6

0.5 0.1 0.4 −0.6 0.4 −0.2 −0.1 0.0

0.0 0.1 0.1 0.1 0.3 0.0 0.1 −0.1

0.3 0.2 0.3 0.2 0.0 0.2 0.2 0.2

0.4 0.4 0.4 0.4 0.4 0.4 0.4 0.4

3.9 4.7 2.0 2.8 3.1 2.9 2.1 −1.0

0.035 0.022 0.024 0.042 0.028 0.045 0.024 0.022

−15.7 −21.2 −79.5 −77.7 −68.0 −75.9 −85.7 −93.8 −80.6 −65.0 10.3

a

These transition states are barrierless. Their partition functions were computed using VRC-TST. bThese barrier heights are based upon multireference calculations and not the coupled-cluster methods outlined above. For details, see the subsections below. E

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

the oxygen. At infinite separation, the CH2CHO + O complex is 3-fold degenerate. As the two fragments approach the saddle point, one of the configurations becomes strongly repulsive, and the transition state has a 2-fold, near degeneracy. In the dynamic bottleneck region, the near degeneracy increases to a difference of 1−2 kcal/mol. Once the fragments have cleared the saddle point, the lowest energy surface leads to the A″ adduct, and the higher energy surface leads to the A′ adduct. Both the ground state and the first excited state are included in the transition state. However, the model does not explicitly include the A′ adduct, as it is assumed to equilibrate with the A″ state. The CCOO dihedral angle is approximately 100° in the bottleneck region. As the fragments approach the CH2CHOO equilibrium geometry, this dihedral angle flattens to 180°. The low-lying excited state in the saddle point corresponds to the formation of A′ CH2CHOO, which is still attractive with respect to CH2CHO + O, Figure 3.

Figure 2. Potential for the C2H3 + O2 reaction. The broken lines correspond to the potentials computed using multireference methods exclusively. The solid lines with symbols correspond to the potentials computed using the new compound method, E3. Blue represents MRCI+Q, red CASPT2. All the potentials are extrapolated to the basis set limit.

optimization. In the second set of calculations, coupled-cluster calculations were performed. Although restricted Hartree−Fock calculations for the doublet state will fail for bond fission reactions sufficiently far from the equilibrium geometry, the same methods will yield a properly converged reference wave function for the quartet state. Thus, high-accuracy energies can be computed for the quartet surface using the coupled-cluster based methodology described in eq E1. The new, more accurate doublet potential is obtained by subtracting the multireference doublet/quartet splitting from the high-accuracy single-reference quartet energies:

Figure 3. CH2CHO + O potential energy surface, relative to CH2CH + O2. The inset shows the potential around the saddle point.

The ground-state potential was computed using the coupledcluster methods described by eq E1; the results are listed under TS 2 in Table 2. The excited-state potential was computed by first calculating the difference in energy between the two states using multireference methods and then adding this difference to the high-accuracy, single-reference ground-state energy. For the splitting between the surfaces, CASPT2(7,6)/cc-pVTZ was used for geometry optimization and single-point energies. The active space of 7 electrons in 6 orbitals corresponds to the 2 π, π* orbitals and the radical orbital in CH2CHO, plus the 3 2p orbitals of the O atom. These calculations, and all subsequent single-point calculations, were averaged over three states to facilitate convergence at long separations. For the single-point calculations that included the vinoxy lone pair, six-state averaging was required to ensure that the lone pair was consistently included in both the transition state and the products. All multistate calculations were performed using a dynamical weighting factor of 3.0.55 The T1-diagnostic for the CCSD(T)/cc-pVTZ optimized saddle point is 0.029, which indicates that single-reference methods should be appropriate for an open-shell system transition state. However, given that the transition state is still quite loose and involves what amounts to a bond fission process, it is nonetheless worth considering whether or not the single-reference method is consistent with multireference results. In addition to the (7,6) minimum active space, the following three combinations of orbitals were considered. The (9,7) active space included the lone pair of the vinoxy oxygen (averaged over six states), the (9,8) active space included the C−O σ, σ* orbitals in vinoxy (averaged over three states), and the (11,9) active space included both the C−O σ, σ* and the

ΔEtotal,doublet = ΔECC,total,quartet − (ΔEmultireference,quartet − ΔEmultireference,doublet)

(E3)

The multireference methods considered in eq E3 were CASPT2(13,11)/cc-pV∞Z//CASPT2(9,7)/cc-pVTZ and MRCI+Q(11,9)/cc-pV∞Z//MRCI+Q(9,7)/cc-pVTZ. The new potentials are shown as the solid lines in Figure 2. As seen in the figure, the CASPT2 and MRCI+Q potentials now agree to within 0.5 kcal/mol in the key dynamic bottleneck region, which reduces the uncertainty in the interaction potential by a factor of 2. The percent deviation in the interaction potential roughly maps to the percent deviation in the corresponding prediction of the rate constant, and the new method reduces the percent deviation of the latter from ∼40% to ∼20%. Similar reductions in the discrepancy between CASPT2 and MRCI+Q results recently were obtained for C3H3 + O2.54 For the VRC-TST calculations, the reference energy was computed using CASPT2(9,7)/cc-pVTZ. Several corrections to the reference energy were applied, including the corrections described by eq E3, as well as the more traditional corrections for basis set and relaxation energy, detailed in section 3.2.1. 3.1.2. TS 2: CH2CHOO → CH2CHO + O. The decomposition of vinylperoxy to CH2CHO + O is complicated by the presence of a low-lying excited state and by the presence of electronic degeneracies in O and resonances and near degeneracies in CH2CHO. CH2CHO (vinoxy) is resonantly stabilized, and the dominant resonance configuration has the radical located on the terminal carbon rather than the oxygen. The reverse process (addition) must break the resonance and confine the radical to F

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

atom) and the orbital degeneracy of the OH radical. The energies along the minimum energy pathway (MEP) were explored at the CASPT2 level with a (7,6) active space consisting of 5 orbitals that correlate with the π spaces of 3 CHCHO and OH, and the other radical orbital of 3CHCHO. The CASSCF part of the calculations include a state averaging over the two lowest states (correlating with the degenerate states of OH). The transition-state partition function for this reaction was obtained with variable reaction coordinate TST. The interaction energies for this analysis were directly determined at the CASPT2(7,6)/aug-cc-pVDZ level. The calculations incorporated a one-dimensional (in ROO) correction based on CASPT2(7,6)/cc-pVTZ evaluations of the relaxation energy and CASPT2(7,6)/cc-pV∞Z basis set corrections, obtained from extrapolation of QZ and 5Z evaluations. An infinite potential was assumed to separate the addition of the O atom to the O-side or the C-side of 3CHCHO, and only the addition to the O-side was considered. Two pivot points for the 3 CHCHO fragment were then placed above and below the O atom. For OH, the pivot point was placed at its center-of-mass. For large separations, a center-of-mass pivot point was also considered for the 3CHCHO fragment. The dominant resonance configuration of 3CHCHO has the π-radical on the C atom rather than the O atom. Thus, the reverse recombination process must first break the resonance in 3 CHCHO and localize the electron density on the less favorable site. This fact, coupled with the presence of strong hydrogen-bonding at long-range, results in the presence of a fairly well-defined saddle-point for the addition process. This saddle point, which lies at ROO = 2.1 Å and has a zero-point corrected energy of −1.7 kcal/mol relative to separated fragments, provides a fairly tight bottleneck to reaction. Correspondingly, the calculated reverse recombination rate, which is roughly constant at about 4 × 10−12 cm3 molecule−1 s−1 for temperatures of 500 K and higher, is considerably smaller than is typical for other radical−radical recombination reactions. The electronic degeneracy of OH, which contributes to the fragment partition function but not the transition-state one, also contributes to this small magnitude for the rate constant. 3.1.4. TS 10: c-CH(CH2)OO → c-CH(O)CH2O. Although the isomerization from c-CH(CH2)OO (dioxiranyl-methyl) to cCH(O)CH2O (oxiranyloxy) proceeds through a tight transition state, the CCSD(T)/cc-pVTZ geometry has a T1 diagnostic of 0.075, which is sufficiently large that single-reference methods are unsuitable. As noted in ref 32, accurate determination of this transition state is complicated by the presence of another, similar transition state, which corresponds to the flipping between enantiomers of oxiranyloxy. The active space used in refs 32 and 33 was (23,15), which corresponds to the Hartree− Fock valence space plus three unspecified antibonding molecular orbitals. On the basis of the present work, one can speculate that the three antibonding molecular orbitals probably would have been the O−O and 2 C−O σ* orbitals. Preliminary investigation into the appropriate active space began with the O−O σ, σ* and the carbon-centered radical. Although this (3,3) active space captures most of the multireference nature of the transition state, it is not sufficient. Further inspection revealed that the 2 C−O σ, σ* orbitals are important. Of still greater significance is the mixing between the lone pair of the oxyranyl oxygen and the O−O σ, σ* orbitals.

vinoxy lone pair (averaged over six states). The largest active space considered was (13,11), which includes the C−C σ, σ*, but this active space was computationally prohibitive for multistate calculations and reasonable basis sets. The results of the single reference versus multireference comparison are shown in Table 3. Table 3. Comparison between Coupled Cluster and CASPT2 Results for CH2CHOO → CH2CHO + O method ΔECBS + ΔET/Q ΔECC,total CASPT2(7,6) CASPT2(9,7) CASPT2(9,8) CASPT2(11,9) Mebel et al.:21 G2M(RCC,MP2)

barrier heighta (kcal/mol) −0.7 −0.5 −0.3 −0.9 −0.2 −0.2 1.7

excited state − ground state (kcal/mol)

1.5 1.7 1.5 1.7

a

Calculated relative to CH2CHO + O. Includes zero-point correction. In each case, the basis set was extrapolated to the complete basis set limit.

To facilitate a meaningful comparison between the coupledcluster and multireference calculations, the post-CCSD(T)/ccpVTZ corrections are separated into two rows. The first row includes only the basis-set extrapolation and the correction for higher-order excitations, because both CCSDT(Q) and CASPT2 are in some sense attempts to converge to the full CI. The ΔECBS + ΔET/Q barrier height is roughly in the middle of the CASPT2/cc-pV∞Z energies. The second row includes the three other corrections of eq E1, which increase the barrier height by 0.2 kcal/mol. The results in Table 3 demonstrate remarkably good agreement between the single-reference and multireference methods. On the basis of the present results, the zero-point corrected barrier height is assumed to be −0.5 kcal/mol, relative to CH2CHO + O. We estimate the corresponding 2σ uncertainty as ±0.5 kcal/mol. The CCSD(T)/cc-pVTZ saddle point was located at an O− O distance of 1.974 Å. The partition function for this transition state was computed as the sum of the contributions from both the ground and excited states, and the individual partition functions were computed using variational transition-state theory over an O−O bond distance of 1.70 to 2.05. A subset of these results is shown in the inset of Figure 3. The excited state contributed roughly 20% of the total flux for the C2H3 + O2 → CH2CHO + O rate constant. Although there must be some local minimum for r > 2.05, this portion of the potential was not determined, as it is irrelevant to the thermal kinetics for the high temperatures of interest here. 3.1.3. TS 9: CHCHOOH → 3CHCHO + OH. The cleavage of the O−O bond in 2-hydroperoxylvinyl (CHCHOOH) is one of the main product channels for HO2 + C2H2, which is an important reaction in acetylene combustion. As will be seen later, it also happens to be a significant channel for C2H3 + O2 at high temperatures. This O−O bond fission produces two radicals and thus its kinetic treatment requires the use of multireference wave functions. As with the CH2CHO + O pathway, this channel is further complicated by the presence of a resonance in the electronic structure for the 3CHCHO fragment (with a π-radical on either the terminal C or the O G

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

two C−O σ, σ* orbitals, the O−O σ, σ* orbitals, one oxygen lone-pair, and the carbon radical. Preliminary calculations confirmed that this channel has negligible influence on the final rate coefficients, so further exploration of the active space was not considered. 3.1.6. TS 20: C2H3 + O2 → C2H2 + HO2. The transition state for the direct abstraction reaction has a CCSD(T)/cc-pVTZ T1 diagnostic of 0.058. Preliminary calculations using CASPT2(9,7)/cc-pVTZ determined a barrier height of 10.3 kcal/mol. The corresponding rate constant, calculated using conventional transition-state theory, was 7% of the net addition rate constant at 2000 K and lower at lower temperatures. Therefore, further exploration of this transition state was deemed unnecessary. 3.2. RRKM/ME. 3.2.1. Total Rate Constant. The results for the C2H3 + O2 capture rate constant are shown in Figure 4. For

Thus, the transition-state geometry was optimized using CASPT2(9,8)/cc-pVTZ. Subsequent single-point calculations included the second lone pair and the C−C σ, σ* orbitals. Owing to the challenge of obtaining the same orbitals for the reactant and transition state, the smallest active space that would consistently encompass the essential orbitals was (13,11). These CASPT2(13,11) calculations were used to compute the barrier height relative to the reactants. Although the single-determinant wave function is ill-suited to describe the transition state, the same doublet/quartet splitting approach described in eq E2 in section 3.1.1 can be applied to this transition state. After proper rotation of the orbitals in the RHF calculation for the quartet surface, the CCSD(T)/ccpVTZ T1 diagnostic was 0.030 for the quartet. The doublet/ quartet splitting was computed using CASPT2(9,8), MRCI +Q(9,8), and CASPT2(13,11). The results are summarized in Table 4. Table 4. Barrier Height for the Transition State cCH(CH2)OO → c-CH(O)CH2O method

barrier heighta (kcal/mol)

CASPT2(13,11) CCSD(T) - 4−2CASPT2(9,8) 4 CCSD(T) - 4−2CASPT2(13,11) 4 CCSD(T) − 4−2MRCI+Q(9,8) 4 [ΔECBS + ΔET/Q] 4 [ΔECC, total] 4 [CASPT2(13,11)] Mebel et al.:21 G2M(RCC,MP2) Carpenter:32 CASPT2(23,15) Mebel and Kislov:33 CCSD(T) MRCI+Q

17.7 17.0 18.0 16.8 43.1 42.8 42.3 24.0 12.5

4

Figure 4. Computed capture rate constant for C2H3 + O2. Symbols correspond to experimental measurements, and the colored lines are the present work. Red corresponds to CASPT2 potentials, and blue to MRCI+Q. The dashed lines correspond to the interaction potential determined purely from the doublet surface, whereas the solid lines are the new method, based upon multireference evaluation of the doublet quartet splitting and CCSD(T) evaluation of the quartet energy, eq E3. For the finite pressures of the experimental data, back-dissociation is expected to be negligible, so the net rate of C2H3 disappearance can be correlated with the capture rate.

17.2 16.2

a

Calculated relative to c-CH(CH2)OO. Includes zero-point correction. In each case, the basis set was extrapolated to the complete basis set limit.

The two approaches in the present work are in remarkably good agreement, especially when compared to the previous attempts. Rows 5, 6, and 7 contain the difference in energy between the transition-state energy on the quartet surface and the reactant energy on the doublet surface. Comparing rows 5 and 7, we see that the T/Q corrected energy and the CASPT2 energy differ by 0.8 kcal/mol. Based upon these findings, the best estimate for the zero-point corrected barrier height is 17.3 kcal/mol. We estimate the 2σ uncertainty of these values as ±1 kcal/mol. 3.1.5. TS 11: c-CHCH2O2 → OCH2CHO. As with the isomerization from c-CH(CH2)OO to c-CH(O)CH2O, the transition state for the isomerization from c-CHCH2OO (dioxetanyl) to OCH2CHO (formyl-methoxy) has a high T1 diagnostic, 0.117 (Table 5). The transition-state geometry was optimized using CASPT2(9,8)/cc-pVTZ, which included the

the finite temperatures and pressures shown in the figure, the reverse rate is expected to be negligible, and the capture rate is an appropriate comparison for these conditions. The predicted rate constants are in excellent agreement with all the published data. In all calculations, the reference energy was the CASPT2(9,7)/cc-pVTZ. The dashed red line represents the inclusion of corrections for CASPT2(9,7)/cc-pVTZ geometry relaxation and CASPT2(13,11) energies extrapolated to the basis set limit. The dashed blue line represents the inclusion of corrections for MRCI+Q(9,7)/cc-pVTZ geometry relaxation and MRCI+Q(11,9) energies extrapolated to the basis set limit. Of particular significance is the effect of the doublet/quartet splitting method, eq E3, which is represented by the two solid lines. This new approach reduces the percent deviation from ∼50% to ∼25%, or a factor of 2 improvement. For all subsequent rate coefficient calculations, the interaction potential used in the VRC-TST calculations was the CASPT2(13,11)/cc-pV∞Z with the doublet/quartet correction. 3.2.2. Product Formation. The rate constants for the stabilization of the initial adduct, CH2CHOO, and the eight bimolecular product channels at 0.01, 1, and 100 atm are shown in Figure 5. At 0.01 atm, the rate of CH2CHOO stabilization is

Table 5. Barrier Height for the Ring-Opening Transition state c-CHCH2OO → OCH2CHO

a

method

barrier heighta (kcal/mol)

CASPT2(9,8) MRCI+Q(9,8)

1.4 2.5

Calculated relative to c-CHCH2OO. Includes zero-point correction. H

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 5. Arrhenius plots for the individual product channels at three different pressures for the addition of C2H3 + O2.

were taken for vibrationally excited C2H3 in pure O2, so a direct comparison is not possible. This discrepancy suggests there may be some role for nonstatistical effects in the product distributions, which would not be too surprising given the magnitude of the microcanonical rates (e.g., 1 × 1012 s−1) for the right-hand side of the PES. The dominant product, CH2O + HCO, is formed from a simple β-scission reaction with a fairly loose transition state, and nonequilibrated energy distributions at high energies would still favor this pathway. As a first-order approximation to gauge the impact of nonstatistical behavior, we modified PAPER to impose an upper limit of 1 × 1012 s−1 for all microcanonical rate coefficients. The effect was roughly 10% for the four exothermic bimolecular product channels. Overall, at temperatures between 1000 and 2000 K, all eight bimolecular product channels are formed at an appreciable rate, including the products observed by Wang and co-workers,25,26 with approximately 2 orders of magnitude difference between the largest and smallest rate constants. At pressures above 10 atm, CH2CHOO is also formed at an appreciable rate. A direct comparison of the predicted pressure dependence of CH2CHOO stabilization at room temperature with the Chishima et al. data is not possible. If we assume that the effects of pressure broadening on the absorption cross sections for CH2CHOO and C2H3 roughly cancel, then we can compare the relative increase in CH2CHOO yield as a function of pressure. Between 50 and 120 Torr, Chishima et al. suggest an approximately 3-fold increase in CH2CHOO stabilization,28 whereas the present work predicts a 45-fold increase (i.e., a branching fraction increase from 8.5 × 10−6 to 3.9 × 10−4). However, the predicted rate of collisional stabilization under these low stabilization probability conditions will be highly sensitive to the details of the energy transfer model. The Chishima data could suggest the possibility of supercollisions, which are not included in the present work. These supercollisions will be less significant at higher pressures, where the branching fraction is significant. It should be mentioned that Matsugi and Miyoshi were unable to reproduce the CH2CHOO absorption signal reported by Chishima et al.28 They attribute the signal at 446.4 nm to interference from side products, not CH2CHOO. There are two distinct pathways to form CH2CO (ketene) + OH, but only the lower energy one is considered in the present analysis, so the rate constant for the formation of CH2CO + OH presented here is a lower bound. The pathway considered

negligible; indeed, above 1300 K, CH2CHOO ceases to exist as a distinct species. By 1 atm, CH2CHOO stabilization is a minor channel at low temperatures, accounting for 1−10% of the total branching ratio between 300 and 1000 K; it ceases to exist above 1600 K. At 100 atm, however, CH2CHOO stabilization is the dominant product channel until 1500 K; it persists as chemically distinct until 2100 K. The dynamics through the lower-energy portion of the potential energy surface (i.e., the four intermediates without an O−O bond that lead to the formation of OCHCHO + H, CH2O + HCO, CO + CH3O, and CO2 + CH3) is effectively collisionless. Due to the large exothermicity relative to reactants, no stabilization of these four intermediates occurs, even at the highest pressures considered. The relative branching ratio among these four bimolecular product channels (e.g., the ratio of OCHCHO + OH to CH2O + HCO) is invariant with pressure. Indeed, separate calculations performed using a 2D master equation in the low-pressure limit confirmed the same relative branching ratios for these channels. This finding also demonstrates the adequacy of the 1D master equation for these branching ratios. As expected, CH2O + HCO and CH2CHO + O are the dominant bimolecular product channel at lower and higher temperatures, respectively. At room temperature and 20 Torr, the branching fraction toward the combined CH2O + HCO and CH2O + H + CO channel is 0.79. Given the assumption that 70% of the nascent HCO goes directly to H + CO,30 the predicted HCO yield is 0.24, which is in excellent agreement with the measured result of Matsugi and Miyoshi, 0.22 ± 0.03. The branching fraction for CH2CHO + O at room temperature and pressures between 10 and 100 Torr is 10%, which is within the stated uncertainties of the Oguchi data.27 Perhaps surprisingly, the next most significant bimolecular channel is OCHCHO (glyoxal) + H, followed by 3CHCHO + OH. Although the previous mechanisms assumed that C2H2 + HO2 is the third most significant channel, the present work suggests that this channel is minor until all but the highest temperatures. Inclusion of the direct abstraction pathway, TS 20, will increase the yield somewhat at the highest temperatures, but its contribution is otherwise negligible. The two most exothermic products, CO + CH3O, and CO2 + CH3 are between 0.1% and 1% of the net rate constant. If we assume that 70% of the HCO promptly dissociates, then our present results have a CO:CO2 ratio of 95:1, or roughly a factor of 5 greater than the estimate from ref 26. However, it should be noted that the data in ref 26 I

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 6. Branching fractions for CH2CHOO. The dashed line is for isomerization to the three-membered ring, c-CH(CH2)OO; the solid lines correspond to the different bimolecular channels.

Figure 7. Branching fractions for the nine products of C2H3 + O2. The two left panes are at 30 Torr, and the two right panes are at 20 atm. The top and bottom panes are the same data but with linear and logarithmic y-axes, respectively, to highlight the low branching channels.

is a four-membered ring internal H-abstraction, CH2CHOO → CH2COOH, TS 4, represented as a solid green line in Figure 1. The resulting α-hydroperoxyl alkenyl radical is not stable and promptly dissociates to CH2CO + OH. The barrier height for TS 4 is 6.5 kcal/mol above C2H3 + O2. The second pathway is a roaming pathway from the longrange CH2CHO···O complex that is part of the CH2CHOO → CH2CHO + O channel, represented as a dashed green line in Figure 1. Because the saddle point for CH2CHOO →

CH2CHO···O is 7.5 kcal/mol below C2H3 + O2 and 0.5 kcal/mol below CH2CHO + O, the vast majority of flux that clears the saddle point will simply continue on to CH2CHO + O instead of roaming to CH2CO + OH.56 The roaming pathway is expected to be no more than 2−5% of the CH2CHOO → CH2CHO + O pathway. Although roaming will have a negligible effect on the rate constant for C2H3 + O2 → CH2CHO + O, it could contribute to an increase in the rate constant for C2H3 + O2 → CH2CO + OH. Nonetheless, a J

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 8. Temperature at which CH2CHO + O becomes the dominant product channel. The thick solid line is for the nominal PES values, and the thin dash-dot line corresponds to TS 2 ± 0.5 and TS 10 ∓ 1.0 kcal/mol, as in Figure 7.

the uncertainty in the present methodology, the barrier heights were adjusted up and down by the 2σ uncertainty estimates. These results are shown in Figure 7. The two pressures, 0.04 atm (or 30 Torr) and 20 atm, were chosen to represent the upper and lower values that are commonly used in most combustion experiments, such as laser flash photolysis reactors and laminar flames at low pressures and flow tubes and shock tubes at high pressures. The linear scale (Figure 7a,b) illustrates the effect of the individual barrier heights on the main products. Figure 7c,d highlights the complex nature of the change in the product distribution. The temperature at which CH2CHO + O overtakes CH2O + HCO as the dominant product channel is 1400 K, as seen in Figure 8. In contrast, this temperature was 800 K for Mebel et al.,21 2900 K for Carriere et al.,31 and 2000 K for Lopez et al.23 Of perhaps greater kinetic significance is the temperature at which CH2CHO + O is greater than the sum of all other channels, which is shown in red. Above 1700 K, C2H3 + O2 is a predominantly chain branching reaction. Above 10 atm, the collisional stabilization of CH2CHOO competes with both CH2O + HCO and CH2CHO + O. The temperature at which CH2CHO + O overtakes all other channels increases with increasing pressure. Consequently, for a fixed temperature, the amount of chain branching decreases with increasing pressure. This decrease in chain branching is amplified by the fact that CH2CHOO decomposes primarily to CH2O + HCO and not CH2CHO + O (Figure 6). Thus, the more CH2CHOO is stabilized, the less chain branching there will be. The uncertainties in the two key transition states, TS 2 and TS 10, have a dramatic impact on the switching temperature. On the basis of the variation in these two barrier heights, the temperature at which CH2CHO + O overtakes all other channels can range from 1300 to 2200 K. The error bars in Figure 8 are based upon the 2σ uncertainties in TS 2 and TS 10, as depicted in Figure 7. The uncertainty in the temperature at which CH2CHO + O overtakes all other channels decreases at higher pressures, from ±400 K at 1 atm to ±100 K at 100 atm. For the temperature at which CH2CHO + O overtakes just CH2O + HCO, the opposite is true. This behavior can be explained by Boltzmann factors.

roaming analysis was not considered to be worth the effort. In the current analysis, this pathway is 0.1% at 1500 K and 0.6% at 2000 K. The complete set of pressure-dependent rate constants are provided in the Supporting Information. As discussed in the Introduction, the HCO that is formed with CH2O will be vibrationally excited, and most of this HCO will dissociate on a time scale faster than collisional stabilization. In the accompanying rate constants, the C2H3 + O2 → CH2O + HCO rate coefficients are split between two channels, with 30% going to CH2O + HCO and 70% going to CH2O + H + CO. The same partitioning is used for CH2CHOO decomposition as well. 3.2.3. CH2CHOO Decomposition. The branching fractions for the thermal decomposition of vinylperoxy at 1, 10, and 100 atm are shown in Figure 6. Because the rate of collisional stabilization is negligible below 1 atm, those values are not reported. At low temperatures, most CH2CHOO will isomerize to the neighboring three-member ring, c-CH(CH2)OO. This isomer, in turn, either decomposes to CH2O + HCO or isomerizes back to CH2CHOO. In contrast to the chemically activated C2H3 + O2 reaction, CH2O + HCO is the most significant bimolecular product channel at all temperatures. The next most significant channels are CH2CHO + O, C2H3 + O2, and OCHCHO + H, with the relative ranking a strong function of temperature and pressure. Because the CH2CHO + O product channel is never more than ∼10% of the total decomposition rate, the collisional stabilization of CH2CHOO shifts the ultimate product distribution away from CH2CHO + O and toward CH2O + HCO. Note that c-CH(CH2)OO ceases to be a distinct chemical species at 800 K and 1 atm, 900 K at 10 atm, and 1100 K at 100 atm, which is why the dashed black line terminates before the solid (bimolecular) lines in Figure 6. Similarly, vinylperoxy itself is no longer chemically stable at 1600 K and 1 atm, 1800 K and 10 atm, and 2100 K and 100 atm. 3.2.4. Sensitivity Analysis for TS 2 and TS 10. As mentioned in the Introduction, the transition states CH2CHOO → CH2CHO + O (TS 2) and c-CH(CH2)OO → c-CH(O)CH2O (TS 10) are responsible for the branching between the main product channels. Unfortunately, these two transition states are among the most difficult to characterize. To test the impact of K

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A

Figure 9. Comparison of the present work with four prior mechanisms at four different pressures. The rate constants from (a) Carriere et al.31 and Lopez et al.,23 (b) Mebel et al.,21 Lopez et al.,23 and Matsugi et al.,29 (c) Carriere et al.,31 and (d) Lopez et al.23 and Matsugi et al.29 To simplify the comparison, the HCO and H + CO channels for the present work and for Matsugi have been lumped into a single channel (shown in black).

3.2.5. Comparison with Previous Values. The four previously published rate constants for this system are Mebel et al.,21 Carriere et al.,31 Lopez et al.,23 and Matsugi and Miyoshi.29 These rate constants are shown in Figure 9, along with the present results. In most instances, the previous work is substantially different from the present work, with Matsugi and Miyoshi demonstrating the best agreement and Carriere et al. the worst. Although Matsugi and Miyoshi are in close agreement with the present work, they include only four product channels: CH2CHOO, CH2CHO + O, CH2O + HCO, and CH2O + H + CO (the latter two being lumped into a single channel in Figure 9). On the basis of the product distribution, it is expected that the present work will accelerate the rate of fuel consumption in combustion models, at least at lower pressures. First, the new results suggest that the rate constant for CH2CHOO → CH2CHO + O is larger than most previous estimates, and this pathway is chain branching. Second, whereas previous models suggested that C2H2 + HO2 was the third most significant product channel, the present work predicts that OCHCHO + H and 3CHCHO + OH have larger rate constants, and H and OH are more aggressive than HO2 at fuel oxidation. By increasing the relative yield of O, H, and OH at the expense of HCO and HO2, we expect the products of C2H3 + O2 to increase the net rate of oxidation.

4. CONCLUSIONS The most accurate analysis to date of the C2H3O2 potential energy surface and kinetics is presented. A new method is described for computing the interaction potential for R + O2 reactions. The method, which combines accurate determination of the potential on the quartet surface with multireference calculations of the doublet/quartet splitting, decreases the uncertainty in the potential and thence the rate constants by more than a factor of 2. The main bimolecular product channels are CH2O + HCO at lower temperatures and CH2CO + O at higher temperatures. Above 10 atm, the collisional stabilization of CH2CHOO directly competes with these two product channels. CH2CHOO decomposes primarily to CH2O + HCO. C2H3 + O2 will be predominantly chain branching above 1700 K. The next two most significant bimolecular products are OCHCHO + H and 3CHCHO + OH, and not C2H2 + HO2. The two most important transition states remain a challenge for electronic structure methods. The uncertainties in these two barrier heights result in a significant uncertainty in the temperature at which OCH2CHO + O overtakes all other product channels. L

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A



(13) Cooke, D. F.; Williams, A. Shock-Tube Studies of the Ignition and Combustion of Ethane and Slightly Rich Methane Mixtures with Oxygen. Symp. Combust., [Proc.] 1971, 13, 757−766. (14) Baldwin, R. R.; Walker, R. W. Elementary Reactions in the Oxidation of Alkenes. Symp. Combust., [Proc.] 1981, 18, 819−829. (15) Park, J.-Y.; Heaven, M. C.; Gutman, D. Kinetics and Mechanism of the Reaction of Vinyl Radical with Molecular-Oxygen. Chem. Phys. Lett. 1984, 104, 469−474. (16) Slagle, I. R.; Park, J.-Y.; Heaven, M. C.; Gutman, D. Kinetics of Polyatomci Free-Radicals Produced by Laser Photolysis. 3. Reaction of Vinyl Radicals with Molecular Oxygen. J. Am. Chem. Soc. 1984, 106, 4356−4361. (17) Westmoreland, P. R. Thermochemistry and Kinetics of C2H3 + O2 Reactions. Combust. Sci. Technol. 1992, 82, 151−168. (18) Bozzelli, J. W.; Dean, A. M. Hydrocarbon Radical Reactions with O2 − Comparison of Allyl, Formyl, and Vinyl to Ethyl. J. Phys. Chem. 1993, 97, 4427−4441. (19) Carpenter, B. K. Computational Prediction of New Mechanisms for the Reactions of Vinyl and Phenyl Radicals with MolecularOxygen. J. Am. Chem. Soc. 1993, 115, 9806−9807. (20) Carpenter, B. K. Ab Initio Computation of Combustion Kinetics. 1. Vinyl Radical + O2. J. Phys. Chem. 1995, 99, 9801−9810. (21) Mebel, A. M.; Diau, E. W. G.; Lin, M. C.; Morokuma, K. Ab Initio and RRKM Calculations for Multichannel Rate Constants of the C2H3 + O2 Reaction. J. Am. Chem. Soc. 1996, 118, 9759−9771. (22) Chang, A. Y.; Bozzelli, J. W.; Dean, A. M. Kinetic Analysis of Complex Chemical Activation and Unimolecular Dissociation Reactions Using QRRK Theory and the Modified Strong Collision Approximation. Z. Phys. Chem. 2000, 214, 1533−1568. (23) Lopez, J. G.; Rasmussen, C. L.; Alzueta, M. U.; Gao, Y.; Marshall, P.; Glarborg, P. Experimental and Kinetic Modeling Study of C2H4 Oxidation at High Pressure. Proc. Combust. Inst. 2009, 32, 367− 375. (24) Eskola, A. J.; Timonen, R. S. Kinetics of the Reactions of Vinyl Radicals with Molecular Oxygen and Chlorine at Temperatures 200− 362 K. Phys. Chem. Chem. Phys. 2003, 5, 2557−2561. (25) Wang, H.; Wang, B.; He, Y.; Kong, F. The Gaseous Reaction of Vinyl Radical with Oxygen. J. Chem. Phys. 2001, 115, 1742−1746. (26) Feng, W.; Wang, B. The Reaction of Vinyl Radical with O2 Studied by Time-Resolved Infrared Emission Spectroscopy. Chem. Phys. Lett. 2002, 356, 505−510. (27) Oguchi, T.; Sato, Y.; Matsui, H. The CH2CHO + O Channel of the Reaction of Vinyl Radical with O2. Chem. Phys. Lett. 2009, 472, 181−184. (28) Chishima, H.; Koshi, M.; Tonokura, K. Pressure Dependence of Vinylperoxyl Radical Formation in the Reaction of Vinyl Radical with Molecular Oxygen. Chem. Lett. 2009, 38, 1150−1151. (29) Matsugi, A.; Miyoshi, A. Yield of Formyl Radical from the Vinyl + O2 Reaction. Int. J. Chem. Kinet. 2014, 46, 260−274. (30) Klippenstein, S. J.; Georgievskii, Y.; Miller, J. A.; Nummela, J. A.; Carpenter, B. K.; Westmoreland, P. R. Vinyl + O2: A Complete Theoretical Treatment. 3rd Joint Meeting of the U.S Sections of The Combustion Institute, Chicago, IL, 2003. (31) Carriere, T.; Westmoreland, P. R.; Kazakov, A.; Stein, Y. S.; Dryer, F. L. Modeling Ethylene Combustion from Low to High Pressure. Proc. Combust. Inst. 2002, 29, 1257−1266. (32) Carpenter, B. K. Ring Opening of Dioxiranylmethyl Radical: A Caution on the Use of G2-Type Ab Initio MO Methods for Mechanistic Analysis. J. Phys. Chem. A 2001, 105, 4585−4588. (33) Mebel, A. M.; Kislov, V. V. The C2H3 + O2 Reaction Revisited: Is Multireference Treatment of the Wave Function Really Critical. J. Phys. Chem. A 2005, 109, 6993−6997. (34) Tajti, A.; Szalay, P. G.; Csaszar, A. G.; Kallay, M.; Gauss, J.; Valeev, E. F.; Flowers, B. A.; Vazquez, J.; Stanton, J. F. HEAT: High Accuracy Extrapolated Ab Initio Thermochemistry. J. Chem. Phys. 2004, 121, 11599−11613. (35) Karton, A.; Rabinovich, E.; Martin, J. M. L.; Ruscic, B. W4 Theory for Computational Thermochemistry: In Pursuit of Confident Sub-kJ/mol Predictions. J. Chem. Phys. 2006, 125, 144108.

ASSOCIATED CONTENT

S Supporting Information *

Cartesian coordinates are provided for all the stationary points in the potential energy surface. The temperature- and pressuredependent rate coefficients are provided in a CHEMKINcompatible PLOG format. The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpca.5b01088.



AUTHOR INFORMATION

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Basic Energy Sciences, Division of Chemical Sciences, Geosciences, and Biosciences under Contract No. DE-AC02-06CH11357 as part of the Argonne-Sandia Consortium on High-Pressure Combustion Chemistry (FWP # 59044). C.F.G. gratefully acknowledges financial support from the Argonne Director’s Postdoctoral Fellowship and from Brown University.



REFERENCES

(1) Frenklach, M. Reaction Mechanism of Soot Formation in Flames. Phys. Chem. Chem. Phys. 2002, 4, 2028−2037. (2) Richter, H.; Howard, J. B. Formation and Consumption of Single-Ring Aromatic Hydrocarbons and Their Precursors in Premixed Acetylene, Ethylene, and Benzene Flames. Phys. Chem. Chem. Phys. 2002, 4, 2038−2055. (3) Knyazev, V. D.; Slagle, I. R. Experimental and Theoretical Study of the C2H3 = H + C2H2 Reaction. Tunneling and the Shape of Falloff Curves. J. Phys. Chem. 1996, 100, 16899−16911. (4) Miller, J. A.; Klippenstein, S. J. The H + C2H2 (+M) = C2H3 (+M) and H + C2H2 (+M) = C2H5 (+M) Reactions: Electronic Structure, Variational Transition State Theory, and Solutions to a Two-Dimensional Master Equation. Phys. Chem. Chem. Phys. 2004, 6, 1192−1202. (5) Jasper, A. W.; Pelzer, K. M.; Miller, J. A.; Kamarchik, E.; Harding, L. B.; Klippenstein, S. J. Predictive A Priori Pressure-Dependent Kinetics. Science 2014, 346, 1212−1215. (6) Ismail, H.; Goldsmith, C. F.; Abel, P. R.; Howe, P. T.; Fahr, A.; Halpern, J. B.; Jusinski, L. E.; Georgievskii, Y.; Taatjes, C. A.; Green, W. H. Pressure and Temperature Dependence of the Reaction of Vinyl Radical with Ethylene. J. Phys. Chem. A 2007, 111, 6843−6851. (7) Goldsmith, C. F.; Ismail, H.; Abel, P. R.; Green, W. H. Pressure and Temperature Dependence of the Reaction of Vinyl Radical with Alkenes II: Measured Rates and Predicted Product Distributions for Vinyl + Propene. Proc. Combust. Inst. 2009, 32, 139−148. (8) Goldsmith, C. F.; Ismail, H.; Green, W. H. Pressure and Temperature Dependence of the Reaction of Vinyl Radical with Alkenes III: Measured Rates and Predicted Product Distributions for Vinyl + Butene. J. Phys. Chem. A 2009, 113, 13357−13371. (9) Klippenstein, S. J.; Miller, J. A.; Harding, L. B. Resolving the Mystery of Prompt CO2: The HCCO + O2 Reaction. Proc. Combust. Inst. 2002, 29, 1209−1217. (10) Krueger, H.; Weitz, E. Diode-Laser Probes of Vinyl Radical Kinetics − The Reaction of C2H3 with HCl and DCl. J. Chem. Phys. 1988, 88, 1608−1616. (11) Fahr, A.; Laufer, A. H. Ultraviolet-Absorption of the Vinyl Radial and Reaction with Oxygen. J. Phys. Chem. 1988, 92, 7229− 7232. (12) Knyazev, V. D.; Slagle, I. R. Kinetics of the Reaction of Vinyl Radical with Molecular Oxygen. J. Phys. Chem. 1995, 99, 2247−2249. M

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry A (36) Martin, J. M. L. Ab Initio Total Atomization Energies of Small Molecules − Towards the Basis Set Limit. Chem. Phys. Lett. 1996, 259, 669−678. (37) Feller, D.; Dixon, D. A. Extended Benchmark Studies of Coupled Cluster Theory Through Triple Excitations. J. Chem. Phys. 2001, 115, 3484−3496. (38) Wolf, A.; Reiher, M.; Hess; Bernd, A. The generalized DouglasKroll transformation. J. Chem. Phys. 2002, 117, 9215−9226. (39) Barone, V. Anharmonic vibrational properties by a fully automated second-order perturbative approach. J. Chem. Phys. 2005, 122, 014108. (40) Werner, H.-J.; et al. MOLPRO, version 2012.1; see http://www. molpro.net. (41) Frisch, M. J.; et al. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (42) Klippenstein, S. J. Variational Optimizations in the RiceRamsberger-Kassel-Marcus Theory Calculations for Unimolecular Dissociations with No Reverse Barrier. J. Chem. Phys. 1992, 96, 367−371. (43) Georgievskii, Y.; Klippenstein, S. J. Transition State Theory for Multichannel Addition Reactions: Multifaceted Dividing Surfaces. J. Phys. Chem. A 2003, 107, 9776−9781. (44) Truhlar, D. G.; Klippenstein, S. J. Current Status of TransitionState Theory. J. Phys. Chem. 1996, 100, 12771−12800. (45) Klippenstein, S. J. RRKM Theory and its Implementations. Comprehensive Chemical Kinetics; Elsevier: Amsterdam, 2003; Vol. 39, Chapter 2, 55−103. (46) Georgievskii, Y.; Miller, J. A.; Burke, M. P.; Klippenstein, S. J. Reformulation and Solution of the Master Equation for Multiple-Well Chemical Reactions. J. Phys. Chem. A 2013, 117, 12146−12154. (47) Georgievskii, Y.; Jasper, A. W.; Zádor, J.; Miller, J. A.; Burke, M. P.; Goldsmith, C. F.; Klippenstein, S. J. PAPER v1 (2014). (48) Goldsmith, C. F.; Green, W. H.; Klippenstein, S. J. Role of O2 + QOOH in Low-Temperature Ignition of Propane. 1. Temperature and Pressure Dependent Rate Coefficients. J. Phys. Chem. A 2012, 116, 3325−3346. (49) Goldsmith, C. F.; Tomlin, A. S.; Klippenstein, S. J. Uncertainty Propagation in the Derivation of Phenomenological Rate Coefficients from Theory: A Case Study of n-Propyl Radical Oxidation. Proc. Combust. Inst. 2013, 34, 177−185. (50) Hippler, H.; Troe, J.; Wendelken, H. J. Collisional Deactivation of Vibrationally Highly Excited Polyatomic Molecules 2. Direct Observations for Excited Toluene. J. Chem. Phys. 1983, 78, 6709− 6717. (51) Lee, T. J.; Taylor, P. R. A Diagnostic for Determining the Quality of Single-Reference Electron Correlation Methods. Int. J. Quantum Chem. 1989, 23, 199−207. (52) Ghigo, G.; Roos, B. O.; Malmquvist, P.-Å. A Modified Definition of the Zeroth-Order Hamiltonian in Multiconfigurational Perturbation Theory (CASPT2). Chem. Phys. Lett. 2004, 396, 142− 149. (53) Harding, L. B.; Klippenstein, S. J.; Lischka, H.; Shepard, R. Comparison of Multireference Configuration Interaction Potential Energy Surfaces for H + O2 → HO2: The Effect of Internal Contraction. Theor. Chem. Acc. 2014, 133, 1429−1436. (54) Moradi, C. P.; Morrison, A. M.; Klippenstein, S. J.; Goldsmith, C. F.; Douberly, G. E. Propargyl + O2 Reaction in Helium Droplets: Entrance Channel Barrier or Not? J. Phys. Chem. A 2013, 117, 13626− 13635. (55) Deskevich, M. P.; Nesbitt, D. J.; Werner, H. J. Dynamically Weighted Multiconfigurational Self-Consistent Field: Multistate Calculations for F + H2O → HF + OH Reaction Paths. J. Chem. Phys. 2004, 120, 7281−7289. (56) Klippenstein, S. J.; Georgievskii, Y.; Harding, L. B. Statistical Theory for the Kinetics and Dynamics of Roaming Reactions. J. Phys. Chem. A 2011, 115, 14370−14381.

N

DOI: 10.1021/acs.jpca.5b01088 J. Phys. Chem. A XXXX, XXX, XXX−XXX