Pressure-dependent kinetics of peroxy radicals formed in isobutanol combustion

, Bio-derived isobutanol has been approved as a gasoline additive in the US, but our understanding of its combustion chemistry still has significant uncertainties. Detailed quantum calculations could improve model accuracy leading to better estimation of isobutanol’s combustion properties and its environmental impacts. This work examines 47 molecules and 38 reactions involved in the first oxygen addition to isobutanol’s three alkyl radicals located a , b , and g to the hydroxide. Quantum calculations are mostly done at CCSD(T)-F12/cc-pVTZ-F12//B3LYP/CBSB7, with 1-D hindered rotor corrections obtained at B3LYP/6-31G(d). The resulting potential energy surfaces are the most comprehensive isobutanol peroxy networks published to date. Canonical transition state theory and a 1-D microcanonical master equation are used to derive high-pressure-limit and pressure-dependent rate coeﬃcients, respectively. At all conditions studied, the recombination of g -isobutanol radical with O 2 forms HO 2 + isobutanal. The recombination of b -isobutanol radical with O 2 forms a stabilized hydroperoxy alkyl radical below 400 K, water + an alkoxy radical at higher temperatures, and HO 2 + an alkene above 1200 K. The recombination of b -isobutanol radical with O 2 results in a mixture of products between 700–1100 K, forming acetone + formaldehyde + OH at lower temperatures and forming HO 2 + alkenes at higher temperatures. The barrier heights, high-pressure-limit rates, and pressure-dependent kinetics generally agree with the results from previous quantum chemistry calculations. Six reaction rates in this work deviate by over three orders of magnitude from kinetics in detailed models of isobutanol combustion, suggesting the rates calculated here can help improve modeling of isobutanol combustion and its environmental fate.


Introduction
In June 2018, the U.S. Environmental Protection Agency approved the large scale blending of bio-derived isobutanol into gasoline at concentrations up to 16%. 1,2 This newly approved oxygenate has multiple advantages over ethanol, such as lower volatility, 3 lower hydroscopicity, 4 and higher energy density, 3 while providing similar knock resistance. 5 Relative to larger alcohol additives, isobutanol also shows longer ignition delay times with lower CO and NO x emissions, 6 indicating it could help boost the octane numbers of fuel while lowering other emissions. With this approval, isobutanol could become a major component of transportation fuel in the United States.
While extensive experimental testing was necessary prior to the approval, 4 detailed isobutanol combustion models are still inaccurate at engine-relevant conditions, indicating a gap in our understanding of isobutanol combustion. Detailed isobutanol mechanisms typically deviate more from experiments than do models of n-butanol, tert-butanol, and sec-butanol. 7,8 Multiple inter-model comparisons using detailed isobutanol models have showed that at lean, low-temperature conditions the models yield a wide range of ignition timings that deviate from experimental ignition delay measurements in homogeneous charge compression 9 and shock tube setups. 10 One mechanism shows the largest deviations from experiments below 850 K and at fuel lean conditions, indicating that the dependence on O 2 might be incorrectly described. 11 Such shortcomings indicate a need for improving our understanding of isobutanol combustion in combustion regimes, especially where peroxy radical chemistry is important. 12 The dominant degradation pathway of peroxy radicals at various temperatures and pressures impacts overall ignition properties such as the negative temperature coefficient region, in which ignition slows with increasing temperature.
While several aspects of isobutanol combustion chemistry have been studied in detail, such as the initial hydrogen abstraction by OH 13 and the oxidation of its smaller oxidation products, 14 accurate pressure-dependent kinetics of peroxy radicals formed in isobutanol combustion has not been fully developed nor incorporated into detailed models. This is in spite of some isobutanol studies showing that ignition delay is sensitive to peroxy radical reactions. 8,10 Given the sensitivity of this fuel to equivalence ratio and pressure, 5 the formation and fate of peroxy radicals can significantly impact ignition timing and multi-stage heat release. 12,15 Peroxy radical chemistry is also expected to control the environmental fate of isobutanol emitted into the environment because peroxy radicals are also key intermediates in atmospheric degradation processes. Fig. 1 shows the peroxy radicals studied in this work, which correspond to the three predominant hydrogen abstractions that isobutanol can undergo. The alkoxy radical formed from hydrogen abstraction of isobutanol by OH is expected to be much less common 13 and thus not included in this work. For each of the three isomers studied, various unimolecular reactions can occur, as shown in Fig. 2. The peroxy radicals (RO 2 ) can isomerize through intramolecular hydrogen abstractions to form various alkyl hydroperoxide radicals (QOOH), or can unimolecularly decompose to form HO 2 and an alcohol with a double bond. QOOH can decompose in a variety of ways. Some decomposition pathways are highly structure-specific requiring additional functional groups, such as the pathway forming H 2 O and an ketoalkoxy radical. 16 b-scission pathways are more general and can lead QOOH to form the same products as HO 2 elimination (if a C-O bond breaks) or an alkyl radical and double bond (if a C-C bond breaks). 17 Multiple simultaneous b-scissions are also possible, forming three or more fragments. Cyclic ethers and OH can also result from QOOH decomposition. Like the alkyl radicals, QOOH can also undergo O 2 addition, which is important for fast, lowtemperature combustion. However, these pathways were not examined because calculations in this work suggest stabilized QOOH does not appear in appreciable quantities above 500 K.
Previous quantum calculations of these important radicals have mostly focused on bRO 2 since it is also formed in the OH-initiated oxidation of isobutene. Sun et al. and Lizardo-Huerta et al. report high-pressure-limit rates obtained from CBSQ//B3LYP/6-31G(d,p) and CBS-QB3, respectively. 18,19 Similar work has not been conducted for the other two isomers, aRO 2 and gRO 2 , but studies on similar species have been done. Zádor et al. have studied the fate of the a-peroxy radical derived from ethanol, though the full surface of aRO 2 derived from isobutanol has not, to our knowledge, been evaluated. 20 Welz et al. calculated a potential energy surface for the gRO 2 isomer, but did not include reactions of all isomers nor report reaction rates. 16 Given the errors for isobutanol models at low temperature and the lack of quantum calculations for the peroxy intermediates, we present thermodynamic and kinetic parameters for isomers, reactants, and products on the a-, b-, and g-isobutanol peroxy surfaces. Since there is potential for some for these reactions to be dependent on pressure, 21 both pressure-dependent and high-pressure-limit rate coefficients are reported. Comparison to other work and sensitivity analysis are conducted to determine the robustness of the rates.

