Pressure- and temperature-dependent combustion reactions

Combustion and Flame 158 (2011) 602–617

Contents lists available at ScienceDirect

Combustion and Flame j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m b u s t fl a m e

Pressure- and temperature-dependent combustion reactions David M. Golden a,⇑, John R. Barker b,⇑⇑ a b

Department of Mechanical Engineering, Stanford University, Stanford, CA 94305-3032, United States Department of Atmospheric, Oceanic and Space Sciences, University of Michigan, Ann Arbor, MI 48109-2143, United States

a r t i c l e

i n f o

Article history: Received 5 May 2010 Received in revised form 17 July 2010 Accepted 24 August 2010 Available online 28 September 2010 Keywords: Master equation Pressure-dependent RRKM Energy transfer MultiWell

a b s t r a c t In combustion systems, many reactions are simple thermal unimolecular isomerizations or dissociations, or the reverse thereof. It is well understood that these reactions typically depend on temperature, pressure and the nature of the bath gas. These kinds of reactions are a subset of the more general behavior that can be described as free radical association reactions that produce highly energized intermediates, which can isomerize or dissociate via multiple chemical pathways. Each reaction rate depends on excitation energy and all of the competing reactions occur in competition with collisional activation and deactivation. These complicated multi-well, multi-channel reaction systems can only be simulated accurately by using master equation techniques. In this paper, master equation calculations are discussed for several examples of reactions important in combustion (and atmospheric chemistry). Current master equation codes are based on statistical RRKM reaction rate constants (including quantum mechanical tunneling) and simplified models for collisional energy transfer. A pragmatic semi-empirical approach is adopted in order to compensate for limited knowledge. The reaction energies needed for RRKM calculations are usually obtained from quantum chemistry calculations, which are often of limited accuracy and may be adjusted empirically. Energy transfer cannot be predicted accurately and must be parameterized by fitting experimental data. For combustion modeling, the master equation results are usually expressed as chemical reactions with rate constants fitted to empirical algebraic equations. However, the results may be expressed more accurately by interpolating from look-up tables. Several current research issues are also mentioned, including the effects of angular momentum conservation, vibrational anharmonicity, slow intramolecular vibrational energy redistribution, and assumptions surrounding the details of collisional energy transfer. Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction Understanding the chemistry of reaction systems (e.g. combustion and atmospheric chemistry), reactants, products and the rates of interchange among them is often addressed by creating a ‘‘model” of a given reacting system with which to compare to experiments. The model typically consists of chemical species, their thermochemistry and ordinary differential equations (ODE’s) that mathematically describe the chemical processes. The parameters of these ODE’s are rate coefficients that may be functions of temperature, pressure and the nature of the medium. The chemistry is best described as the ‘‘mechanism” for the process. ‘‘Modeling” includes the act of assigning specific parameters to individual processes in the mechanism. In practice these appellations are used interchangeably. The present paper addresses issues related to accurate descriptions of the rate constants for individual ‘‘pressure-dependent” ⇑ Corresponding author. Fax: +1 650 723 1748. ⇑⇑ Corresponding author. E-mail addresses: [email protected] (D.M. Golden), jrbarker@umich. edu (J.R. Barker).

reactions and how they vary as functions of pressure and temperature. The discussion here focuses on free radical association reactions that produce highly energized intermediates, which can isomerize among multiple isomers (‘‘wells”) or dissociate via multiple chemical pathways. Each reaction pathway (isomerization and dissociation) is considered to be a reaction ‘‘channel”. Simple dissociations and recombinations are a subset of this class, which is often labeled with the understandable yet classic misnomer ‘‘unimolecular” reactions. Such single-channel, single-well reaction systems were first described using RRKM theory, which is described below. Later, the more accurate master equation methods, which are also described below, were applied to simple reactions and then extended to multiple-channel, multiple-well systems. Master equation (ME) simulations have proved to be of great use in interpreting, interpolating, and extrapolating empirical kinetics data on complex reactions initiated by photoactivation, chemical activation, or recombination. There have been quite a few papers in the combustion literature describing RRKM and ME calculations. Most of these are erudite descriptions of the physics of these processes with applications to some specific reaction systems. Often ab initio calculations of relevant potential

0010-2180/$ - see front matter Ó 2010 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustflame.2010.08.011

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

energy surfaces (PES’s) are an integral part of these discussions. The present paper places less emphasis on the technical computation and more emphasis on pragmatic methods for accurately describing the chemical reactions. In some cases, PES’s have not even been computed for the reactions, but rational heuristic methods are presented instead for estimating the necessary parameters.

1.1. Models and mechanisms Often, for combustion, the experiments that the model is attempting to describe involve the pyrolysis and/or oxidation of a single fuel. Of course, of more general interest would be a real fuel mixture, which may consist of hundreds of chemical compounds. (The analog in atmospheric chemistry research might be the introduction of a single pollutant into a smog chamber.) The existing computer codes for solving these coupled ODE’s are capable of identifying the reactions to which the measured property is most sensitive, for a given set of input parameters. The input parameters, specific reactions and the corresponding thermochemistry and rate coefficients are found in various ways. Sometimes mechanisms and models can be generated from a combination of chemical intuition and specific knowledge gleaned from earlier studies. This is not always possible for more complex systems and the development of mechanisms and models becomes a difficult task. In most important applied problems, there are many reactions and intermediates for which, the rate constants and thermochemistry are not known and are often difficult to measure. Advances in computational chemistry, numerical methods, and computer hardware have made it possible to construct and solve kinetic models capable of reliably modeling such complicated systems. It is now possible to generate mechanisms using codes that have been developed for this purpose. This software must first be able to generate all possible combinations of chemical species and then be able to identify and eliminate those that are not chemically possible. In addition, the software must have built-in methods for estimating thermochemistry and rate constants. An example of mechanism generation software [1] can be found at: http://rmg.sourceforge.net/. See also Carstensen and Dean [2] for Rate Constant Rules for the Automated Generation of Gas-Phase Reaction Mechanisms. Once the most sensitive reactions (and/or thermochemical values) for a specific target experiment have been identified, it is possible to modify the corresponding rate constants and thermochemical values to fit the target experiment. This is sometimes referred to as ‘‘tuning” the model. Of course, tuning will improve the calculation for the particular target experiment, but the modifications may result in a poorer fit to some other target experiment that is sensitive to the same parameter. A more sophisticated ‘‘optimization” procedure can be employed. An example and discussion of the procedure is found in relation to the mechanism developed for natural gas combustion known as ‘‘GRI-Mech” [3] (http://www.me.berkeley.edu/gri_mech/). As explained on the above web site, the systematic optimization procedure used to produce the GRI-Mech mechanisms involves the following steps: 1. Assemble a reaction mechanism consisting of a complete set of elementary chemical reactions. 2. Assign values to their rate constants from the literature, by quantum chemistry calculation, or by judicious estimation. Treat temperature and pressure dependences in a proper and consistent manner. Also evaluate error limits, and the thermodynamics used for the equilibrium reverse rate constants. This creates a mechanism based on prior knowledge.

603

3. Search the literature for reliable experiments (‘‘targets”) that are relevant to the combustion that the model will be designed to describe. These experiments should include key combustion properties that the mechanism is to predict. 4. Use a computer code to solve the reaction mechanism kinetics, computing values for the observables of these ‘‘target” experiments. Also apply sensitivity analysis to determine how the model input rate constants affect the result. Compare computed results with data. Typically, these will not match many of the targets. 5. The choice of experimental targets should be sensitive to a representative sample of the rate parameters, under a representative set of conditions. Many parameters will significantly affect more than one target. Select, according to sensitivities and uncertainties, those parameters making the largest impacts on a given target. These are the potential optimization parameters. 6. Map the model response by repeating computations of the target observables for a minimum subset of combinations of possible values for the optimization parameters – within their appropriate error limits – according to a central composite factorial design. Try to include all ‘‘active” (i.e., significant) parameters. 7. Create, using the results of these factorial-design-directed calculations, the polynomial functions (the response surfaces) that mimic the results of the computer simulation for each target. This solution mapping technique in essence creates a representation of the predicted target values for the set of possible mechanisms within stated error estimates. The rate constants are normalized on a logarithmic scale, and a second order polynomial expression is typically created by regression analysis of the appropriate factorial design calculation set. Details of the factorial design and response surface fitting are given in Frenklach et al. [4]. Use the response surfaces to calculate target values that are then compared to measured values in an error function known as the objective function, which is then minimized. The values of the variables thus obtained are statistically the best values resulting from the universe of data considered as targets. The present paper is mostly concerned with Item #2 in this list. The above procedure will have generated a test set or model that can be used to predict other related experimental findings. This step is called ‘‘validation”. When new information is obtained, the optimization procedure should be repeated. At this point the model can be used to predict yet to be measured events. However, the uncertainty in the prediction is governed by the uncorrelated uncertainty and sensitivity of the above input parameters. Frenklach, Packard and coworkers have developed a methodology called ‘‘Data Collaboration” [5,6]. They point out that the space of input parameters and their uncorrelated uncertainties can be thought of as a hypercube. Using methods of robust control theory, they demonstrate that considering the universe of available data as a whole, the space of input data and their uncertainties is reduced from a hypercube to a much smaller volume defined by correlated uncertainties which they term ‘‘the feasible set”. These are values of the range of input parameters that mutually describe the universe of data. With this methodology they can describe model uncertainty and thus make predictions that contain uncertainty bounds constrained by the complete set of experimental data. 1.2. About master equation calculations ME simulations are often very successful in modeling experimental data [7–16]. Sophisticated Transition State Theory (TST) formulations for computing rate constants are also very successful

604

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