Methods
To quantify the rates of reaction of isobutanol peroxy radicals, we conducted quantum calculations of the lowest-energy stable conformers and transition states and used those calculations in deriving rate coefficients with both transition state theory (TST) and the 1-D master equation.

Quantum calculations
To efficiently obtain accurate rate coefficients necessary for our analysis, we conducted geometry optimization, hindered rotor scans, and frequency calculations with B3LYP/CBSB7 or BMK/ 6-31g(d) using Gaussian 03. 22 Energy was calculated with CCSD(T)-F12/cc-pVTZ-F12 using either Molpro 2012.1.21 or Molpro 2015.1.37, including scaling the triples energy to the ratio of correlation energies from F12 and MP2. 23 For each molecule and transition state, we first conducted a geometry optimization of an initial guess. Transition states were verified to have one imaginary frequency. Internal Reaction Coordinate (IRC) calculations were conducted at B3LYP/CBSB7 level theory and visualized to ensure Fig. 1 The three peroxy isomers derived from isobutanol correspond to the three main sites for hydrogen abstraction from isobutanol. the transition states correspond with the correct reactants and products (see ESI, † for trajectories of each IRC calculation). For each rotatable dihedral, we calculated relaxed hindered rotor scans using 10-degree increments. When hindered rotor scans ended in a different conformer than the lowest-energy conformer, hindered rotor scans were rerun with the fixed dihedral angles and/or bond lengths to ensure an the scan proceeded as intended. While most geometry, hindered rotor and frequency calculations used the B3LYP method, calculations for compounds with a hydroperoxy group on the a carbon used BMK because B3LYP proved difficult to converge on these structures.
At each stationary point geometry, single-point CCSD(T)-F12a/cc-pVTZ-F12 were performed. Quantum chemistry outputs necessary to reproduce the thermodynamic and kinetic calculations are available in the ESI. † Of the 85 calculations performed, over a third had T 1 diagnostic values above 0.02, indicating that many of these structures likely have multireference character. 24 Three of these had T 1 diagnostic values above 0.05, corresponding to the formation of alkoxy radicals from RO 2 . The D 1 values also show similar effects with over a third of structures having a value greater than 0.1. The potential for multireference behavior, as indicated by these diagnostics, is thought to increase the uncertainties in the final energies to approximately 21 kJ mol À1 . 20 Based on work that compared single reference and multireference energies for a very similar system (oxygen addition to hydroxyethyl radical), 20 the difference in energies did not have any convincing systematic trend.

Thermodynamic and kinetic calculations
For thermodynamic and kinetic parameters, this work used Arkane, which is freely available as part of the Reaction Mechanism Generator software suite. 25 The partition function for each molecule and transition state was obtained using the rigid-rotor harmonic oscillator approximation with independent 1-D hindered rotor corrections. The data from the scans were fit using a Fourier series and separated from vibrations following Goldsmith et al. 26 In a few cases, relaxed scans broke hydrogen bonds but did not reform them during the scan and rigid scans were not accurate. The Fourier fits from these calculations gave unrealistic negative energies, so a cosine fit was used instead. Moments of inertia for each rotation were estimated using the lowest-energy conformer geometry. Vibrations were treated harmonically and their frequencies were scaled differently based on the usage. 27 Scale factors of 0.99 28 and 1.004 29 were applied to the frequencies for zero-point energy and the vibrational partition function, respectively.
Standard enthalpies, entropies, and heat capacity were obtained from the partition function, with energy adjustments based on atoms, 30 bonds, 30 and zero-point energies. These were fit over the range 10 to 3000 K using two NASA polynomials in the standard CHEMKIN II format. 31 High-pressure-limit kinetics used canonical transition state theory with 1-D Eckart tunneling corrections. Unlike the thermodynamic calculations, atom and bond corrections were not applied when determining kinetic parameters. Twenty data points from 180 to 1500 K were fit to the modified Arrhenius form (k = AT n exp(ÀE a /RT)).
Pressure-dependent networks were constructed for a-, b-, and g-isobutanol radicals. Pressure-dependent rate coefficients for each network were obtained from the master equation at 20 temperatures from 180 K to 1500 K and 20 pressures from 0.01 bar to 100 bar. The reservoir state method was used to compute the phenomenological rate coefficients, 32 based on the implementation described in Allen et al. 33 However, since Eckart tunneling was used in this work, the method was modified slightly; the cutoff between Boltzmann-distributed reservoir and the excited states was taken at the energy where the microcanonical reaction rate was 1% of the rate at the transition state's zero-point energy. This allows the transition point between excited states and Boltzmann-distributed states to adjust with the amount of tunneling. The maximum energy separation was 2 kJ mol À1 with at least 200 energy grains per well. Lennard-Jones parameters of s = 4.64 Å and e = 318 cm À1 were calculated using the 1-D minimization method with gRO 2 and N 2 as a bath gas. 34 Collisional energy transfer parameters of hDE down i = 250 Â (T/300 K) 0.85 cm À1 were used, a value that has been used for a similar system of 2-butanol radicals in helium. 35 For the three barrierless entrance channels of R + O 2 , inverse Laplace transform (ILT) of analogous reaction rates was used to obtain k(E) from k(T). 33 For aR + O 2 , we used the high-pressurelimit rate of the analogous ethanol reaction, CH 3 CHOH + O 2 , from Zádor et al. 20 Rates for oxygen addition to bR and gR were taken from the rate rules developed by Miyoshi for tertiary and primary additions, respectively. 36 Pressure-dependent rates were fit at each temperature and pressure to Chebyshev polynomials. Thermodynamic parameters, pressure-dependent and high-pressure-limit rate coefficients are available in CHEMKIN II format in the ESI. †

Monte Carlo uncertainty analysis
Previous work has highlighted important sources of uncertainty in similar systems. Uncertainty analysis for oxygen addition to the propyl radical in Goldsmith et al. showed collisional energy transfer rates, zero-point corrected energies, and location of barrierless reactions to be important. 37 Xing et al. also found the Lennard-Jones parameter and tunnelling frequency important in their study of ethanol decomposition. 38 Allen et al. showed that the method of solving for phenomenological rates can also impact pressure-dependent rates. 33 To bound the accuracy of the results in this work, multiple input parameters into the pressure-dependent solver were varied for each of the networks using a Monte Carlo approach. From these calculations, distributions of net rates of oxygen addition and branching ratios of the hydroxyalkly + O 2 reaction and the unimolecular RO 2 reaction can be computed.
In this work, seven categories of parameters were varied: zero point energies, rate of reactions input into ILT, negative frequencies of transition states, preexponential factor of collisional energy transfer, temperature exponent of collisional energy transfer. Lennard-Jones parameter, and the method of solving for phenomenological rates. The zero point energies of each stationary point were treated with a normal distribution with standard deviation of 10 kJ mol À1 , which gives a 95% confidence interval similar to the error reported by Zádor et al. for highly multi-reference energies. The rate of each reaction input into the ILT algorithm was varied with a lognormal distribution with a standard deviation of two, given the rate used was for a similar barrierless reaction. Negative frequencies of transition states, which impacts tunnelling, were treated as a lognormal distribution with standard deviation of 0.2. The preexponential factor of the collisional energy transfer was taken as a lognormal distribution with standard deviation of 0.5, and its exponent was a normal distribution with standard deviation of 0.15. The Lennard-Jones s parameter was assumed to be a lognormal distribution with a standard deviation of 0.2. The last four uncertainties described were used in Goldsmith et al. and Xing et al. 37,38 Two methods to solve phenomenological rates were used: modified strong collision and reservoir state. Reservoir state method was weighted with three times the probability of modified strong collision because it was seen as potentially more accurate. 33

Species and reaction naming
In total, this study involves 47 species and 38 transition states, so there is a need to develop a succinct and descriptive nomenclature. To ease discussion, a, b, and g always refer to the carbon location relative to the alcohol group, even if multiple functional groups are present.
All compound names are listed next to their corresponding structure on the three networks in Fig. 3, 6 and 9. Stable compound names are italicized throughout the text.
For contrast, reaction names are bolded. Each reaction starts with the Greek letter of the network to which it belongs. All reaction names are listed on the potential energy surfaces in Fig. 3, 6 and 9.

Results and discussion
For each peroxy radical, the pathways shown in Fig. 2 were evaluated. A QOOH species was found for each unique abstractable hydrogen. For each of these intermediates, the lowest-energy b-scission reaction, cyclic ether formation, and, if available, water formation pathways were found. In total 38 transition states were obtained.
In the following sections we describe major results from each peroxy radical individually, showing potential energy surfaces, microcanonical rates, pressure-dependent branching ratios, and comparisons to other quantum calculations. We then discuss the uncertainty of the pressure-dependent branching ratios and overall rates. We end by comparing the rates to a detailed combustion mechanism.