even when experimental data are not available [17–20]. Thermal dissociation, recombination, isomerization, chemical activation, and photoactivation, which are closely related, are described using the same ME theoretical framework [13]. Statistical rate theory (i.e. RRKM Theory) has been highly developed, making it possible to compute k(E, J) from first principles with good accuracy. However, there are times when ME simulations fail to describe the experimental data, or the empirical parameters appear to be un-physical even when a good description is achieved. ME calculations require microcanonical rate constants, k(E), and energy transfer parameters. To calculate k(E), workers in this field usually compute the ab initio properties of the critical points (i.e. local minima and saddle points) on the PES on a potential energy surface, and then use the corresponding vibrational frequencies, anharmonicities, and moments of inertia to compute microcanonical rate constants from Rice–Ramsperger–Kassel–Marcus (RRKM) theory [10–12,20]. Activation energies and energy transfer parameters are rarely known with sufficient accuracy [21] and instead are often adjusted to fit experimental data. The resulting semiempirical ME model can then be used to interpolate and extrapolate the experimental data. Together, total internal energy (E) and total angular momentum (with quantum number J) control the microcanonical rate constant, k(E, J), for recombination reactions and their reverse, unimolecular dissociation reactions [9–12]. These ‘‘barrier-less” reactions are particularly sensitive to angular momentum because they do not have intrinsic energy barriers (i.e. local maxima separating reactants and products on the potential energy surface) and therefore centrifugal barriers play a dominant role. In the non-equilibrium chemical reaction systems important in planetary atmospheres (including Earth), in combustion, in photochemical reactions, and other non-equilibrium systems, collisional energy transfer affects the energy and angular momentum (rotational) distributions and therefore plays a key role in determining the rates of chemical reactions. Collisional energy transfer competes with the chemical reactions and the net reaction rates and product distributions depend on temperature, on pressure, the nature of the collision partners, on how the reactions are initiated, and on other properties of the reaction system. To accurately model such systems requires formulation of a master equation (ME), which can account for all of the factors just mentioned. The current state of ME development and many applications have been reviewed recently by Pilling and Robertson [15], by Miller and Klippenstein [16], and by Fernandez-Ramos et al. [18]. In the present work, all computations of pressure and temperature dependences are performed using the Multiwell suite of codes [22,23]. Barker and coworkers have specialized in developing and applying stochastic methods [24,25] for carrying out such ME simulations [13,26–29]. The Multiwell website (http://aoss-research.engin.umich.edu/multiwell/index.php) describes these programs and lists papers citing their use.

2. Theory In this section, the properties of unimolecular reactions are described and then the basic RRKM theory is presented. Although the original RRKM theory has been largely superceded by master equation simulations, the theory is of considerable pedagogical interest and it serves to introduce a number of fundamental concepts. The master equation is then described with particular emphasis placed on the problems associated with its reduction from two dimensions (dependence on E, J) to one dimension (dependence on E). These problems will persist until practical methods are available for solving the two-dimensional ME. Until then, the two-dimensional ME will continue to be reduced (approximately) to the

one-dimensional ME, which can be solved using established methods. 2.1. RRKM Theory for single-channel, single-well reactions The thermal decomposition of simple molecules has been long noted to depend on both pressure and temperature. In fact, the nature of the bath gas, which is generally responsible for most of the pressure in the chemical system, is also a variable. At sufficiently high pressures the rate is first order and the rate constant is independent of pressure, while at sufficiently low pressures the rate is proportional to pressure. Studies of such systems have clarified the fact that these reactions proceed as a result of the competition between the decomposition of the thermally activated species and collisions with the bath gas that deactivates them. This competition also occurs when the molecule is photo-activated or activated by its production by a chemical reaction (i.e. chemical activation). Likewise many reactions proceed via a coupled set of transformations along a path on a potential energy surface that contains multiple reactant wells separated by barriers between them. In all of these cases the transformations among species are controlled by the competition between the reaction of an activated molecule and its collisional deactivation. The competition between collisional activation/deactivation and the decomposition step is embodied in the Lindemann mechanism [30], which can be written as follows:

A þ M  A þ M

ða; -aÞ

A  B þ C

ðb; -bÞ

where A is the unexcited reactant, M is an energy transfer partner (i.e. bath gas), the asterisk denotes vibrational excitation, and the products are B and C. Reactions (a) and (–a) are for energy transfer, which maintains the Boltzmann energy distribution of the reactant. Reaction (b) is the actual decomposition step, and Reaction (–b) is the recombination reaction. By invoking the pseudo-steady-state approximation for [A], the concentration of A, one obtains expressions for the rate of the global unimolecular decomposition reaction (kuni) and its reverse, the recombination reaction (krec). The global rate constants kuni and krec are related through the equilibrium constant Keq = kuni/krec and both depend on [M], the concentration of the energy transfer partner M, which represents all chemical species, since all are energy transfer agents. The rate constant for kuni takes the following form:

kuni ¼

k0 k1 ½M k0 ½M þ k1

ð1Þ

where k0 and k1 are parameters that depend on temperature and k0 also depends on the nature of the bath gas. At very low values of [M], the rate constant reduces to the ‘‘low-pressure limit” (klow  k0[M]), which is proportional to pressure. At very high pressures kuni approaches the ‘‘high pressure limit” (k1), which is independent of pressure. In the intermediate pressure range, the rate constant is in the ‘‘pressure fall-off” regime, where kuni falls below the magnitudes of both limiting rate constants. These general properties of kuni also apply to krec, because krec is proportional to kuni. For a more realistic description of such a reaction system, all of the rovibrational states must be considered and the dependence of the decomposition rate constant on energy (E) and angular momentum (quantum number J) must be taken into account: k(E, J). Note that E is the total internal energy (vibrations and rotations, but omitting translation of the center of mass; E = e + EJ in Fig. 1) measured from the zero point energy of the reactant molecule. RRKM theory includes the effects of the huge numbers

605

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

Fig. 1. Energy scheme for a ‘‘barrier-less” unimolecular reaction, including the effects of angular momentum conservation.

of rovibrational states and uses a microcanonical version of transition state theory (TST) to describe k(E, J). The RRKM expression for kuni is based on the Lindemann mechanism and the energy scheme shown in Fig. 1, which represents the ‘‘barrier-less” dissociation of a generic reactant molecule while conserving total angular momentum. The RRKM rate constant takes the following form (for the derivation, see standard monographs [10–12]):

kuni ¼

Q þ1 expðE0 =kTÞ hQ 1 Q 2

þ

þ

þ

Z

1

eþ ¼0

Gðeþ Þ expðeþ =kTÞdeþ 1 þ kðeþ Þ=kc bc ½M

kðe Þ ¼ Gðe Þ=hqðe þ Eo þ DEJ Þ

 DEJ ¼ kTð1  Iþ =IÞ +

k1 ¼

Q þ1 Q þ2 kT expðE0 =kTÞ Q 1Q 2 h

ð4Þ

ð2aÞ

which is identical to canonical transition state theory. At the lowpressure limit, kuni becomes

ð2bÞ

k0 ½M ¼

where E0 is the critical (threshold) energy for reaction (i.e. E0J for J = 0 in Fig. 1); Q1 and Q þ 1 are the partition functions for the adiabatic rotations in the excited molecule and in the transition state, respectively; Q2 is the partition function for the active modes in the reactant molecule (i.e. vibrations, internal rotations and nonadiabatic external rotations, which can rapidly exchange energy among themselves), e and e+ are the energies in the active modes of the molecule and the transition state, respectively (see Fig. 1). The ‘‘active” modes freely interchange energy, while the adiabatic modes are constrained in that regard. In particular, the rotational energy associated with conservation of angular momentum cannot be exchanged freely. The effect of angular momentum conservation is most important when there is no intrinsic energy barrier along the reaction path. For a given angular momentum quantum number J, the adiabatic rotational energy in the molecule is EJ and the adiabatic rotational energy at the transition state is EJ+. The difference between the rotational energies is due to the differences in structure, which result in different moments of inertia and rotational constants. The rotational energy difference is defined as DEJ ¼ EþJ  EJ . (Note that for bond fission reactions the moment of inertia at the transition state is larger than that of the molecule, so that DEJ is negative. Energy conservation requires this energy to be released during bond fission, which aids the dissociation process.) A reasonably accurate simplification is achieved by replacing DEJ in Eq. (2b) with its average value [12]:



The microcanonical rate constant k(e) in Eq. (2) is the rate, or frequency of passing through the transition state at energy e, G(e+) is the sum of states in the transition state, and q(e) is the density of active vibrational and internal rotational states in the molecule at e. Both G(e+) and q(e) are computed by including all of the vibrations, internal rotations, and the non-adiabatic external ‘‘Krotor”, which is assumed to be an active degree of freedom. Bimolecular rate constant kc is for bimolecular collisions, [M] is the total molecular concentration, and bc is a dimensionless collision efficiency, which depends on the identity of the bath gas. The strong-collision approximation corresponds to bc = 1. If e exceeds the energy threshold for reaction, reaction occurs, or energy is removed by collisions to levels lower than e. When bc < 1, this becomes the pseudo- or modified-strong-collision approximation. As discussed above, krec is related to kuni via the equilibrium constant at all pressures [31]. The essence of the competition between the reaction of an activated molecule and its collisional deactivation is seen in the second term in the denominator of the integral in Eq. (2a) (note the similarity to Eq. (1)). Since the temperature dependence of the collisional term is usually much less than of the chemical term, for a given molecule the rate constant falls-off from its high-pressure limit more and more as temperature increases. On the other hand, for a fixed temperature the pressure fall-off is reduced for larger molecules. These have more degrees of freedom, which increases the density of states and makes k(e) smaller. At the high-pressure limit, it is easily shown that kuni in Eq. (2) reduces to

ð3Þ

where I and I are the moments of inertia for the adiabatic rotations in the transition state and in the unimolecular reactant.

Q þ1 kc ½Mbc Q 1Q 2

Z

1

qðeÞ expðe=kTÞde

ð5Þ

e¼E0 þhDEJ i

The above equations for kuni can easily be solved numerically, but reformulation in terms of the master equation leads to more accurate results and also allows the theory to be applied when multiple reaction channels must be considered. 2.2. The master equation for multiple-channel, multiple-well reaction systems The master equation describes the population distributions and/or concentrations of reactants, intermediates, and products as the chemical system evolves with time. In complex unimolecular reaction systems, a highly excited molecule can isomerize to any number of intermediate structures (local minima on the potential energy surface) and any or all of these intermediates may isomerize further and dissociate via any number of possible channels. A full ME formulation for such systems can be written down, but its solution is a difficult problem because of the numerical challenges and because critical physical information is usually lacking. As a result, approximations must be adopted that are especially important in two areas: the treatment of angular momentum conservation and the treatment of collisional energy transfer. These issues are addressed in the following sections. 2.2.1. The two-dimensional master equation The time-dependent population distributions are functions of total energy and angular momentum, which are in principle described by a 2-D master Eq. (2DME). The 2DME is almost always solved by first reducing it to a 1-D master Eq. (1DME), which depends only on energy and which can be solved by any of several

606

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

methods. The 2DME can be written as follows for each chemical species in the reaction system [10,11,15,16]:

XZ 1 dNðE0 ; J 0 ; tÞ ¼ FðE0 ; J 0 ; tÞ þ RðE0 ; J 0 ; E; JÞNðE; J; tÞdE dt 0 J XZ 1  RðE; J; E0 ; J 0 ÞNðE0 ; J 0 ; tÞdE 0

J



channels X

ki ðE0 ; J 0 ÞNðE0 ; J 0 ; tÞ

ð6Þ

i¼1

In Eq. (6), N(E0 , J0 , t)dE0 is the time-dependent concentration of a chemical species with energy in the range E0 to E0 + dE0 and angular momentum quantum number J0 ; R(E, J; E0 , J0 ) is the (pseudo-first-order) rate coefficient for collisional energy transfer from initial energy E0 to energy E and quantum number J0 to J; F(E0 , J0 , t)dE0 is a source term (e.g. thermal, chemical, or photo activation, or isomerization); and ki(E0 , J0 ) is the unimolecular reaction rate constant for molecules at energy E0 reacting via the ith reaction channel. Other terms, such as those involving radiative emission and absorption, have been omitted because they are of little importance in combustion. The energy transfer kernel R(E, J; E0 , J0 ) is conventionally written as the product of the bimolecular collision rate constant kc(E0 , J0 ), the bath gas concentration [M], and the ‘‘collision step-size distribution”, P(E, J; E0 , J0 ), which expresses the probability density that a molecule initially with initial energy E0 will undergo an inelastic transition to the energy range E to E + dE, while simultaneously undergoing the transition from J0 to J: 0

0

RðE; J; E ; J ÞdE ¼

XZ J

0

1 0

(

)   R E; J; E0 ; J 0 dE R P RðE; J; E ; J ÞdE ð7aÞ 1 0 0 J 0 RðE; J; E ; J ÞdE 0

0

¼ kc ðE ; J Þ½MPðE; J; E0 ; J 0 ÞdE

0

ð7bÞ

The first factor on the right hand side of Eq. (7a), the integral over the rates of all inelastic transitions from the initial (E0 , J0 ) point, is the total frequency of inelastic collisions, kc(E0 , J0 ) [M], and the second factor (in curly brackets) is P(E, J; E0 , J0 )dE. Note that P(E, J; E0 , J0 ) is normalized:

XZ J

1

PðE; J; E0 ; J 0 ÞdE ¼ 1

ð8Þ

0

It is important to emphasize that the factorization of R(E, J; E0 , J0 ) in Eq. (2) is merely for convenience and that kc(E0 , J0 ) [M] and P(E, J; E0 , J0 ) never occur independently of one another. By considering detailed balance, we obtain the relationship between the probability densities for up- and down-collisions:

  PðE; J; E0 ; J 0 Þ kc ðE; JÞ qðE; JÞ ¼ exp ðE  E0 Þ=kB T PðE0 ; J 0 ; E; JÞ kc ðE0 ; J 0 Þ qðE0 ; J 0 Þ

ð9Þ

It is common practice to construct a normalized collision step-size distribution (the probability density) by specifying an arbitrary but physically plausible non-normalized function f(E, J; E0 , J0 ) and normalizing:

PðE; J; E0 ; J 0 Þ ¼

for up-steps instead. The most common choice in 1DMEs is the familiar exponential-down energy transfer model, which will be described later. Instead of applying detailed balance in a separate step, it can also be incorporated directly in the model function if desired [32]. Up to this point, no serious approximations have been introduced, but numerical solution of the explicit 2DME is very challenging and extensive ancillary information is needed before quantitatively accurate solutions can be obtained. For example, the collision rate constants and the function f(E, J; E0 , J0 ) must be supplied, as well as the densities of states and microcanonical rate constants. Because the requisite ancillary information is not known with certainty, several pragmatic assumptions have been adopted in the past. Over time, some of these assumptions have been replaced with better approximations based on specific knowledge. For example, the once-common rigid rotor and harmonic oscillator approximations have been eliminated in recent years by the development of much better methods for treating anharmonic and nonseparable vibrations [33–35]. Hindered internal rotors are now routinely treated by computing the exact eigenstates, although they are usually still treated as separable from other (coupled) degrees of freedom [35]. It has been conventionally assumed that the bimolecular collision rate constant kc(E, J) is independent of E, but recent work suggests that an energy-dependent version produces more accurate results, especially at low energies [29]. Developments in statistical reaction rate theories have enabled researchers to predict k(E, J) with a useful level of accuracy [17–19]. Promising treatments of non-statistical effects are still in their infancy [36– 39], but such effects are relatively small for many reactions and can often be neglected.

f ðE; J; E0 ; J 0 Þ NðE0 ; J 0 Þ

ð10Þ

where N(E0 , J0 ) is the normalization constant. After incorporating Eq. (10), the detailed balance expression Eq. (9) takes the following form:

  f ðE; J; E0 ; J 0 Þ NðE0 ; J 0 Þ kc ðE; JÞ qðE; JÞ ¼ exp ðE  E0 Þ=kB T NðE; JÞ kc ðE0 ; J 0 Þ qðE0 ; J 0 Þ f ðE0 ; J 0 ; E; JÞ

ð11Þ

For convenience, the non-normalized function f(E, J; E0 , J0 ) is usually specified for down-steps, but one could choose to specify a function

2.2.2. Reduction to one-dimension Because of the difficulties encountered in solving the 2DME, it is usually reduced to 1-dimension (E) by making simplifying assumptions, and then solving the 1DME numerically by one of the standard methods [15,16,40,41]. The first approach to this problem originated with Marcus [42], who invoked the concept of adiabatic rotation to make ‘‘centrifugal corrections”, as described above. This topic has been reviewed by Waage and Rabinovitch [43] and is found in standard monographs [10–12]. In this approach, the RRKM specific rate constant k(E, J) has been reduced to the approximate expression given by Eq. (2b), with DEJ given by Eq. (3), and all explicit reference to angular momentum quantum number J has been removed. By invoking the strong-collision assumption, which assumes that the Boltzmann distribution is always maintained, k(E, J) is reduced to k(E). A similar approach can be used to reduce the 2DME to a 1DME. Centrifugal corrections are very useful for pragmatic semiempirical models, where the expense of elaborate predictive calculations is not warranted. Long experience has demonstrated that this approach can give reasonable agreement with experimental data. For predictive calculations, however, practical solutions to the 2DME have been sought that can take better advantage of the improved level of skill in predicting k(E, J) for certain classes of reactions. Along these lines, an approach to reducing the 2DME to a 1DME was found by Smith and Gilbert [11,44] and extended by Miller, Klippenstein, and coworkers [16,40]. In place of the independent variables E, J, Smith and Gilbert expressed the 2DME in terms of e, J, where the ‘‘active energy”, e, of the energized species is defined as e = E  EJ, where EJ = BeJ(J + 1) and Be is the rotational constant of the species at its equilibrium geometry (see Fig. 1). The 2DME can be written equivalently in terms of either set of independent variables: E, J or e, J. Smith and Gilbert sought a steady-state solution of the master equation in order to obtain an expression for kuni, the unimolecular

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

rate constant. They discovered a simplifying assumption: if the energy transfer rate is expressed as the separable product of functions for total energy transfer and for angular momentum transfer, the 2DME can be reduced approximately to a 1DME with an effective value for k(e) that incorporates information about angular momentum (see their paper and book for more details [11,44]). Miller, Klippenstein, and Raffy [16,40] modified and extended this approach and made comparisons with a few sets of experimental rate data [16,40]. Their results showed that both the Smith–Gilbert approach and their extended version could give reasonable agreement with experiment, but their extended version seemed to perform slightly better. In both cases, the key assumption is that the rotational energy distribution after a collision is independent of the initial rotational energy. Although the 1DME obtained by applying this approach has had good success in empirically modeling experimental reaction rate data, it should be emphasized that the key assumption has not been verified by comparisons with direct energy transfer measurements or with explicit 2DME calculations. In fact, recent trajectory calculations on model systems strongly suggest that the key assumption is not accurate [35]. Clearly, this is still an active area of fundamental research. In summary, neither method for reducing a 2DME to a 1DME is exact, but both methods result in 1DMEs that give reasonable agreement with experimental data. The ‘‘active energy” (e) master equation is functionally similar to the total energy (E) master equation, but the density of states and specific rate constants differ in detail. Because of these differences, the energy transfer parameters (see below), which are usually adjusted empirically in order to fit experimental data, are different for the two approaches even when fitting the same experimental data. This difference prevents easy transferability of energy transfer parameters among calculation methods and can lead to confusion! After reduction to one-dimension, the master equation and detailed balance expression take the following forms:

Z 1 dNðx0 ; tÞ ¼ Fðx0 ; tÞ þ Rðx0 ; xÞNðx; tÞdx dt 0 Z 1 channels X  Rðx; x0 ÞNðx0 ; tÞdx  ki ðx0 ÞNðx0 ; tÞ 0

ð12Þ

i¼1

f ðx; x0 Þ Nðx0 Þ kc ðxÞ qðxÞ ¼ exp ½ðx  x0 Þ=kB T  f ðx0 ; xÞ NðxÞ kc ðx0 Þ qðx0 Þ

ð13Þ

where x designates either E (the total energy), or e (the active energy) and the primed and un-primed quantities correspond to the initial and final values, respectively. The other symbols retain the definitions given earlier. The solution of this integro-differential master equation with the detailed balance constraint depends on the particular application, but ancillary information is always required for densities of states, reaction rate constants, thermochemistry, and energy transfer. 2.2.3. Solving the 1-D master equation The 1DME can be solved in a number of ways, depending on the application. Before turning to specific solutions however, it is useful to consider the qualitative physical process for a single-channel reaction system. Depending on the application, one may assume an initial finite quantity of reactant, or one may assume that the reactant is continuously supplied (i.e. F(x0 ; t) = 0). Initially (at t = 0), the reactant is assumed to be present in a specified initial energy distribution. The MultiWell Program Suite [27,29,35], for example, allows users to choose among a delta function at a specified energy, a thermal Boltzmann distribution at a specified temperature, a chemical activation distribution produced at a specified temperature and from a specified reaction, or a user-supplied energy distribution.

607

As the process begins, both chemical reaction and collisional energy transfer can take place simultaneously. For reactants with low internal energy, collisional energy transfer is usually much more rapid than chemical reaction and the energy transfer proceeds toward the establishment of a steady-state energy distribution. In the absence of reaction, the steady-state distribution is just the thermal Boltzmann distribution. Reactants with energies higher than the reaction critical (threshold) energy E0 may react more rapidly than collisional thermalization, resulting in a depletion of population at high energies. Eventually, a steady-state energy distribution is established in which the rates of chemical reaction and of collisional activation and deactivation are balanced at all energies. Prior to establishment of the steady-state, the energy distribution evolves with time and it is not possible to define a proper rate constant (which must be independent of time). Only after the fast transients have subsided and steady-state has been established is it possible to extract a rate constant. In many cases of practical interest, steady-state is never achieved, or the approach to steady-state is itself of interest. Such is the case, for example, in measurements of the ‘‘incubation time” in shock tube experiments (see below). When steady-state cannot be assumed, the time-dependent 1DME (Eq. (12)) must be solved. Except for special, idealized models, the 1DME cannot be solved in closed form and one must resort to numerical solutions. A common point of departure is to formulate the master equation as an energy grained master equation. 2.2.3.1. Energy-grained master equation. Historically, the 1DME is solved numerically by first recasting it in the form of an ‘‘energy grained master equation” (EGME). The EGME is constructed by dividing the energy scale into a large number of ‘‘bins” or ‘‘energy grains” and expressing the 1DME as a finite set of coupled ordinary differential equations describing the rate of change of the population in each grain. When assuming steady-state, the coupled differential equations reduce to a set of coupled algebraic equations, one equation for each energy grain. In principle, the solutions become independent of grain size and converge to the exact solutions when the energy grains are made small enough, but convergence is often thwarted because computer time and numerical difficulties increase as the number of energy grains increases. The EGME can be written in matrix notation [45]:

dy ¼ Jy dt

ð14Þ

where y is a column vector of the time-dependent populations and the square matrix J (described as negative semi-definite) includes energy transfer and reaction [46]. The solution can be found using an eigenvector-eigenvalue analysis: one eigenvalue is identically equal to zero and all of the others are <0 (i.e negative). When the eigenvalues are well separated, the non-zero eigenvalue with the smallest magnitude is identified with the reaction rate constant. Whether or not the eigenvalues are well separated, the solution of the master equation can be expressed as an eigenvalue expansion. The details of the eigenvector-eigenvalue analysis can be found elsewhere [8,10–12,47–49]. The steady-state master equation is appropriate for computing rate constants under steady-state conditions. It is implicitly assumed that the population distributions achieve steady-state on a time scale that is very short relative to that for chemical reaction. If this condition is not satisfied, the rate constants cannot be defined unambiguously and the rate equations will not produce good simulations under the experimental conditions. Good separation of time scales is achieved for relatively slow reactions, but caution is needed when reactions are proceeding at nearly the same rate as collisional activation and deactivation. Klippenstein and Miller

608

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

have described a useful method for obtaining phenomenological rate constants from the eigenvalue-eigenvector solution of the time-dependent multi-well master equation when the eigenvalues are well separated [50].

simulation of the experiment. In other words, the stochastic simulation is essentially identical to the numerical solution of a standard energy-grained time-dependent master equation. 2.3. Factors in the 1-D master equation

2.2.3.2. Stochastic simulation. Of greatest interest to us is the timedependent solution of the master equation. Instead of utilizing the master equation directly, it is possible to obtain time-dependent solutions by carrying out stochastic simulations. As long as the stochastic simulations are based on the same physics as the master equation, the results are identical in the limit of very small grain size (for the master equation) and a large number of stochastic trials (for the stochastic simulation) [24,25,41,51]. Chemical reaction systems are intrinsically stochastic. This is because any finite volume contains an integer number of particles and because chemical reactions occur more or less randomly within the volume. The usual deterministic rate equation expressions of chemical kinetics are actually approximations that are accurate for large numbers of particles, but perform poorly when the numbers are small. A great strength of stochastic methods is that extremely complicated systems can be simulated in a straightforward way. Coupled sets of stiff non-linear differential equations, for example, often require special techniques for solution, but stochastic simulations are in no way special. In the most recent implementations, highly complicated multi-well, multi-channel unimolecular reaction systems can be modeled with relative ease [13,27,35,41,52]. A weakness of stochastic simulation is that it requires significant computer time. The results of a stochastic simulation are ‘‘noisy”. In the limit of an infinite number of stochastic trials the results converge to the exact solution of the master equation, but for a finite number N of stochastic trials, the stochastic error (relative standard deviation) in the results is proportional to N1/2. Thus many trials are needed for precise results, as discussed elsewhere [27]. The application of Gillespie’s stochastic simulation algorithm (SSA) [24,25,51] to master equations has been described elsewhere [13,26,27,41]. It is the basis for the MultiWell master equation code [35]. The basic idea for applying Gillespie’s SSA is to carry out a stochastic random walk as a function of energy. The time duration of each step is determined by Monte Carlo selection from the proper probability density distribution of time steps. Monte Carlo selection is also used in each step to choose the outcome: collisional activation, collisional deactivation, or reaction via a selected path. The initial energy distribution and the size of every individual energy transfer step are also determined by Monte Carlo selection. The MultiWell master equation code implements this algorithm in a way that maintains detailed balance, allows treatment of complicated multi-well reaction schemes, offers multiple choices for collisional energy transfer models, enables calculation of k(E) microcanonical rate constants, and a long list of other capabilities [27,35]. In particular, it enables simulations of unimolecular reaction systems that involve thermal activation, photo activation, and chemical activation. A list of literature citations can be found on the MultiWell web site, where the free-source computer codes are available for down-loading. The results of a MultiWell stochastic simulation are presented in a group of output files. The central result is a set of time-dependent fractional abundances for all of the chemical species. After the reaction has gone to completion, these ‘‘fractions” correspond to chemical branching ratios. The time-dependent variation of the abundances can be used to extract phenomenological rate constants when steady-state population distributions have been achieved. For simple recombination or unimolecular reactions, the abundances are used to extract the corresponding rate constants. Otherwise, the fractions constitute a direct time-dependent

Whether the problem is solved by using an eigenvalue-eigenvector solution of the energy grained master equation or by stochastic simulation, the chemical physics of the problem is the same. In this section, we summarize the principal factors that must be included in an ‘‘active energy” one-dimensional (e) master equation, which is the basis for the MultiWell master equation code. 2.3.1. Sums and densities of states: G(e) and q(e) Densities of states q(e) are commonly computed by assuming that the vibrations, rotations, and internal rotations are separable. Stein and Rabinovitch [34] showed how the Beyer–Swinehart algorithm [33] can be used to carry out efficient ‘‘exact counts” of states for separable degrees of freedom. Recent developments have produced practical methods for computing sums and densities of states for non-separable degrees of freedom [35,53,54]. The new algorithms produce results that are more realistic, but at the expense of more time-consuming calculations. The ‘‘active” degrees of freedom include the vibrations, internal rotations, and the ‘‘K-rotor”, which is the rotation associated with the symmetry axis (the K quantum number) of a symmetric top. Most chemical species can be approximated as symmetric tops, even if they are unsymmetrical [10]. It is assumed that energy can randomize freely among the active degrees of freedom in both the excited unimolecular reactant and the transition state. As discussed above, the ‘‘adiabatic” degrees of freedom are associated with quantum number J, which is assumed to be conserved [10] and which is used in the approximate centrifugal corrections described above. For a rigid rotor, the range of quantum number K is restricted to the 2J + 1 values between J and +J. However, the usual pragmatic approximation is that K is restricted only by conservation of energy and its range is not otherwise restricted [10]. Aubanel et al. have commented that because real chemical species are not rigid rotors, this pragmatic approximation is possibly more realistic than restricting the range of K [55]. Thus this pragmatic approximation is a standard practice. 2.3.2. Microcanonical rate constants: k(e) Rate constants are usually computed from RRKM theory, as described above, or by using the inverse Laplace transform (ILT) of the high pressure limit rate constant. The latter technique is specifically intended for semi-empirical applications, since it incorporates the Arrhenius parameters of the k1, the unimolecular rate constant at the high-pressure limit [10–12,49]. The ILT rate constant is written

kILT ðeÞ ¼ A1

qðe  E1 Þ qðeÞ

ð15Þ

where A1 and E1 are the Arrhenius parameters and q(e) is the density of states of the reactant (including the K-rotor); q(e  E1) is also the density of states of the reactant, but shifted in energy by an amount equal to E1. The ILT method has also been extended to non-Arrhenius rate constants [56,57]. Various strategies can be used to estimate k(e) from computed or experimental parameters. When the transition state is associated with an intrinsic energy barrier (i.e. a local maximum on the PES (including zero point energy)), electronic structure calculations can be used to compute the vibrational frequencies, and structure of the transition state. These parameters can be used to

609

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

compute the sums and densities of states directly. When no local maximum exists on the PES (i.e. a ‘‘barrier-less reaction”), as in many radical + radical recombination reactions, variational transition state theory (VTST) may be used to obtain the rate constant. Microcanonical VTST can be used to compute k(e, J) (or k(E, J)) directly by means of a suitably flexible model for the reaction [17,19,58–63]. From k(E, J), one can use the Smith–Gilbert [44] or Miller et al. [16,40] simplifying assumptions. Alternatively, a common approximation is to assume that the distribution of J is thermally distributed and k(e) is the thermal average (k(e) = hk(e, J)i) or corresponds to the average J value (k(e) = k(e, hJi)). It is easier to use canonical VTST to compute k1(T) as a function of temperature. From the computed k1(T) it is not difficult to obtain Arrhenius parameters (or more complex functions), which can be used with the ILT method to obtain k(e). The energies calculated using electronic structure methods are in most cases subject to significant uncertainties, which have a strong effect on reaction rate constants at low temperatures; the effect is not as large at combustion temperatures. Thus a computed value for E1 is often adjusted empirically to fit experimental data for k1(T). Although E1 is not identical to E0, the critical energy for reaction, replacement of E1 by E0 in the ILT expression can lead to more accurate k(e) near the reaction threshold. One can obtain a semi-empirical estimate of k(e) by using computed or experimental k1(T) in combination with simple models. The hindered, or restricted, Gorin model is a prime example. 2.3.2.1. Hindered Gorin model. When the high-pressure rate constant is known, it is more convenient to use a ‘‘restricted” Gorin Model with a ‘‘hindrance parameter” selected to reproduce the known rate constant [64]. According to the Gorin model, the two molecular fragments rotate independently of one another while separated at the distance corresponding to the centrifugal maximum (rmax) of the effective potential of the bond being broken [65]. For those systems with loose transition states, the rotations of both fragments and the over-all transition state rotations are treated approximately as symmetric tops. The over-all transition state is assumed to have a 2-D external adiabatic rotation with moment of inertia given by Iz2D ¼ lr 2max , where l is the reduced mass of the two fragments, and a 1-D external rotation (the ‘‘K-rotor”) with moment of inertia Ik. The K-rotor is not adiabatic and is assumed, as discussed above, to mix energy freely with the active vibrations. When non-linear, the fragments A and B, which make up the transition state, are each characterized by 2-D rotations with moments of inertia Ia and Ib, respectively. In addition, the remaining rotation in each fragment is combined in the transition state to produce an internal rotation (with reduced moment of inertia Ir), corresponding to counterrotation of the two fragments. Given the uncertainties in many quantities, the K-rotor moment of inertia for the transition state is often assumed to be identical to that of the unexcited adduct molecule. The Gorin model is essentially the same as Phase Space Theory. The restricted or hindered-Gorin model [66,67] is a semi-empirical model that is parameterized by fitting to experimental data. According to this model, it is assumed that the two fragments interfere sterically with each other and thus cannot rotate freely. The effect is to reduce the available phase space and hence reduce the sum of states. Operationally, a ‘‘hindrance” parameter g is defined, which can vary from zero (free rotation) to unity (completely hindered). The 2-D moments of inertia Ia and Ib (for fragments A and B) are multiplied by the factor (1  g)1/2 to obtain the effective 2-D moments of inertia used for calculating the sum of states. Rate constants are then computed and the hindrance adjusted until agreement is achieved. The hindrance parameter can be computed at each temperature and separate calculations are performed at

that temperature. All the examples herein were performed that way. The fitted model enables the extrapolation and interpolation of the limited experimental data. In particular, it provides an estimate for k1 at each temperature. For many reactions of practical import, the potential function describing the breaking bond is not known. For simplicity and because it has the long range dependence on r6 expected for the long range potential between neutral non-polar species, the Lennard–Jones potential [68] has often been assumed. For the Lennard–Jones potential, the moment of inertia for the two-dimensional adiabatic external rotation is given by Iz2D ¼ lr 2e ð6De =RTÞ1=3 , where re is the equilibrium bond distance, l is the reduced mass, and De = D0  DEz, where D0 is the bond dissociation enthalpy at 0 K, DEz is the zero point energy difference between products and reactants, and R is the gas law constant. The Morse [69] and Varshni [70] potential functions are more accurate near the equilibrium bond distance. An alternative calculation of Iz2D uses the distance between the bonding atoms and computes the moments of inertia for this transition state. If another potential is used in place of the Lennard–Jones potential, the empirical fitting parameters needed for the restricted Gorin model are somewhat different, but the final results are essentially unchanged, since whatever model is chosen must be parameterized by fitting to the same experiments. 2.3.3. Collision rate constant for energy transfer: kc The bimolecular collision rate constant kc(x0 ) (see Eq. (7b)) is usually assumed to be independent of energy and to be given by kLJ, the rate constant for collisions between particles that interact according to a Lennard–Jones intermolecular potential[10–12] Recent work, however, has shown that more accurate results can be obtained when an energy-dependent rate constant is assumed [29]. This new approach has been incorporated in MultiWell. In this approach, it is assumed that kc(x0 ) = kLJ at a reference energy (see Ref. [29] for details). Thus the Lennard–Jones parameters are required as input parameters for computing kLJ. 2.3.4. Energy transfer model Energy transfer continues to be an area of significant uncertainty [21]. For practical calculations, the exponential-down model [7,8,71–73] is usually employed, unless there is specific information that enables the use of a biexponential model [74,75] or of Luther’s model, which is based on kinetically controlled selective ionization (KCSI) experiments [76–80]. Luther’s model applies to deactivating collisions (energy transfer in activating collisions is obtained by using detailed balance):

c 

 0 ðe  eÞ for f ðe; e0 Þ ¼ exp  0 aðe Þ

e 6 e0

ð16Þ

where a(e0 ) is an energy-dependent parameter which is usually either assumed to be a constant, or is expressed by a low-order polynomial obtained from experimental data [21]. Luther’s model allows the parameter c to vary. Typical values of c range from about 0.6 (for weak colliders) to 1.5 (for strong colliders). The exponentialdown model corresponds to c = 1. The energy transfer parameter a(e0 ) is closely related to hDEidown, the average energy transferred in deactivation collisions: 0

hDEðe Þidown ¼



h 0 ic ðe  e0 Þ exp  ðae ðe0eÞÞ de h 0 ic R e0 exp  ðae ðe0eÞÞ de 0

R e0 0

ð17Þ

hDEidown must be distinguished from hDEiall, which is obtained by averaging over both up- and down-collisions: 0

hDEðe Þiall ¼

R1 0

ðe  e0 ÞPðe; e0 Þde R1 Pðe; e0 Þde 0

ð18Þ

610

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

where P(e, e0 ) is the step-size distribution. The parameter a(e0 ) is usually adjusted empirically to fit experimental data, but it is sometimes estimated by analogy with other energy transfer collision partners for which experimental data are available. Classical trajectory calculations have proved to be quite effective for modeling energy transfer involving large molecules [81– 96], but the requisite potential energy functions must be obtained using electronic structure calculations, which are often not accurate enough for predictive purposes. A 2-D energy transfer probability density was recently obtained by using classical trajectories [32]. 2.3.5. Rate constants for radiationless electronic transitions Chemical species with different spin quantum numbers and in different electronic states can easily be defined in ME models. Intersystem crossing (ISC) can occur when PESs of differing spin (e.g. singlet PES and triplet PES) intersect one another. Internal conversion (IC) can take place via conical intersections between PESs of the same spin (e.g. between a ground state singlet PES and an excited state singlet PES). It is not trivial to compute the rate constants for ISC and IC from first principles, but sometimes suitable rate constants are known from experiments. It is beyond the scope of this paper to discuss theoretical methods for computing rate constants for radiationless transitions and reader is referred to a few papers which may be of assistance [97–99]. Good experimental data exists for ISC in pyrazine [100–104]. In this series of papers, collisional energy transfer was investigated by monitoring the decay of phosphorescence from pyrazine initially excited to the triplet state. The ISC rate is slower at lower energies and thus the phosphorescence decay time lengthens as the vibrationally excited triplet state molecules undergo vibrational deactivation. Although never published, stochastic ME simulation was successful in modeling the competition between ISC and collisional deactivation [L.A. Yoder, J.R. Barker, unpublished (2000)]. The calculations were carried out by defining kISC(E), the energy-dependent rate constant for ISC from the vibrationally excited triplet state molecules to the ground singlet state. The energy dependence of kISC(E) was obtained by least-squares analysis of the experimental ISC rate constants (functions of excitation wavelength) and collisional energy transfer was described by using the conventional exponential model with energy transfer parameter a(E0 ). By adjusting a(E0 ) it was possible to simulate the experiments reported by Weisman and coworkers [100–104]. 2.3.6. Initial energy distributions Master equation simulations are initiated by specifying the chemical species that is initially excited (in a multi-well system) and the initial (active) energy distribution. Several standard initial energy distributions are described in this section. All of them are available as selectable options in the MultiWell master equation code [13,27,35].

2.3.6.2. Photo-activation. For photoactivation, it is often assumed that the absorption cross section is independent of any initial vibrational excitation residing in the absorber. Thus if the absorber is initially at thermal equilibrium, the population of photoexcited molecules will consist of the initial Boltzmann distribution, shifted to higher energy by an amount equal to the photon energy ht:

Photo-activation : P phot ðeÞ

0

qðeÞ exp½e=kT qðeÞ exp½e=kTde

Thermal :BðeÞ ¼ R 1 0

ð19Þ

where T is the ambient (bath) temperature. For shock tube simulations, T is the ambient temperature prior to arrival of the shock (e.g. 300 K); after passage of the shock, the translational temperature may be much higher (e.g. 1500 K). In MultiWell, two temperatures are used as input data: Tvib, the equilibrium vibrational temperature prior to arrival of the shock and Ttrans, the translational temperature after the shock has arrived.

ð20Þ

2.3.6.3. Chemical activation. In recombination reactions, as in all other adduct-forming reactions, the nascent distribution of a new product species is located entirely above its dissociation energy. It can re-dissociate via the association pathway. Since the reactants are assumed to have been at thermal equilibrium prior to reaction, they carry with them some thermal energy, which will reside in the newly formed adduct. The chemical activation distribution describes the nascent reactive flux originating from reactants. By the principle of detailed balance, it is also the energy distribution of the reactive flux originating from dissociation of the adduct, when it is in equilibrium with the reactants. Thus it can be expressed in terms of either the forward reaction or the reverse. For convenience, it is usually expressed in terms of the dissociation reaction:

kðeÞqðeÞ exp½e=kT kðeÞqðeÞ exp½e=kTde E0r

Chemical activation : Pca ðeÞ ¼ R 1

ð21Þ

where k(e) is the microcanonical rate constant for dissociation, q(e) is the density of states of the adduct, and the range of integration for barrier-less reactions is from the critical energy E0r for the dissociation reaction (modified by centrifugal corrections) to 1, as shown. Quantum mechanical tunneling can be important when there is an intrinsic energy barrier. In that case, the range of integration is from the thermochemical energy of the reactants (higher than the adduct energy but lower than the top of the intrinsic energy barrier) to e = 1. 2.4. Empirical formulas It has been the practice to codify the results of master equation calculations using variations of semi-empirical expressions devised by Troe [105,106].

kðM; TÞ

 1 k0 ðTÞ½M f1þ½log ðk0 ðTÞ½M=k1 ðTÞÞ=ð0:751:27 logðF c ÞÞ2 g ð22Þ Fc ¼ 1 þ ðk0 ðTÞ½M=k1 ðTÞÞ With the rate constants k0(T) and k1(T), represented as functions of temperature in modified Arrhenius form, Arrhenius form or as exponential functions of T. The parameter Fc is represented by equations of the form

F c ¼ ð1  aÞ expðT=T  Þ þ a expðT=T  Þ þ expðT  =TÞ 

2.3.6.1. Thermal activation. For a thermal unimolecular reaction, one must specify a Boltzmann probability density function:

qðe  htÞ exp ½ðe  htÞ=kT  qðe  htÞ exp ½ðe  htÞ=kT de

¼ R1



ð23Þ



where a, T , T and T are constants. Data and evaluations in this form for combustion are found in Baulch et al. [107] to which readers are directed for many references therein. The IUPAC Evaluation [108] uses the equation above, but treats Fc as temperature-independent, but varying for different reactions. For the NASA/JPL Evaluation [109], values of k0 and k1 were chosen to best describe the data according to:

kðM; TÞ ¼

 2 1 k0 ðTÞ½M 0:6f1þ½log ðk0 ðTÞ½M=k1 ðTÞÞ g 1 þ ðk0 ðTÞ½M=k1 ðTÞÞ

ð24Þ

Both the evaluations usually describe extant data adequately. They can and do disagree when extrapolated out of the measured data range.

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

These expressions are only useful for single channel reactions and cannot be used for multiple-channel or multiple-well systems. In particular, chemical activation processes, wherein a higher energy recombination channel produces an intermediate with a reactive pathway at lower energy, require the master equation approach. Look-up tables may be the best way to represent such processes in a pragmatic combustion chemical model. 3. Example applications of master equation simulations In order to create chemical models for practical systems such as the chemistry of combustion or of the atmosphere, it is necessary to list the relevant fundamental chemical steps, their thermochemistry and their rate coefficients. For some reactions, this latter information must include not only values as a function of temperature, but as a function of pressure and nature of possible colliding species as well. Typically reactions that are truly bimolecular may be represented as a function of temperature by an Arrhenius equation. (It is best to use at least three parameters, but often only two are used.) In the more complex case of unimolecular reactions and association reactions or chemical activated processes, parameterization is accomplished by using a ME method, fitting to a data set (which is often limited in comparison with practical needs) and extrapolating to the regions of interest. Current practice usually involves using some form of RRKM theory and standard assumptions about energy transfer. It is often the case, that when comparing calculated results to the data that the choice of the high-pressure-limiting rate constant determines the fit at higher pressures and the choice of energy transfer parameters determines it at the lower end of the pressure scale. The current implementations of rate theory, enhanced as it is by quantum calculations of molecular and transition state parameters, yields reasonable values for the highpressure rate constant, while the exponential-down collisional energy transfer a is fitted using experimental data. Often, it is found that the empirical function a(T) = a300(T/300 K)c provides a satisfactory fit. The model also entails use of the Lennard–Jones collision frequency, the parameters of which are often estimated. (The Multiwell Manual offers a list of such parameters.) The energy transfer parameter a is known to vary with the nature of the collision partners [21,110]. In the following, several examples pertaining to both combustion and atmospherically important processes are presented; all were addressed using the Multiwell suite of codes. 3.1. Methane formation/dissociation An enormous number of papers have been devoted to various aspects of methane dissociation and the reverse. Most are referenced in a review by Baulch et al. [107], which gives recommendations that are largely based on studies due to Cobos and Troe [107,111]. Harding et al. [112] presents a theoretical study of the high-pressure limit. Methane being the simplest hydrocarbon fuel, any models describing its combustion require adequate knowledge of the rate parameters for this system. A theoretically consistent codification that includes the rate parameters as functions of temperature, pressure and nature of the bath gas for the association reaction along with the parameters for the temperature dependence of the equilibrium constant was presented by Golden [113]. To avoid potentially serious inconsistencies among the rate constants and the equilibrium constant, it is strongly recommended that rate parameters for the association rate constant be used along with the equilibrium constant based on thermochemical data (the dissociation reaction rate constant is calculated using the association rate constant and the equilibrium constant). Given

611

the uncertainties in both quantities required for theoretical calculations and those inherent in an experimental study, it is strongly suggested that the results of this codification, along with the appropriate uncertainty limits, are to be taken as a starting point for optimization in a given model such as exemplified in the study known as ‘‘GRI-Mech” [3]. Given current capabilities using codes such as Multiwell [23,27] to compute rate constants from fundamental information, this codification can be in the form of lookup tables that allow interpolation, without introducing errors due to fitting the results to empirical expressions, although the latter may be simpler to use On-the-fly calculations using MultiWell or some other master equation code are also feasible. Golden [113] used the hindered-Gorin model for the transition state and the parameters in Table A1 (Appendix), to reconcile most of the extant data on the decomposition of methane and the reverse combination of a methyl radical with a hydrogen atom. Using the values in Table A1 for methane and the critical energy and literature values for structural parameters for CH3 and H, the equilibrium constant can be calculated. Table A2 shows that even for these well-known species computed values of the equilibrium constant can vary by more than 20%. Table A2 also illustrates some inconsistencies in the literature concerning values for the equilibrium constant obtained from ratios of the rate constants at limiting pressures. As is almost always the case in these types of calculations, the value of DEd(T), needed for the description of the energy transfer, is used as a fitting parameter, and anharmonicities are neglected in the calculation of densities and sums of states as well. The uncertainties that arise from all the above combined with experimental uncertainties is part of the argument for optimization and data collaboration advanced earlier. In this case, it was found that a very large segment of the extant data could be fit using the master equation/RRKM method embodied in the Multiwell suite [23,27], with values of the energy transfer parameter that varied with temperature according to a(T) = a300(T/300). Values at 300 K were 100, 150, 400, 500 cm1 for He, Ar, CH4, C2H6 respectively, in good agreement with the expected ratios [11,114]. 3.2. Methyl chloride formation/dissociation Any examination of the literature will discover values for DEd that vary over wide ranges. In fact in a series of studies, where the object was to fit pressure dependent data in atmospherically relevant temperature ranges, for the combination of Cl + NO2[115], BrO + NO2[116], IO + NO2[117] and ClO + ClO [118], it was found that in some cases even very high values of DEd = 5000 cm1 were not sufficient to explain the low-pressure end of the data curves. However, in the study [113] of methane decomposition and/or the combination of CH3 + H to give methane discussed above, it was found that ‘‘rational” values could be used to fit data in He, Ar, CH4, and C2H6 over the temperature range 200 < T/K < 2000. In that study the model and its parameters are clearly stated (See Table A1). That is, a hindered-Gorin model is used to obtain the high-pressure limit and the usual exponential-down model is used to describe energy transfer. The hindered-Gorin model may be only an approximation, but there is no doubt that, as used, it can be made to reproduce the high-pressure-limiting rate coefficients obtained in more sophisticated studies [112]. Armed with this information and curious as to why the earlier mentioned studies were so different, the CH3Cl system was chosen as a logical next step to explore the effect of molecular structure on average energy transferred. (HCl elimination from CH3Cl was not considered, as the barrier is expected to be high and the A-factor low.) In contrast with the satisfyingly consistent results in the methane study [113], the data fit for the methyl chloride case [119] was disconcerting. For example, the value for DEd for the He collider in the association study is 400 cm1 for all three temperatures

612

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

studied. This compares with a value of 100 cm1 in the methane case. In addition, the usual assumption of Lennard–Jones collision frequency was augmented by 30%. Neither of these is necessarily incorrect, but they, together with fits to the higher temperature dissociation studies, which found (DEd)T = 4  (T/300)0.4 where the collision partner was Ar or CH4 rather than He, prompted calculations analogous to those carried out in the methane study [113]. The calculations for CH3Cl, using the parameters shown in Table A3 together with the Multiwell suite of codes [120], found that the value of energy transfer parameter a needed to describe the data for CH3 + Cl in He at temperatures between 200 and 300 K is about four times higher than the value needed to fit data for the CH3 + H reaction. The data for the dissociation reaction CH3Cl ? CH3 + Cl in Ar are also higher than the value needed to explain the data for the analogous reaction CH4 ? CH3 + H in Ar. While the parameters used in [120] may fit the extant data, extrapolation outside the experimental limits (a goal of all exercises of this type!) is more problematical, especially since the conservation of angular momentum links the parameters of the high-pressure limit to those of the low-pressure limit (see Eqs. (2b) and (3)). This means that the values of the energy transfer parameter are linked to the structure of the transition state, making an understanding of the values of DEd more difficult. This is a general problem that reflects the need for a deeper understanding of energy transfer. 3.3. Incubation in shock waves ”Incubation” is the non-steady-state process during which collisional thermalization takes place simultaneously with chemical reactions. When incubation is essentially instantaneous on the time scale of the reaction, the arrival of the shock wave merely starts the unimolecular reaction process. However, if significant reaction occurs while collisional thermalization is taking place, the reaction yields are affected and reaction rate constants cannot be defined unambiguously. Prior to passage of a shock wave, the active energy of a unimolecular reactant is described by the ambient thermal distribution (typically near 300 K). The adiabatic compression associated with passage of the shock wave produces a high translational temperature on a very short time scale. Translations are thermalized within a few collisions, rotations require a few more, and vibrations typically require a large number [121]. Thus the vibrational energy relaxes slowly, long after translations and rotations have been thermalized at the new temperature. In this context, ‘‘relaxation” refers to evolution of the energy distribution from the initial Boltzmann distribution at ambient temperature to the final steady-state distribution, which is subject to the shock temperature and to the possible influence of chemical reactions. The transient internal energy relaxation subsequent to the shock can be followed by monitoring the average vibrational energy and by measuring the time-dependent yield of the unimolecular reaction. As the relaxation takes place, the average energy increases until a new steady-state energy distribution has been established. The ‘‘vibrational relaxation time” (svib) characterizes the transition from the initial to the final energy distribution. In the absence of unimolecular reaction, the final energy distribution is a Boltzmann at the same temperature as the bath. When the unimolecular reaction rate can compete with collisional activation, the final steady-state energy distribution at higher energies is preferentially depleted, relative to the thermal distribution. This steady-state energy distribution is the familiar ”fall-off” distribution, which gives rise to the pressure fall-off of the unimolecular reaction rate coefficient [10–12]. While the relaxation is taking place, there is a time delay before the reactant molecules are activated above the reaction threshold

and unimolecular reaction can occur. The ‘‘incubation time” (sinc) is the delay between the passage of the shock and the onset of unimolecular reaction. Only after the final steady-state population distribution has been established can time-independent rate constants be defined unambiguously. Several shock-initiated reaction systems have been studied in order to better understand the incubation period and energy transfer. Kiefer and coworkers have developed laser schlieren methods with exceptionally fast time resolution for measuring sinc and svib in a shock tube [122–124] Three systems investigated by Kiefer and coworkers have been modeled by Barker and coworkers, who employed a linear master equation with numerical solutions obtained using Gillespie’s stochastic simulation algorithm (see above). These are the decomposition reactions of cyclohexene [123,125], norbornene [126,127], and 1,1,1-trifluoroethane [128, 129]. The linear master equation used for the shock tube simulations assumes the reactant molecule is present in very low concentrations in an inert bath gas and neglects possible ‘‘energy pooling” collisions between two excited reactant molecules. Non-linear master equation simulations of energy pooling have been reported by Davis et al. [130]. The norbornene decomposition reaction [126,127] provides a good illustration of the principles involved. This system was modeled originally by using a master equation code that was similar to MultiWell. The new approach to collisional energy transfer that was recently incorporated in MultiWell requires somewhat different empirical parameters for fitting the experimental data, but in the end the master equation results are practically indistinguishable [29]. The norbornene thermal decomposition yields 1,3-cyclopentadiene and ethylene as stable molecular products.

The unimolecular reaction rate data [126,131–133] are consistent with a retro-Diels-Alder elimination. Kiefer et al. studied the reaction over wide ranges of both temperature and pressure (542– 1480 K and 34–416 Torr of Kr bath gas). The master equation calculations were carried out using a transition state model, which was based in part on one reported of Kiefer et al. [126], but which was adjusted [127] to fit the experimental data for kuni. Because the reaction has a tight transition state, centrifugal corrections are minimal. The conventional exponential-down model was employed (Eq. (16) with c = 1). LennardJones parameters were taken from Kiefer et al. and the energy transfer parameter a(E) was fitted by using the experimental data for both kuni and sinc. Full details can be found in the original reference [127]. The empirical parameter a(E) was revised [29] for use with MultiWell (mentioned above) and is represented by a linear function:

aðEÞ ¼ 70 þ 0:003E

ð25Þ

where a and E are expressed in cm1 units. All of the parameters needed for the calculations are given in Appendix. In the model simulations, the translational temperature is assumed to increase instantaneously to the shock temperature. The vibrational energy distribution subsequently evolves due to collisions with the hot krypton bath gas, as shown by the simulation results in Fig. 2. At time t = 0, the vibrational distribution is the 300 K Boltzmann. After 30 ls, the vibrational energy distribution is approaching the new steady-state, corresponding to the fall-off distribution at 1197 K. The average vibrational energy hEvibi has

613

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

Probability Density (arbitrary)

1.0

Norbornene

Norbornene

T(initial) = 300 K T(shock) = 1197 K Kr = 48.3 Torr

N / N0

T(initial) = 300 K T(shock) = 1197 K Kr = 48.3 Torr

0 µs 2.5 µs

τ

inc

5 µs 10 µs

0

5000

25 µs

10000

Evib,

15000 cm-1

20000

0.1

25000

increased from 600 cm1 at 300 K to 14,000 cm1 (see Fig. 3) and the width of the distribution has greatly increased. The hEvibi for the steady-state distribution is slightly lower than that of the Boltzmann (hEvibi  15,500 cm1) because the chemical reaction selectively depletes the population at the highest energies. The increase in hEvibi can be fitted to an exponential, defining the vibrational relaxation time svib. As the high-energy portion of the energy distribution evolves to energies above the reaction threshold energy (15,000 cm1), reaction begins to occur and the concentration of norbornene decreases. When the steady-state vibrational energy distribution has been established after 30 ls the decay is first order: a straight line on the plot. The extrapolation of the straight line back to where N/N0 = 1 defines the incubation time sinc, as shown in Fig. 4. The most interesting part of this process is the energy relaxation and its effect on the reaction rate. Until steady-state has been established, the derivative d(ln N)/dt is not constant and therefore the rate constant cannot be defined. At higher pressures, sinc becomes smaller: at sufficiently high pressures, sinc becomes negligible and the rate constant can be defined over essentially all of the decay. Measurements of sinc and svib provide crucial information about the rates of energy transfer, as do measurements of pressure falloff of kuni. However, sinc and svib provide information about energy transfer at energies below the reaction energy threshold, while kuni is sensitive mostly to energies above the threshold. Information of this type is very useful for building accurate ME models, but the experiments are quite challenging.

< Evib >, cm-1

15000

10000

Norbornene

T(initial) = 300 K T(shock) = 1197 K Kr = 48.3 Torr

0

0

10

20

40

60

80

Time (µs)

Fig. 2. Evolution of the vibrational energy distribution following arrival of a shock.

5000

0

20

30

Time (µs) Fig. 3. Average vibrational energy during the evolution of the population distribution.

Fig. 4. Semi-log plot of the fraction of norbornene as a function of time following the shock. The Incubation time sinc is defined as the time at which the N/N0 = 1 for the extrapolated first order decay, as shown by the dashed lines.

3.4. Multiwell, multi-channel simulations Many recombination and association reactions produce highly rovibrationally excited species, which may isomerize to produce one or more isomers and may also decompose to produce several different sets of products... all in competition with collisional energy transfer. Each isomer corresponds to a ‘‘well” on the global potential energy surface and each isomerization and decomposition reaction constitutes a ‘‘reaction channel”. Since the reaction rate constant for each reaction channel depends on internal energy and collisional energy transfer is constantly affecting the energy (and angular momentum) distributions, such reaction systems can be modeled accurately only by using master equation simulation methods. The first steps in the oxidation of acetyl free radical (CH3C(O)) provide is a relatively simple example of a multi-well, multi-channel reaction system [28,134]. The acetyl radical is formed in the atmosphere and in combustion by abstraction of the aldehydic H-atom from acetaldehyde. Acetaldehyde is ubiquitous as a product of the atmospheric photo-oxidation of larger organic compounds and as a product of incomplete combustion in, for example, diesel engines. When acetyl radical combines with O2, it produces a vibrationally excited AcO2 free radical (the asterisk denotes vibrational excitation), which may undergo further reactions, or be thermalized by collisions. The over-all reaction has been shown to produce OH radical under some conditions[135,136] Many aspects of this reaction system have been studied experimentally [137–143], by quantum chemical calculations, [28,134,144] and by master equation simulations [28]. The potential energy surface for the first few steps in the oxidation of acetyl free radical are shown in Fig. 5. The initial reaction between acetyl radical and O2 is strongly exothermic and the acetylperoxyl free radical (A1) is initially highly excited. At low temperatures and high pressures (e.g. in the atmosphere), the initially excited A1 is quickly deactivated by collisional energy transfer. At lower pressures, collisional energy transfer is slower and A1 retains its excitation energy long enough to isomerize to A2. This ‘‘tail-biting” isomerization proceeds through a cyclic transition state; alkyl and alkoxyl free radicals undergo similar reactions. Subsequent isomerization to A20 , another conformer of A2, is followed by a concerted elimination of OH radical and its cyclic co-product C3 (an a-lactone). The relative yield of OH radical (and C3) depends on temperature and total pressure. Master equation simulations of this reaction system were carried out recently by Maranzana et al., [28] who also carried out

614

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617

Fig. 5. Potential energy diagram of the first steps in the oxidation of acetyl free radical. See Maranzana et al. [28] for details.

0

A2 ! HO2 þ Ketene

ð26Þ

The time-dependent concentrations of products from the present simulations are shown in Fig. 6, relative to the initial concentration of nascent acetyl peroxyl radical (A1). The initial energy distribution of A1 (the chemical activation distribution) extends from the energy threshold for decomposition of A1 up to considerably higher energies. The A1 molecules with the highest energies react very quickly by re-dissociating to acetyl radical (Ac) and O2. The A1 molecules with lower energies have longer lifetimes and can experience deactivating collisions. Even after a few collisions they are still energized, however, and continue to react mostly via the lower energy ‘‘tail-biting” pathway from A1 to A2 (see Fig. 5). After isomerization, the newly formed A2 retains its excitation energy and can react by returning to A1 via the same pathway by which it was formed, or it can react via the much lower energy pathway to form A20 . The A20 isomer has three reaction channels: return to A2, concerted elimination of OH radical to produce the a-lactone (C3), and bond fission to produce HO2 + Ketene. Fig. 6 shows that the concentration of A1 decreases monotonically to only a few percent of the initial concentration. Most has reacted, but on this time scale (and the somewhat longer time scale of Fig. 7) some has been collisionally stabilized (thermalized). Initially, the isomers A2 and A20 are not present, but they quickly Ac + O2

Relative Concentration

100

10-1 A1 OH + C3

10-2

Ac + O2

100

Relative Concentration

chemical structure calculations, which are in good agreement with previous work [134,144] Fig. 5 shows a subset of the system modeled by Maranzana et al., who presented complete details of the calculations. In addition to the reactions depicted in Fig. 5, the following reaction was also included:

10-1 OH + C3

10-2

A1

10-3

HO2 + CH2CO

10-4

A2' A2

10-5 400

500

600

700

800

900

1000

Temperature (K) Fig. 7. Final yields of species relative to the initial concentration of acetyl peroxyl radical (A1) vs. temperature at 1 bar of N2.

build up to maxima and then decay to final relative concentrations of 0.1%. The yields of OH, HO2, ketene and C3 asymptotically approach limiting values. In the end, the concentrations of OH radical and the acetyl radical are quite comparable at 700 K. Fig. 7 shows that the ‘‘final” relative yields of OH and ketene depend only weakly on temperature. The figures show that this relatively simple system exhibits complicated behavior. Experiments would be expected to behave in much the same way. The energy distributions of the ‘‘wells” (i.e. A1, A2, and A20 ) approach steady-state only after a large number of collisions. Even at steady-state, the energy distributions are non-Boltzmann: they would be Boltzmann only if no reactions were taking place. Before steady-state has been achieved, reaction rate constants cannot be defined, since the reaction flux coefficients (non-steady-state rate coefficients) vary with time. Even at steady-state, the reaction rate constants are non-thermal. What a mess! But these are the situations where master equation simulations excel, since they are designed to deal with them. They are invaluable tools for modeling real-world reaction systems.

A2'

10-3

10-4

A2 HO2 + CH2CO

0

5 10-10

1 10-9

1.5 10-9

2 10-9

2.5 10-9

Time (s) Fig. 6. Simulated acetyl +O2 reaction for 1 bar of N2 and 700 K. The concentration of species relative to the initial concentration of acetyl peroxyl radical (A1) vs. time.

3.5. Other systems In a study of the decomposition of CH3 radicals to CH2 + H and CH + H2[145], values of DEd at 2800 K of 150 cm1 were used in a fit to the data for both reaction channels with Ar as the bath gas. In a study involving CH3OH dissociation [146], using a very different method to solve the master equation than the use of

615

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617 Table A1 Parameters for multiwell calculations.

Table A3 Parameters for multiwell calculations.

CH4 Critical energy at 0 (K/kcal mol1) Vibrational frequencies (cm1) (J-rotor) Adiabatic moments of inertia (AMU A2) (K-rotor) Active external rotor (AMU A2) Symmetry; electronic degeneracy; optical isomers H–CH3 (Gorin transition state) Frequencies (cm1)

CH3Cl Critical energy at 0 (K/kcal mol1) Vibrational Frequencies (cm1)

103.3 3019 (3), 2917, 1534 (2), 1306 (3) 3.19 3.19 12; 1; 1

81.9 3042(2), 2966, 1455(2), 1355, 1015(2), 732 38.0 3.31 3; 1; 1

(J-rotor) Adiabatic moments of inertia (AMU Å2) (K-rotor) Active external rotor (AMU Å2) Symmetry; electronic degeneracy; optical isomers Cl—CH3 (Gorin transition state) Frequencies (cm1)

3184 (2), 3002, 1383 (2), 580

3184(2), 3002, 1383(2), 580

2

(J-rotor) Adiabatic moments of inertia (AMU A ) 300 K, 500 K, 1000 K, 1500 K, 2000 K

2

(J-rotor) Adiabatic moments of inertia (AMU Å ) 200 K, 250 K, 300 K, 1600 K, 2100 K

27.9, 24.9, 21.1, 19.0, 17.6 3.60 1.75

(K-rotor) Active external rotor (AMU A2) Moments of inertia active CH3 2-D rotor (AMU A2) Hindrance (to match k1 (cm3 mol1 s1) = 3.5  1010) 300 K, 500 K, 1000 K, 1500 K, 2000 K 78%, 79%, 80%, 82%, 83% Symmetry; electronic degeneracy; 3; 1; 1 optical isomers 3.33; 94.9 Collisions: (r (A2); e (K)) CH4 Ar 3.75; 98.3 He 2.55; 10.0 hDEi300 (cm1) 100(He), 150(Ar), 400(CH4), 500(C2H6) hDEdiT (cm1) = hDEdi300 (cm1) (T/300)

Multiwell, a value of DEd = 133(T/300)0.8 in Ar was used. (The values for the rate constant recommended in that study were confirmed experimentally in [147].) The bottom line seems to be that the lack of real understanding of collisional energy transfer requires, at a minimum, a balance between the resources and time spent computing detailed potential energy surfaces and efforts to come to grips with this lack of understanding. If energy transfer is to be treated as other than a parameter to explain or extrapolate data for unimolecular reactions (or their reverse), focus on these processes is required.

456.2, 435.9, 419.6, 283.0, 263.0 3.31 1.75

(K-rotor) Active external rotor (AMU Å2) Moments of inertia active CH32-D Rotor (AMU Å2) Hindrance (to match k1 from Parker et al. [119]) 200 K, 250 K, 300 K, 1600 K, 2100 K Symmetry; electronic degeneracy; optical isomers Collisions: (r (Å); e (K)) CH3Cl Ar He hDEi200 (cm1), hDEi250 (cm1), hDEi300 (cm1) hDEdi1600 (cm1), hDEdi2100 (cm1)

74%, 75%, 77%, 70%, 70% 3; 1; 1 4.07;373 3.75; 98.3 2.55; 10.0 400(He) 1500(Ar), 1700(Ar),

Table A4 Simulation parameters.

Norbornene Krypton Transition state

rLJ

eLJ

(Å)

(K)

5.5 3.655 –

330 178.9 –

Symmetry

Electronic degeneracy

Optical isomers

1 – 1

1 – 1

1 – 1

Table A5 Rotational constantsa.

Acknowledgments At Stanford University, this effort has been supported a by Grant CHE-0535555 from The National Science Foundation. At The University of Michigan, this work has been supported, in part, by The National Science Foundation (Atmospheric and Geospace Sciences) and the National Aeronautics and Space Administration (Upper Atmosphere Research Program).

K-rotor (1-D rotor) J-rotor (2-D rotor) a

Norbornene (cm1) [126]

Transition state (cm1) [127]

0.146

0.146

0.107

0.107

Centrifugal corrections were employed.

Table A2 Comparison of equilibrium constants (Keq/cm3 mol1) CH3 + H = CH4. T (K)

Keq (Ba)

k1-assn

k1-dissn

Keq–HiP

k0-assn

k0-dissn

Keq–LoP

Keq/HiP

Keq/LoP

Keq(MW)

Ba/MW

200 300 400 500 600 700 800 900 1000 1200 1500 2000 2500

1.31E+89 6.12E+50 4.25E+31 1.38E+20 3.03E+12 1.04E+07 8.25E+02 5.37E01 1.52E03 2.30E07 3.52E11 5.45E15 2.85E17

3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10 3.50E10

5.33E99 8.80E61 1.13E41 3.30E30 1.45E22 4.19E17 5.21E13 7.97E10 2.81E07 1.87E03 1.24E+01 8.22E+04 1.61E+07

6.57E+88 3.98E+50 3.10E+31 1.06E+20 2.41E+12 8.36E+06 6.72E+02 4.39E01 1.24E03 1.87E07 2.83E11 4.26E15 2.17E17

1.23E28 5.91E29 3.52E29 2.36E29 1.70E29 1.29E29 1.01E29 8.18E30 6.77E30 4.87E30 3.26E30 1.94E30 1.30E30

4.35E106 5.22E73 1.81E56 1.52E46 6.26E40 3.33E35 1.16E31 6.65E29 1.07E26 2.17E23 5.10E20 9.30E17 5.56E15

2.82E+77 1.13E+44 1.95E+27 1.56E+17 2.71E+10 3.87E+05 8.69E+01 1.23E01 6.35E04 2.25E07 6.39E11 2.09E14 2.34E16

2.00 1.54 1.37 1.30 1.26 1.24 1.23 1.22 1.22 1.23 1.25 1.28 1.32

4.66E+11 5.40E+06 2.18E+04 8.85E+02 1.12E+02 26.77 9.50 4.36 2.40 1.02 0.55 0.26 0.12

3.26E+88 4.03E+50 3.83E+31 1.36E+20 3.02E+12 1.01E+07 7.73E+02 4.85E01 1.33E03 1.91E07 2.79E11 4.26E15 2.31E17

4.02 1.52 1.11 1.01 1.01 1.03 1.07 1.11 1.14 1.21 1.26 1.28 1.24

Keq(Ba); k1-assn; k1-dissn; k0-assn; k0-dissn are values recommended by Baulch et al.[107]. Keq–HiP; Keq–LoP are values of the equilibrium constant obtained from the ratio of k’s. Keq/HiP; and Keq/LoP are ratios of Keq(Ba) to Keq–HiP and Keq–LoP respectively. Keq(MW) is the value computed with the Thermo code in Multiwell and Ba/MW is the ratio of Keq(Ba) to Keq(MW). Bold values are in the recommended temperature range. k0-dissn had two different recommendations depending on the temperature range. The values in italics are from the lower temperature range.

616

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617 Table A6 Vibrational frequencies (cm1). [4] Norbornene [126]

Transition state [127]

3090.6 2997.2 2980.8 2975.5 2932.8 2926.6 1574.6 1477.7 1457.7 1299.5 1284.1 1167.7 1126.6 1093.2 1021.3 964.2 938.6 906.7 873.5 809.9 769.4 709.8 471.7 381.0 3063.2 2991.0 2959.5 2916.2 1455.6 1339.9 1285.8 1270.9 1254.2 1206.2 1177.2 1115.2 1035.1 950.7 915.0 898.9 833.1 794.1 664.2 494.5 258.0

3091 3105 3105 3102 3026 2960 1580 1500 1500 1500 1300 1300 1126 1202 1123 1060 1032 997 900 940 320 435 590 395 3075 2988 3043 2886 1500 900 1350 1350 1206 1294 1226 1138 950 1006 988 833 664 220 129 325

[5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23]

[24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40]

Appendix A A.1. Norbornene: shock-initiated isomerization simulations All parameters are discussed in detail elsewhere [127]. Notes on simulations: (1) Energy transfer model: exponential-down with a(e)/ cm1 = 70 + 0.003e. (2) Reaction critical energy: E0 = 182.8 kJ mol1. (3) Simulations were initiated using a thermal energy distribution with Tvib = 300 K. (4) Ttrans = shock translational temperature.

[41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54] [55] [56]

Tables A4, A5, A6 References [1] W.H. Green Reaction Mechanism Generator. . [2] H.-H. Carstensen, A.M. Dean, J. Phys. Chem. A 113 (2) (2009) 367–380. [3] G.P. Smith, D.M. Golden, M. Frenklach, N.W. Moriarty, B. Eiteneer, M. Goldenberg, C.T. Bowman, R.K. Hanson, S. Song, W.C. Gardiner Jr., V.V.

[57] [58] [59] [60] [61] [62]

Lissianski, Z. Qin, GRI-Mech 3.0. . M. Frenklach, H. Wang, M.J. Rabinowitz, Prog. Energy Combust. Sci. 18 (1) (1992) 47–73. M. Frenklach, A. Packard, R. Feeley, in: W.C. Robert (Ed.), Comprehensive Chemical Kinetics, vol. 42, Elsevier, 2007. T. Russi, A. Packard, R. Feeley, M. Frenklach, J. Phys. Chem. A 112 (12) (2008) 2579–2588. D.C. Tardy, C.W. Larson, B.S. Rabinovitch, Can. J. Chem. 46 (1968) 341–343. D.C. Tardy, B.S. Rabinovitch, J. Chem. Phys. 45 (10) (1966) 3720–3730. T. Baer, W.L. Hase, Unimolecular Reaction Dynamics. Theory and Experiments, Oxford University Press, New York, 1996. W. Forst, Unimolecular Reactions. A Concise Introduction, Cambridge University Press, Cambridge, 2003. R.G. Gilbert, S.C. Smith, Theory of Unimolecular and Recombination Reactions, Blackwell Scientific, Oxford, 1990. K.A. Holbrook, M.J. Pilling, S.H. Robertson, Unimolecular Reactions, Wiley, Chichester, 1996. J.R. Barker, Chem. Phys. 77 (2) (1983) 301–318. J.R. Barker, D.M. Golden, Chem. Rev. 103 (12) (2003) 4577–4591. M.J. Pilling, S.H. Robertson, Annu. Rev. Phys. Chem. 54 (2003) 245–275. J.A. Miller, S.J. Klippenstein, J. Phys. Chem. A 110 (2006) 10528–10544. S.J. Klippenstein, Y. Georgievskii, L.B. Harding, Phys. Chem. Chem. Phys. 8 (2006) 1133–1147. A. Fernandez-Ramos, J.A. Miller, S.J. Klippenstein, D.G. Truhlar, Chem. Rev. 106 (11) (2006) 4518–4584. S.J. Klippenstein, in: N.J.B. Green (Ed.), Unimolecular Kinetics. Part I. The Reaction Step, Vol. 39, Elsevier, Amsterdam, 2003. L.T. Nguyen, J.F. Stanton, J.R. Barker, Chem. Phys. Lett., accepted for publication. J.R. Barker, L.M. Yoder, K.D. King, J. Phys. Chem. A 105 (5) (2001) 796–809. J.R. Barker, N.F. Ortiz, Int. J. Chem. Kinet. 33 (4) (2001) 246–261. J.R. Barker, N.F. Ortiz, J.M. Preses, L.L. Lohr, A. Maranzana, P.J. Stimac, T.L. Nguyen, Ann Arbor, Michigan, USA, 2009. . D.T. Gillespie, J. Phys. Chem. 81 (25) (1977) 2340–2361. D.T. Gillespie, J. Comp. Phys. 22 (4) (1976) 403–434. J.R. Barker, J. Chem. Phys. 72 (6) (1980) 3686–3694. J.R. Barker, Int. J. Chem. Kinet. 33 (4) (2001) 232–245. A. Maranzana, J.R. Barker, G. Tonachini, Phys. Chem. Chem. Phys. 9 (31) (2007) 4129–4141. J.R. Barker, Int. J. Chem. Kinet. 41 (2009) 748–763. F.A. Lindemann, Trans. Faraday Soc. 17 (1922) 598–599. J.A. Miller, S.J. Klippenstein, S.H. Robertson, M.J. Pilling, N.J.B. Green, Phys. Chem. Chem. Phys. 11 (2009) 1128–1137. J.R. Barker, R.E. Weston Jr., J. Phys. Chem. A, accepted for publication. T. Beyer, D.F. Swinehart, Comm. Assoc. Comput. Machines 16 (6) (1973) 379. S.E. Stein, B.S. Rabinovitch, J. Chem. Phys. 58 (6) (1973) 2438–2445. T.L. Nguyen, J.R. Barker, J. Phys. Chem. A 114 (2010) 3718–3730. S. Nordholm, Chem. Phys. 137 (1989) 109–120. D.M. Leitner, P.G. Wolynes, J. Phys. Chem. A 101 (4) (1997) 541–548. D.M. Leitner, B. Levine, J. Quenneville, T.J. Martinez, P.G. Wolynes, J. Phys. Chem. A 107 (49) (2003) 10706–10716. R.E. Weston Jr., J.R. Barker, J. Phys. Chem. A 110 (2006) 7888–7897. J.A. Miller, S.J. Klippenstein, C. Raffy, J. Phys. Chem. A 106 (19) (2002) 4904– 4913. L. Vereecken, G. Huyberechts, J. Peeters, J. Chem. Phys. 106 (16) (1997) 6564– 6573. R.A. Marcus, J. Chem. Phys. 43 (8) (1965) 2658–2661. E.V. Waage, B.S. Rabinovitch, Chem. Rev. 70 (1970) 377–387. S.C. Smith, R.G. Gilbert, Int. J. Chem. Kinet. 20 (1988) 307–329. E.W. Montroll, K.E. Shuler, J. Chem. Phys. 26 (1957) 454. H.O. Pritchard, The Quantum Theory of Unimolecular Reactions, Cambridge University Press, Cambridge, 1984. F.P. Buff, D.J. Wilson, J. Chem. Phys. 32 (3) (1960) 677–685. M. Hoare, J. Chem. Phys. 38 (7) (1963) 1630. W. Forst, J. Phys. Chem. 76 (3) (1972) 342–348. S.J. Klippenstein, J.A. Miller, J. Phys. Chem. A 106 (40) (2002) 9267–9277. D.T. Gillespie, J. Comp. Phys. 28 (3) (1978) 395–407. P.K. Venkatesh, Physica A (Amsterdam) 289 (18) (2001) 359–376. M. Basire, P. Parneix, F. Calvo, J. Chem. Phys. 129 (2008) 081101. F. Wang, D.P. Landau, Phys. Rev. Lett. 86 (2001) 2050–2053. E.E. Aubanel, D.M. Wardlaw, L. Zhu, W.L. Hase, Int. Rev. Phys. Chem. 10 (3) (1991) 249–286. M.A. Hanning-Lee, N.J.B. Green, M.J. Pilling, S.H. Robertson, J. Phys. Chem. 97 (4) (1993) 860–870. S.H. Robertson, M.J. Pilling, D.L. Baulch, N.J.B. Green, J. Phys. Chem. 99 (36) (1995) 13452–13460. D.W. Wardlaw, R.A. Marcus, J. Chem. Phys. 110 (3) (1984) 230–234. S.J. Klippenstein, L.R. Khundkar, A.H. Zewail, R.A. Marcus, J. Chem. Phys. 89 (8) (1988) 4761–4770. S.J. Klippenstein, Chem. Phys. Lett. 170 (1) (1990) 71–77. S.J. Klippenstein, J. Chem. Phys. 96 (1) (1992) 367–371. S.J. Klippenstein, A.F. Wagner, S.H. Robertson, R. Dunbar, D.M. Wardlaw, 1.0 ed., 1999. .

D.M. Golden, J.R. Barker / Combustion and Flame 158 (2011) 602–617 [63] E.E. Greenwald, S.W. North, Y. Georgievskii, S.J. Klippenstein, J. Phys. Chem. A 111 (25) (2007) 5582–5592. [64] D.M. Golden, Int. J. Chem. Kinet. 37 (10) (2005) 625–632. [65] E. Gorin, Acta Physicochim., URSS 9 (1938) 681. [66] S.W. Benson, Thermochemical Kinetics, Wiley, New York, 1976. [67] G.P. Smith, D.M. Golden, Int. J. Chem. Kinet. 10 (5) (1978) 489–501. [68] J.E. Lennard-Jones, Proc. R. Soc. London Ser. A 106 (1924) 463–477. [69] P.M. Morse, Phys. Rev. 34 (1929) 57–64. [70] Y.P. Varshni, Rev. Mod. Phys. 29 (1957) 664–682. [71] D.C. Tardy, B.S. Rabinovitch, Chem. Rev. 77 (1977) 369–408. [72] J. Troe, Ber. Bunsen-Ges. Phys. Chem. 77 (1973) 665–674. [73] J. Troe, J. Chem. Phys. 66 (1977) 4745–4757. [74] V. Bernshtein, I. Oref, J. Phys. Chem. 97 (1993) 12811–12818. [75] I. Oref, in: J.R. Barker (Ed.), Vibrational Energy Transfer Involving Large and Small Molecules, vol. 2B, JAI Press Inc., Greenwich, CT, 1995. [76] K. Frerichs, M. Hollerbach, T. Lenzer, K. Luther, J. Phys. Chem. A 110 (2006) 3179–3185. [77] U. Hold, T. Lenzer, K. Luther, A.C. Symonds, J. Chem. Phys. 119 (21) (2003) 11192–11211. [78] T. Lenzer, K. Luther, K. Reihs, A.C. Symonds, J. Chem. Phys. 112 (9) (2000) 4090–4110. [79] U. Hold, T. Lenzer, K. Luther, K. Reihs, A.C. Symonds, J. Chem. Phys. 112 (9) (2000) 4076–4089. [80] K. Luther, K. Reihs, Ber. Bunsen-Ges. Phys. Chem. 92 (1988) 442. [81] D.L. Bunker, S.A. Jayich, Chem. Phys. 13 (2) (1976) 129–134. [82] A.J. Stace, J.N. Murrell, J. Chem. Phys. 68 (7) (1978) 3028–3039. [83] N.J. Brown, J.A. Miller, J. Chem. Phys. 80 (1984) 5568. [84] H. Hippler, H.W. Schranz, J. Troe, J. Phys. Chem. 90 (23) (1986) 6158–6167. [85] H.W. Schranz, J. Troe, J. Phys. Chem. 90 (23) (1986) 6168–6175. [86] A.R. Whyte, K.F. Lim, R.G. Gilbert, W.l. Hase, Chem. Phys. Lett. 152 (1988) 377–381. [87] G. Lendvay, G.C. Schatz, J. Phys. Chem. 94 (1990) 8864–8866. [88] K.F. Lim, R.G. Gilbert, J. Phys. Chem. 94 (1990) 77–84. [89] W.L. Hase, C.L. Darling, L. Zhu, J. Chem. Phys. 96 (11) (1992) 8295–8306. [90] G. Lendvay, G.C. Schatz, in: J.R. Barker (Ed.), Vibrational Energy Transfer Involving Large and Small Molecules, vol. 2B, JAI Press Inc., Greenwich, CT, 1995. [91] T. Lenzer, K. Luther, J. Troe, R.G. Gilbert, K.F. Lim, J. Chem. Phys. 103 (2) (1995) 626–641. [92] T. Lenzer, K. Luther, J. Chem. Phys. 104 (9) (1996) 3391–3394. [93] M.J.T. Jordan, D.C. Clary, J. Chem. Phys. 106 (13) (1997) 5439–5453. [94] A.S. Mullin, G.C. Schatz, in: A. Mullin, G.C. Schatz (Eds.), Highly Excited States: Relaxation, Reaction, and Structure, vol. 678, American Chemical Society, Washington, DC, 1997. [95] L.M. Yoder, J.R. Barker, J. Phys. Chem. A 104 (45) (2000) 10184–10193. [96] D. Troya, J. Phys. Chem. A 109 (26) (2005) 5814–5824. [97] J.C. Lorquet, B. Leyh-Nihantt, J. Phys. Chem. 92 (1988) 4784–4787. [98] F. Remacle, D. Dehareng, J.C. Lorquet, J. Phys. Chem. 92 (1988) 4784–4787. [99] J.N. Harvey, M. Aschi, Phys. Chem. Chem. Phys. 1 (1999) 5555–5563. [100] T.J. Bevilacqua, R.B. Weisman, J. Chem. Phys. 98 (8) (1993) 6316. [101] D.R. McDowell, F. Wu, R.B. Weisman, J. Phys. Chem. 101 (29) (1997) 5218. [102] D.R. McDowell, F. Wu, R.B. Weisman, J. Chem. Phys. 108 (22) (1998) 9404. [103] F. Wu, R.B. Weisman, J. Chem. Phys. 110 (11) (1999) 5047. [104] F. Wu, R.B. Weisman, J. Chem. Phys. 112 (23) (2000) 10173–10178. [105] R.G. Gilbert, K. Luther, J. Troe, Ber. Bunsen-Ges. Phys. Chem. 87 (1983) 169. [106] J. Troe, J. Phys. Chem. 83 (1) (1979) 114–126. [107] D.L. Baulch, C.T. Bowman, C.J. Cobos, R.A. Cox, T. Just, J.A. Kerr, M.J. Pilling, D. Stocker, J. Troe, W. Tsang, R.W. Walker, J. Warnatz, J. Phys. Chem. Ref. Data 34 (3) (2005) 757–1397. [108] R. Atkinson, D.L. Baulch, R.A. Cox, J.N. Crowley, R.F. Hampson, R.G. Hynes, M.E. Jenkin, M.J. Rossi, J. Troe, Atmos. Chem. Phys. 7 (2007) 981–1191. [109] S.P. Sander, J. Abbatt, J.R. Barker, J.B. Burkholder, R.R. Friedl, D.M. Golden, R.E. Huie, C.E. Kolb, M.J. Kurylo, G.K. Moortgat, V.L. Orkin, P.H. Wine, in: Jet

[110] [111] [112] [113] [114] [115] [116] [117] [118] [119] [120] [121] [122] [123] [124] [125] [126] [127] [128] [129] [130] [131] [132] [133] [134] [135] [136] [137] [138]

[139] [140] [141] [142] [143] [144] [145] [146] [147]

617

Propulsion Laboratory, Pasadena, JPL 09-31, 2009. . D.M. Golden, Chem. Soc. Rev. 37 (2008) 717–731. C.J. Cobos, J. Troe, Z. Phys. Chem. 167 (1990) 129–149. L.B. Harding, Y. Georgievskii, S.J. Klippenstein, J. Phys. Chem. A 109 (21) (2005) 4646–4656. D.M. Golden, Int. J. Chem. Kinet. 40 (6) (2008) 310–319. K.A. Boering, J.I. Brauman, J. Chem. Phys. 97 (8) (1992) 5439–5450. D.M. Golden, J. Phys. Chem. A 111 (29) (2007) 6772–6780. R. Walsh, D.M. Golden, J. Phys. Chem. A 112 (2008) 3891–3897. D.M. Golden, J. Phys. Chem. A 110 (9) (2006) 2940–2943. D.M. Golden, Int. J. Chem. Kinet. 35 (2003) 206–211. J.K. Parker, W.A. Payne, R.J. Cody, F.L. Nesbitt, L.J. Stief, S.J. Klippenstein, L.B. Harding, J. Phys. Chem. A 111 (6) (2007) 1015. D.M. Golden, Int. J. Chem. Kinet. 41 (4) (2009) 245–254. J.T. Yardley, Introduction to Molecular Energy Transfer, Academic Press, New York, 1980. J.H. Kiefer, in: A. Lifshitz (Ed.), Shock Waves in Chemistry, Marcel Dekker, New York, 1981. J.H. Kiefer, J.N. Shah, J. Phys. Chem. 91 (1987) 3024. J.H. Kiefer, L.L. Buzyna, A. Dib, S. Sundaram, J. Chem. Phys. 113 (1) (2000) 48– 58. J. Shi, J.R. Barker, Int. J. Chem. Kinet. 22 (1990) 187–206. J.H. Kiefer, S.S. Kumaran, S. Sundaram, J. Chem. Phys. 99 (5) (1993) 3531– 3541. J.R. Barker, K.D. King, J. Chem. Phys. 103 (1995) 4953–4966. J.H. Kiefer, C. Katapodis, S. Santhanam, N.K. Srinivasan, R.S. Tranter, J. Phys. Chem. A 108 (2004) 2443–2450. J.R. Barker, P.J. Stimac, K.D. King, D.M. Leitner, J. Phys. Chem. A 110 (2006) 2944–2954. M.J. Davis, J.H. Kiefer, J. Chem. Phys. 116 (18) (2002) 7814–7827. W.C. Herndon, W.B. Cooper, M.J. Chambers, J. Phys. Chem. 68 (1964) 2016. B.C. Roquitte, J. Phys. Chem. 69 (1965) 1351. R. Walsh, J.M. Wells, J. Chem. Soc. Perk T 2 (1) (1976) 52–55. H. Hou, A. Li, H.Y. Hu, Y. Li, H. Li, B.S. Wang, J. Chem. Phys. 122 (22) (2005) 224304. G.S. Tyndall, T.A. Staffelbach, J.J. Orlando, J.G. Calvert, Int. J. Chem. Kinet. 27 (10) (1995) 1009–1020. J.V. Michael, D.G. Keil, R.B. Klemm, J. Chem. Phys. 83 (4) (1985) 1630–1636. N.I. Butkovskaya, A. Kukui, G. Le Bras, J. Phys. Chem. A 108 (7) (2004) 1160– 1168. R.K. Talukdar, M.E. Davis, L. Zhu, A.R. Ravishankara, J.B. Burkholder, in: Determination of the OH Radical Yield in the CH3CO + O2 Reaction, 19th International Symposium on Gas Kinetics, Orleans, France, July 22–27, 2006, Orleans, France, 2006. M.A. Blitz, D.E. Heard, M.J. Pilling, Chem. Phys. Lett. 365 (5–6) (2002) 374– 379. B. D’Anna, V. Bakken, J.A. Beukes, C.J. Nielsen, K. Brudnik, J.T. Jodkowski, Phys. Chem. Chem. Phys. 5 (9) (2003) 1790–1805. P. Devolder, S. Dusanter, B. Lemoine, C. Fittschen, Chem. Phys. Lett. 417 (1–3) (2006) 154–158. G. Kovács, J. Zádor, E. Farkas, R. Nádasdi, I. Szilágyi, S. Dóbé, T. Bérces, F. Márta, G. Lendvay, Phys. Chem. Chem. Phys. 9 (2007) 4142–4154. S.A. Carr, M.T. Baeza-Romero, M.A. Blitz, M.J. Pilling, D.E. Heard, P.W. Seakins, Chem. Phys. Lett. 445 (4–6) (2007) 108–112. J.W. Lee, C.J. Chen, J.W. Bozzelli, J. Phys. Chem. A 106 (31) (2002) 7155–7170. V. Vasudevan, R.K. Hanson, D.M. Golden, C.T. Bowman, D.F. Davidson, J. Phys. Chem. A 111 (19) (2007) 4062–4072. A.W. Jasper, J.A. Miller, J. Phys. Chem. A 113 (2009) 5612–5619. V. Vasudevan, R.D. Cook, R.K. Hanson, C.T. Bowman, D.M. Golden, Int. J. Chem. Kinet. 40 (8) (2008) 488–495.