a-Peroxy network
The lowest-energy pathway in the a-peroxy network, shown in Fig. 3 involves a two-step HO 2 elimination of the peroxy radical. The first step, called aAdductFromRO2, forms a hydrogen bonded complex, called aadduct. This has an activation energy lower in energy than the separated products. The second step, aAdductSplit, involves breaking the hydrogen bond to form isobutanal and HO 2 . This barrierless reaction was included in the master equation by using an ILT of the collision-limited rate of the products. 33,39 The TS for aAdductFromRO2 is a concerted reaction with the SOMO perpendicular to the COO plane; the analogous reaction of other a-hydroxy peroxy radicals has been discussed previously. 20 The a-peroxy network also includes three intramolecular hydrogen abstraction reactions from the peroxy radical to form alkyl or alkoxy radicals. The lowest-energy direct hydrogen abstraction involves the formation of the alkoxy product, aQOOH [O]. Unlike the formation of aQOOHb and aQOOHg, this does not require the breaking of a hydrogen bond, allowing for a lower barrier despite a higher ring strain and larger bond dissociation energy. 40 For the other hydrogen abstraction reactions, the reaction with a six-membered transition state forming aQOOHg has a lower barrier than the one with a fivemembered transition state forming aQOOHb, due differences in ring strain.
This network also includes subsequent isomerization reactions, aAlkoxyIsomFromb and aAlkoxyIsomFromc, forming aQOOH [O]. These subsequent isomerization reactions have even higher barriers than both the direct isomerization from aRO 2 (aAlkoxyIsom) and the lowest-energy decomposition pathways (aEpoxyFrom and aC4EtherFrom). Due to the high barriers of these alkyl-based isomerizations, they were not calculated for the band g-networks.
The products formed from the isomerizations, aQOOHb, aQOOHg, and aQOOH[O], each have a number of decomposition pathways. The alkoxy pathway has three b-scission reactions that break a bond connected to the a-carbon. The aQOOHb behaves similar to peroxy radicals from non-oxygenated fuels, with the formation of an epoxide (in this case trisub_epoxy) + OH, or an alkene (ibutenol) + HO 2 as the main channels. The aQOOHg decomposition can undergo reactions aDoubleb-scission, aC4Et-herFromc, and a-bscissionFromc. The last of these reactions has the highest barrier and forms a methyl radical and alkene. Since this reaction is unlikely to dominate at any realistic condition, reactions forming both alkenes and methyl radicals were excluded for the other peroxy networks. The other b-scission reaction involves the simultaneous breaking of two bonds, similar to the Waddington mechanism described for b-peroxy radicals, 17 but with a higher barrier. aC4EtherFromc has the lowest decomposition barrier for aQOOHg, but it is still significantly higher than the decomposition reactions of aQOOHb and aQOOH [O].
From the quantum calculations, we derived density of states and microcanonical rates, with selected microcanonical rates shown in Fig. 4. The formation of aAdduct is the dominant pathway, as it is faster than other channels by at least an order of magnitude at all energy levels. For aAdduct, shown in Fig. 4b, most energies result in the breaking of the hydrogen bond to form HO 2 and isobutanal. Since this rate is based on the ILT for forming aAdduct from HO 2 and isobutanal through microcanonical equilibrium, the reaction rate reported here has only order of magnitude accuracy. However, given the wide gap in Fig. 4b, one or two orders of magnitude change in the aAdductSplit rate will not shift the major product. Fig. 5a shows the major products from phenomenological rate coefficients. Under practically all conditions studied, aR + O 2 proceeds via a well-skipping reaction to form HO 2 and isobutanal. Even the stabilized aRO 2 formed at low temperature and high pressure forms the same product or the hydrogenbonded adduct, as shown in Fig. 5b. This result is consistent with isobutanol oxidation experiments between 500-700 K, which show isobutanal as a major product. 41 The main products, HO 2 and isobutanal, could also form by direct hydrogen abstraction of aR by O 2 , which was not studied in this work. However, based on the work of an analogous ethanol network by Zádor et al., this channel is unlikely to be important even at temperatures as high as 1000 K. 20 Because the dominant channel from aR + O 2 makes HO 2 , and not a carbon-centered radical, aR does not lead to significant low-T chain-branching via second O 2 addition. (The HO 2 can lead to chain branching via H 2 O 2 , but typically that only becomes important above 750 K). Since the aR is predicted to be the dominant radical from isobutanol combustion, 11,13,15 a large negative temperature coefficient effect is not expected. This is in agreement with previous experiments that found no negative temperature coefficient region in isobutanol combustion.  , commonly refered to as the Waddington mechanism, 17 is often the major sink of bRO 2 radicals. 18,41 The microcanonical rates of the major intermediates are shown in Fig. 7. Of all the pathways available to bRO 2 , three pathways dominate at different energies. At the lowest energies (oÀ65 kJ mol À1 , i.e. below the zero-point energy of bQOOH[O]), b-aQOOHIsom, which occurs through tunneling, is the fastest reaction. After forming bQOOHa at this low energy, it will likely form trisub_epoxy + OH through bEpoxyFroma, as shown in Fig. 7c.
At moderate energies, the fastest reaction involving bRO 2 is bAlkoxyIsom. The rates of the alkoxy intermediate are given in Fig. 7b. The reverse reaction for bQOOH[O] re-forming bRO 2 (bAlkoxyIsom) is favored at most energies and the b-bscission-FromAlkoxy is favored only at high energies. The strong reversibility of the alkoxy isomerization allows for other secondary pathways to contribute to product formation.
The reverse reaction rate re-forming bRO 2 in Fig. 7b is nonmonotonic, decreasing at higher energies. This is because the zero-point-corrected energies of bQOOH[O] and bRO2Isom are similar (differing by only 2 kJ mol À1 ) and because bQOOH[O] has more freedom of motion than bRO2Isom. This causes the ratio of density of states between transition state and reactant to reach a peak around À60 kJ mol À1 . Because this isomerization is so fast, in the energy range from À60 kJ mol À1 through À10 kJ mol À1 , bRO 2 and bQOOH[O] equilibrate on a 10 ps timescale.
Based on Fig. 7a, most high-energy bRO 2 will dissociate back into bR and O 2 , which is essentially non-reactive. The fraction that does react mostly forms alkene, ibutenol, and HO 2 through b-aHO2ElimFromRO2 and b-cHO2ElimFromRO2. Fig. 8a shows the major products from phenomenological rate coefficients. In the case of bR + O 2 , the two HO 2 pathways become dominant at higher temperatures. The Waddington pathway (forming acetone, CH 2 O, and OH) will dominate at moderate temperatures and low pressures. Stabilized bRO 2 forms at the moderate and lower temperatures.
The products of a stabilized bRO 2 look quite different than for the entrance channel due to lower available energy. For example, epoxides only form in substantial amounts at low temperatures and when bRO 2 is stabilized. While epoxides have not been observed in isobutanol oxidation, 41 most studies have focused mostly on temperatures higher than 500 K.
At temperatures between 500 and 1000 K and pressures above 10 5 Pa, Fig. 8 shows that bR + O 2 primarily forms stabilized bRO 2 , and stabilized bRO 2 primarily forms bR + O 2 , indicating a partial equilibrium between bR + O 2 and bRO 2 . Under these conditions, the secondary products, not explicitly shown in Fig. 8, would be the dominant net pathways. At lower temperatures, the calculations in this work suggest that the dominant secondary products are trisub_epoxy and OH at low temperatures and acetone, CH 2 O, and OH at higher temperatures.
With the exception of temperatures below 250 K and above 10 5 Pa, very few stabilized QOOH radicals are predicted to be formed. As with the a-network, the lack of QOOH radicals suggests slow preignition chemistry.
3.2.1 Comparison with previous calculations. There have been more previous quantum calculations of the bRO 2 surface than of the other two isomers' surfaces because bR can also form from OH addition to isobutene. 18,19,42 This section compares this work's results with those of two other quantum studies. They include all the reactions in this work except for bH2OForm, 16 which has not been published for the b-isobutanol surface before.
The barrier heights in this work generally agree with those previously published, though with a few slight differences (see Table 1). 18,19 For bEpoxyFromc, our reaction barrier is 5 kJ mol À1 and 15 kJ mol À1 lower than the barriers computed by Lizardo-Huerta et al. and Sun et al., respectively. This may be due to the presence of intramolecular hydrogen bonding in the TS conformers in this work and that of Lizardo-Huerta et al. but not in Sun et al. 18,19 For b-cHO2elimFromRO2, the barrier in this work is 10 kJ mol À1 lower than the rate computed by Sun et al., but only 3 kJ mol À1 lower than the barrier computed by Lizardo-Huerta et al. The barrier height for bHO2elimFroma determined in this study is between those calculated in the other two studies. The calculated barriers for the rest of the reactions fell within 8 kJ mol À1 of those in the other two studies. A recent study by Li et al. looked at bAlkoxyIsom and b-bscissionFromAlkoxy using two methods, CCSD(T) and DLPNO-CCSD(T), with basis set extrapolation. 42 Relative to the current work, DLPNO-CCSD(T) resulted in a larger difference in barrier heights than CCSD(T) did (see Table 1). For these two reactions larger barrier hight differences exist between the present work and that of Li et al. than between this work and that of Sun et al. 18 Table 1 reduced this difference in rates to an average factor of 15, indicating that energy differences do not entirely account for the difference in rates; this indicates  Fig. 6, with bR + O2 being the reverse of bRO2Form, and bRO2Isom being the reverse of bAlkoxyIsom (for subplot b) and b-aQOOHIsom (for subplot c). The energy values correspond to those in Fig. 6 and are relative to that of bR + O 2 .  that differences in calculated partition functions are also contributing. In the work of Li et al., only the bAlkoxyIsom using the DLPNO-CCSD(T) method, differed by more than a factor of 10 over the same temperature range, due in part to the 16 kJ mol À1 barrier height difference. The major pressure-dependent pathways in this work, shown in Fig. 8, show general agreement with the conclusions drawn in Sun et al., though our calculations predict more alkene formation. Both studies predict that bR + O 2 gives stabilized bRO 2 at higher pressures and lower temperatures, with the switching of the major pathway at around 10 5 Pa and 1000 K. With lower pressures (T = 1000 K), our calculations suggest alkene formation will dominate over the Waddington reaction, whereas those of Sun et al. suggest the reverse. Other major products in this work include ibutenol and trisub_epoxy, which Sun et al. also found as large secondary pathways. 18 For the reactions of stabilized bRO 2 , this work and that of Sun et al. find the reaction forming bR + O 2 to dominate for all pressures at 800 K. For lower pressures (at 800 K), this work and that of Sun et al. predict the main secondary channel to be the Waddington pathway. 18 However, at high pressures, Sun et al. found bQOOH [O] formation to be the main secondary channel, but this work finds it to be one of the slower pathways. This difference may have arisen from the approximations made to get phenomenological rates. If we use the modified strong collision method to obtain pressure-dependent rates from the master equation instead of the reservoir state method, we obtain more rapid formation of bQOOH [O], in better agreement with the results of Sun et al. (see Fig. S3 in the ESI, † for the major products obtained using the modified strong collision method).
The reactions shown in Fig. 8 do not explicitly include the non-reaction path, in which excited bRO 2 reforms bR and O 2 (though this is taken into account implicitly by an overall decrease in reaction rate). This non-reaction can be important when determining how much of bR will b-scission to form isobutene instead of adding O 2 , which the model by Merchant et al. had trouble accurately predicting. 11 Fig. S45 in the ESI, † shows the fraction of excited bRO 2 which decomposes back to the bR. For the b-network at high temperatures, over 95% of excited bRO 2 re-forms bR + O 2 , in agreement with the results in Sun et al. 18 Given the small region showing the formation of acetone, CH 2 O, and OH in Fig. 8a, the Waddington pathway appears less important than has been suggested previously. 41,43 This is due to sensitivity in the kinetic rates and the masking of secondary products in Fig. 8. The fraction going to this path is highly dependent on the barrier heights. A sensitivity study in which the barrier height of b-cHO2FromRO2 was raised by 10 kJ mol À1 , to correspond to that of Sun et al., increased the range of conditions where the Waddington pathway was most dominant (Fig. S1 of the ESI †). In addition, in the region around 800 K and 50 bar, where bR + O 2 and stabilized bRO 2 are essentially reversible, the major secondary pathway is the Waddington path.

c-Peroxy network
As shown in Fig. 9, gRO 2 has four isomerization reactions and a direct path to form HO 2 . The lowest-energy isomerization involves a 6-member transition state to form gQOOHa, which can decompose via a number of pathways: cDoublebscission-Froma, cAldolFroma, cC4EtherFroma, and cH2OForm.
The other isomerization reactions are less important. Despite only slightly higher product energies, the decomposition pathways of gQOOH and gQOOH[O] all have high barriers, making them unlikely to be important at most conditions. The isomerization reaction with the highest barrier, c-bQOOHIsom, forms cQOOHb. This intermediate has two lower-energy decomposition paths forming either an alkene or a cyclic ether. The direct HO 2 elimination pathway from gRO 2 has a similar barrier to c-bQOOHIsom. Fig. 10 shows the microcanonical rates for relevant intermediates. For gRO 2 , the c-aQOOHIsom pathway is dominant at the lower energies and is taken over by the rate to reform gR + O 2 only at higher energies. At high energies, other product channels such as c-cQOOHIsom and cHO2ElimFromRO2 also have competitive rates.
At the lowest energies shown in Fig. 10b, gQOOHa will reform gRO 2 , though this is taken over at slightly higher energies by cH2OForm. At the highest energies, many reactions can effectively compete, as is the case for the microcanonical rates of gRO 2 .
The large number of potential pathways at high energies is shown by the dark shading at high temperatures in Fig. 11a. This includes a mixture of HO 2 + galkene, CH 2 O + OH + propene1ol, and OH + disub_c4ether. Between 800-1100 K, the dominant product is H 2 O + galdoxy, which agrees with the results of Welz et al. 16 Below 600 K, the major product is stabilized gRO 2 .
Stabilized gRO 2 starts with much lower internal energy, which prohibits the higher-energy reactions that are available to the entrance channel gR + O 2 . At the lowest temperatures simulated, shown in Fig. 11b, stabilized gQOOHa radical is formed. Above 400 K, the well-skipping reaction forming water and the alkoxy radical is the major pathway. At the highest temperatures and pressures, dissociation to reform gR + O 2 is the dominant pathway.
Relative to the other two networks, the g-network is most likely to produce a stabilized QOOH radical. This species can undergo unimolecular reaction or reaction with O 2 (via addition of H-atom abstraction). Comparison of these reaction rates (high-pressure unimolecular reactions from this work, bimolecular reactions with O 2 by Sarathy et al. 15 ) indicates that gQOOHa reacts almost entirely through bimolecular H-abstraction channel in air, forming HO 2 and a ketohydroperoxide, over the range of conditions in this study. The ketohydroperoxide will likely decompose into two OH radicals, two carbonyls, and a carbon monoxide, which is the dominant decomposition pathway for other g-separated, aldehyde-containing ketohydroperoxides. 12 Overall, this pathway for isobutanol oxidation can produce one HO 2 and two OH radicals from a single initiation. If important enough, this radical branching could create a negative temperature coefficient region. However, since our calculations show gQOOHa only forming below 500 K, this channel is unlikely to impact combustion significantly.

Comparison with previous calculations. Welz et al.
calculated an extensive potential energy surface for the g-network, but did not calculate the corresponding kinetic rates. The present work includes all but one of the reactions examined by Welz et al. for the g-network and added also includes decomposition pathways for three isomers: gQOOH[O], gQOOHg, and gQOOHb. 16 The reaction not included in this present work, which involved a direct path from gRO 2 to galdol + OH, was instead replaced by a reaction from gQOOH[O] to form the same products: galdol and OH, based on the IRC trajectory (view IRC trajectories in the ESI †). This reaction, named cAldoxyHabs, involves a hydrogen abstraction of the a-hydrogen with a simultaneous breaking of the O-O bond. Table 2 shows the barrier height differences between this work and that of Welz et al. Decent agreement, within 8 kJ mol À1 , occurs for all but three reactions. We report barrier heights for cAlkoxyIsom to be 20 kJ mol À1 lower and that of cAldolFroma to be 20 kJ mol À1 higher than those in Welz et al. The barrier for cC4EtherFroma was calculated to be 12 kJ mol À1 higher than that in Welz et al. A smaller difference of 6 kJ mol À1 occurs for cH2OForm, the dominant pathway at most conditions. These differences may be due to different level of theory used in the two studies, possibly amplified by multireference effects. Using the barriers given by Welz for the cH2OForm reactions does not lead to significant shifts in predicted product distribution (see Fig. S1 in the ESI †).

Uncertainty analysis
We evaluated how the inputs into the pressure-dependent kinetic solver impact the product branching ratios using two methods. The first involved changing an input parameter a specified amount and viewing how it would change the Fig. 9 The potential energy surface for the g-peroxy radical. Energy includes zero-point energy and atom corrections. Species name and structures appear in black near the structure energy level. Reaction names appear in gray above the transition state energy.
branching ratios across the range of temperatures and pressures shown in Fig. 5, 8 and 11. This was done for three input parameters: barrier heights, collisional energy transfer, and the method of solving phenomenological rates. The results appear in Fig. S1-S3 of the ESI. † From these sensitivity studies, the largest change appeared when changing the method used in solving phenomenological rate constants, and the smallest change ocurred when varying the collisional energy transfer. This methodology was particularly useful when evaluating how a difference in barrier heights between this work and previous works impacted the conclusions in this work. This method does not give an overall uncertainty of the branching ratio given the multiple sources of uncertainty.
To evaluate this overall uncertainty, Monte Carlo calculations were conducted varying seven types of parameters: electronic energies, tunneling frequency, Lennard-Jones parameters, collisional energy transfer parameters, rates of reaction used for barrierless reactions, and the method of solving for phenomenological rates. Of the 2000 runs for each network, 1008, 932, and 1194 were successful for the a-, b-, and g-networks, respectively. Unsuccessful calculations typically resulted from initial conditions where a transition state is lower in energy than the corresponding reactants or products, which occurred more frequently in this work than it would in other work 37 due to the larger network size and a higher uncertainty value being placed on single-point energies.
Branching ratios from the resulting models were calculated for bimolecular R + O 2 and unimolecular RO 2 reactions at three  Fig. 6, with cR + O2 being the reverse of cRO2Form, and cRO2Isom being the reverse of c-aQOOHIsom. The energy values correspond to those in Fig. 9 and are relative to that of gR + O 2 .  conditions: 300 K and 1 Â 10 5 Pa, 600 K and 3 Â 10 5 Pa, and 900 K and 1 Â 10 6 Pa. The median value and 90% confidence limits are reported in Tables S1-S18 of the ESI. † Also included are the two inputs parameters that had the highest absolute value of Spearman rank coefficient, indicating they had a large impact on the branching ratio. In general, branching ratio uncertainty was higher in regions where the original model showed competing pathways. For example, the reaction of bR + O 2 at 300 K ant 10 5 Pa in Fig. 8a shows that bRO 2 is the only major product. The 90% confidence interval for the branching ratio of bR + O 2 forming bRO 2 , as determined from Monte Carlo, is 0.96-1.00, indicating there is also little uncertainty in this value. On the other hand, the reaction of gRO 2 at 600 K and 3 Â 10 5 Pa has three major pathways, shown in Fig. 11b. For each of these three pathways, the 90% confidence intervals for the branching ratio spans at least 0.01 to 0.90, indicating that any one of these reactions could reasonably be the dominant product at these conditions. Given this uncertainty, any user of this model should ensure their conclusions are not sensitive to the branching ratios discussed in this work when the original model predicts a branching ratio of the dominant channel less than 90% (demarcated with darker shading in Fig. 5, 8 and 11).
The factors that most contributed to branching ratio uncertainty were single-point energies, the rate used in ILT, the method used to solve phenomenological rates, and negative frequencies of transition states. Single-point energies and the frequency of transition states were more important at lower temperature, whereas the rates used to solve for the ILT were important for the branching ratio of RO 2 at higher temperatures. The method used to solve for phenomenological rates was important where well skipping occurs, including for reactions of RO 2 at high temperatures and reactions of R + O 2 at all three temperatures. Relative to these sources of uncertainty, collisional energy transfer and Lennard-Jones parameters were not large sources of uncertainty when determining branching ratios.
Tables S19-S21 (ESI †) show the uncertainty in overall reaction rate for various pathways and the major factors contributing to them. For the oxygen addition reaction, the only major source of uncertainty comes from the rate used in the ILT. For stabilized RO 2 the uncertainty came from three sources: the electronic energy of the peroxy radical, the electronic energy of the isomerization transition state, and the rate used in the ILT, with the last being more important at higher temperatures.
While the uncertainty analysis here was conducted at only three conditions, the ESI, † contains the models and code to estimate uncertainty anywhere within 180-1500 K and 10 3 -10 7 Pa, enabling such an uncertainty assessment to be conducted for any given set of reaction conditions.

Comparison with a detailed combustion mechanism
The reaction rates from these calculations can be compared to estimates from the detailed isobutanol mechanism by Sarathy et al. 15 Two other isobutanol mechanisms, published by Hui et al. 7 and Merchant et al., 11 do not include explicit representation of peroxy radicals, so these cannot be compared directly to the calculations here.
The Sarathy mechanism uses pressure-independent rate coefficients, which were compared to our high-pressure-limit rate coefficients. Of our 38 reactions, only 25 correspond directly to reactions in the Sarathy mechanism. The Waddington mechanism, which the Sarathy mechanism represents as one step, was also added and is compared to the bAlkoxyIsom reaction in this work. These reactions were compared at 50 temperatures between 500 and 1000 K. At each temperature, a ratio was taken of the rate in this work to rates in the Sarathy mechanism, and the mean of the ratios at the 50 temperatures is designated as the deviation between the two models. Fig. 12 shows the deviation for the 25 reactions as a histogram, with color coding indicating different types of reactions.
Given that the rates in Sarathy et al. consist of estimates, obtaining a rate within two orders of magnitude is reasonable. Good agreement is found for HO 2 elimination reactions from RO 2 and most of the RO 2 isomerization reactions.
Other reaction types deviate more strongly, with Sarathy et al. estimating epoxy formation around a million times slower than this work suggests. In the Sarathy mechanism, these rates originate from rate rules based on alkane fuels; such reactions have an activation energy of around 92 kJ mol À1 , 44 which is substantially higher than activation energies given from quantum calculations for ethane oxidation, between 50-70 kJ mol À1 . 45 In addition, the rates in Sarathy et al. are 2-3 orders of magnitude slower than those recommended by a more recent structureactivity relationship. 46 The barrier heights for epoxy formation reactions from isobutanol in this present work and other studies 18,19 are also approximately 20 kJ mol À1 lower than those reported for alkane and ethanol oxidation. 20,45,46 The large differences between the estimates in Sarathy et al. and those in this work can be attributed to the use of older structure-activity relationships and/or differences between alkane chemistry (on which the Sarathy et al. rate rule is based) and the chemistry of isobutanol oxidation.
The two classes of QOOH decomposition shown in Fig. 12, 'QOOH b-scission' and 'HO 2 from QOOH', are also predicted by Sarathy et al. to be 10-10 4 times slower than determined in this Fig. 12 Histogram of the average ratio of rates calculated in this work and those estimated in the Sarathy mechanism 15 between 500-1000 K, grouped by reaction type.
work. The slower decomposition of QOOH would increase the predicted concentrations of QOOH; this would lead to an increased predicted importance of O 2 addition to QOOH, which can leading to a faster simulated ignition than what would be predicted based on this present work.

Conclusion
The calculations in this work update the chemistry of peroxy radicals formed in isobutanol oxidation using pressuredependent rates for all three peroxy radical isomers. We obtain general agreement with the potential energy surface proposed The three isomers in this work react through different channels. The fate of aR + O 2 almost entirely results in HO 2 + isobutanal, which agrees with the results of Zádor et al. 20 The reaction of gR + O 2 almost entirely forms H 2 O and alkoxy at lower temperatures, as suggested by Welz et al., with a mixture of products at higher temperatures. The b-network predominately proceeds via the Waddington pathway at lower temperatures, though there is much more competition with other pathways at higher temperatures than other studies suggest. The rates from this work can explain the lack of two-stage ignition from isobutanol combustion.
Given the potential use of isobutanol as a biofuel, these rates can be integrated into detailed kinetic mechanisms to help improve the accuracy of combustion models for low-temperature combustion (500-700 K). Refining the oxidation chemistry of key products predicted in this work, such as alkenols and larger aldehydes, might help improve accuracy of detailed models still further. Due to the potential for multireference characteristics, important reactions in this work (such as bAlkoxyIsom, b-bscissionFromAlkoxy, and b-cHO2elimFromRO2) could be better constrained using multireference calculations. Further, an improved understanding of low-temperature isobutanol oxidation may require the inclusion of the O 2 reactions with stabilized QOOH radicals, such as gQOOHa.

Conflicts of interest
This fundamental research was funded by grants to MIT from the NSF and from BP. We note that BP is a part owner of BUTAMAX, a company commercializing the biofuel isobutanol. The authors declare that they have no financial interest in isobutanol commercialization and no conflict of interest.