When Rate Constants Are Not Enough.

Subscriber access provided by SUNY DOWNSTATE

Article

When Rate Constants Are Not Enough John R. Barker, Michael Frenklach, and David M Golden J. Phys. Chem. A, Just Accepted Manuscript • Publication Date (Web): 13 Apr 2015 Downloaded from http://pubs.acs.org on April 13, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry A is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

When Rate Constants Are Not Enough John R. Barker,a* Michael Frenklach,b* and David M. Goldenc* a Department

of Atmospheric, Oceanic and Space Sciences, University of Michigan, Ann Arbor, MI 48109-2143 b Department of Mechanical Engineering, University of California, Berkeley, CA 94720-1740 c Department of Mechanical Engineering, Stanford University, Stanford, CA 94305-3032

Abstract Real-world chemical systems consisting of multiple isomers and multiple reaction channels often react significantly prior to attaining a steady state energy distribution (SED). Detailed elementary reaction models, which implicitly require SED conditions, may be invalid when non-steady state energy distributions (NSED) exist. NSED conditions may result in reaction rates and product yields that are different than those expected for SED conditions, although this problem is to some extent reduced by using phenomenological models and rate constants. The present study defines pragmatic diagnostics useful for identifying NSED conditions in stochastic master equation simulations. A representative example is presented for each of four classes of common combustion species: RO2 radicals, aliphatic hydrocarbons, alkyl radicals, and polyaromatic radicals. An example selected from the seminal work of Tsang et al. demonstrates that stochastic simulations and eigenvalue methods for solving the master equation predict the same NSED effects. NSED effects are common under relatively moderate combustion conditions and accurate simulations require a master equation analysis.

Introduction Chemical reaction models that include all of the relevant reaction steps and account for their temperature and pressure dependences are exceptionally useful for modeling both laboratory experiments and practical chemical systems. With this goal in mind, researchers assemble “detailed reaction models”, comprised of large sets of elementary reactions, for use in modeling combustion and atmospheric chemistry. One of the principal challenges in this effort is obtaining accurate elementary-reaction rate coefficients.

*

Correspondence: [email protected], [email protected], [email protected]

ACS Paragon Plus Environment

1


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 29

Typically, only a few of these rate coefficients have ever been measured in experiments, and the experiments have typically only explored a small fraction of the temperature, pressure, and composition ranges encountered in realistic combustion systems. Hence, the community has turned to theoretical methods, which range from semi-empirical estimation to fully ab initio treatments. One motivation for the present paper is to call renewed attention to the fact that representation of a reacting chemical system by elementary reactions is only applicable under certain conditions, which may not be attained during combustion processes. A second motivation is to describe ways that the suitability of the conditions can be monitored in stochastic simulations of reacting systems. We start by recalling that the “elementary reaction” is a useful abstraction, associated with a population distributed over a huge number of quantum states. The elementary rate constant is defined as an appropriate time independent average over the distributed population. Such averages typically involve distributions over translational, rotational, vibrational, and electronic states. When the distributions vary with time, the averages also vary, thus invalidating the primary condition of the “elementary reaction rate constant” and thereby making its use in detailed models inappropriate. Phenomenological models can overcome this limitation, at least to some extent,1-6 but assessment of their accuracy remains an important task. Energy Distributions and Diagnostics It has been long known that when molecules are placed in a heat bath, their translational energy distribution is collisionally equilibrated very rapidly, but equilibration of the internal degrees of freedom takes longer. When fully equilibrated there exists a steady state energy distribution (SED). If these molecules can undergo unimolecular reactions there will be some conditions of pressure and temperature in which reaction proceeds prior to complete equilibration (i.e., before an SED is established). During this time period, (i.e., in which there exists a nonsteady state energy distribution, or NSED), reactions cannot be described as having time independent rate constants. As specifically pointed out by Tsang, Bedanov, and Zachariah in their landmark papers,7, 8 this behavior is not incorporated into programs such as CHEMKIN9 that are used to simulate chemical processes.

The interplay of collisional energy transfer and chemical reaction can be described by a coupled set of differential equations, called a Master Equation (ME),10 which describes the energy distribution and concentration of a reacting chemical species. When the reacting species is highly diluted in an inert bath gas, energy transfer collisions between two excited reacting molecules can be neglected and the ME is approximately linear. Due to its linear nature, the formal solution of a ME may be expressed in a matrix form, which is commonly solved by an eigenvector-eigenvalue decomposition.11 The ME can also be solved by stochastic simulation,12, 13 which is


2

Page 3 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


discussed below and in the Appendix. The papers of Tsang and co-workers7, 8 offer an excellent description of a matrix solution. (An open source code5 for solving the master equation is available.) The solution to the matrix equation is expressed by the eigenvalues, from which time dependent concentration profiles can be obtained. The eigenvalues tend to fall into two groups termed4 the “internal relaxation eigenvalues" (IEREs) and the “chemically significant eigenvalues" (CSEs), an SED is considered to exist and time independent rate constants can be defined when the magnitudes of the CSEs are much smaller than those of the IEREs. In general, all of the eigenvalues are affected by both energy transfer and chemical transformation, but the IEREs mostly describe the former and the CSEs mostly describe the latter. When the magnitudes of the IEREs and CSEs are sufficiently separated, the energy transfer may be regarded as fast enough so that the chemical rate can be expressed using a time independent phenomenological rate constant.2 Clearly, the existence of an SED is an idealization, which can only be achieved when the IEREs are infinitely faster than the CSEs (i.e. at the high pressure limit), or at infinite time, when the system is fully equilibrated. When this condition does not abide, concentration profiles can still be obtained, but an NSED exists and only phenomenological rate constants can be defined. This shows that elementary rate constants are idealizations, as well, but the system may be close enough to SED conditions so that elementary rate constants are wellapproximated by time independent phenomenological ones, in which the time dependence has been subsumed. A pragmatic diagnostic is needed in order to ascertain whether the energy distribution is sufficiently steady so that an SED may be pragmatically assumed to exist. In the matrix solution of the master equation, the diagnostic for this condition is to observe the separation or lack thereof between the IEREs and the CSEs; the pragmatic criterion is typically a factor of ≳10.5 NSED conditions can cause not only ill-defined elementary rate constants, but chemical species may react away before steady state is achieved. However, as pointed out by Tsang and coworkers7, 8, even when rate constants are apparently "close enough" in value to those at steady state, product branching ratios may still be affected. As mentioned earlier, phenomenological models with phenomenological rate constants can overcome this problem to some extent.1-6 In isolated isomerization reactions, which have been studied extensively,3, 14 chemical equilibration and vibrational relaxation occur simultaneously. At high energies (above the isomerization energy barrier), the population distributions achieve microcanonical equilibrium very rapidly, due to the large forward and reverse microcanonical rate constants (k(E)s). At moderate pressures, this process is far more rapid than collisional vibrational energy transfer, which transports population up and down the energy ladders. When NSEDs exist in one well, they also exist in the other well, due to the microcanonical equilibrium established between the wells at each energy. As a result, vibrational relaxation occurs simultaneously in both wells. An eigenvalue solution of the master equation for this


3


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 29

process produces two CSEs and a set of IEREs describing evolution of the energy distributions due to collisional energy transfer. One of the CSEs is identically zero and the other governs relaxation of the entire system to equilibrium. Except at the high pressure limit, SEDs are only achieved in this process after chemical equilibration is attained, but useful phenomenological rate constants connecting the wells can still be defined.1, 3 The accuracy of the resulting phenomenological kinetics model depends on the separation of the IEREs and CSEs. For isolated isomerization reactions, the phenomenological model serves well, despite the existence of NSED conditions. The situation is not so clear, however, when irreversible reaction paths (e.g. fragmentation reactions) are added to the mix. In this case, the microcanonical equilibrium between the high energy levels in the two wells is disturbed to a greater or lesser degree, which affects the isomerization rate. Furthermore, and perhaps more importantly, the NSED conditions that are present in the two wells during equilibration may have important effects on the rates of the irreversible reactions. A phenomenological model can be readily formulated for this system,1, 4, 5, 15, 16 but its accuracy must be assessed. Some example chemical systems in this category are discussed below. In this paper we offer a quantitative diagnostic for use in stochastic simulations, in which the IEREs and CSEs are never known. In stochastic simulations,10,17 eigenvalues are not obtained. Instead, time dependent energy distributions and concentration profiles are generated by the stochastic simulation algorithm (SSA) defined by Gillespie.13, 18 Gillespie also proved rigorously that the SSA produces the same solutions as a set of coupled ordinary differential equations describing the same physics. Using the SSA to solve the ME raises the question as to how to determine whether the energy distribution is sufficiently steady to be considered an SED. One clear manifestation of NSED conditions is the existence of an incubation period in the time dependent concentration profile of a reactant or product in a unimolecular reaction system. However, other quantities routinely available in master equation simulations are more sensitive. The average active (or vibrational) energy Evib ( t ) has long been used for this purpose.19 As long as Evib ( t ) is varying with time, NSEDs are present, whether or not an incubation period is observed. However, Evib ( t ) is only the first moment of the energy distribution. Quack, in his study of unimolecular isomerization,3 examined the entire energy distribution, as did Tsang and coworkers.7, 8 Although examining the entire energy distribution is fully informative, it is not practical for most simulations. A more practical approach is to use Evib ( t ) and the reaction flux coefficients (see Appendix), which are easily computed by master equation codes. Reaction flux coefficients may be regarded as averages of the energy distribution, weighted by the microcanonical rate constants, ki(E). Thus they are


4

Page 5 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


sensitive to the energy distribution at the energies where reaction takes place. Since each well in a multi-well reaction system may be involved in several reactions, the reaction flux coefficients provide considerable information and are sensitive indicators, as shown in the examples presented below. However, since some of the reactions may have similar energy thresholds, the information may be somewhat redundant. In those cases, we have found it practical to use the sum of the reaction flux coefficients for all of the reactions originating in a given well: Nj

R j ( t ) = ∑ ri, j ( t )

(1)

i=1

where Rj(t) is the desired time dependent sum, Nj is the number of reaction channels originating from the jth well and rj,i(t) is the reaction flux coefficient for the ith channel. The MultiWell master equation code20 reports rj,i(t) for all reactions in output file “multiwell.rate” and Rj(t) is easily calculated from the output. NSEDs can originate in chemical activation reactions, in exothermic fragmentation reactions, as a result of photo reactions, in shock tubes, etc. The stochastic master equation simulation method is very well suited to simulating such reaction systems. For present purposes, we have mostly confined our attention to shock-initiation, but the results of our study apply widely to other ways of generating NSEDs.

Methods All of the simulations reported here were carried out using the MultiWell Program Suite.20 Parameters used in the simulations were obtained from previously published work, as cited in later sections of this paper. The multiwell ME code in this suite uses the same input information as the matrix based solutions and uses Gillespie's SSA. 13, 17 The main result of a calculation (output file "multiwell.out") is the set of species concentrations as a function of the time (and average number of collisions). The time dependent average vibrational (or active) energy in each well is reported in the same output file. Reaction flux coefficients are reported in output file "multiwell.rate". The parameters needed for the calculations are discussed in the original publications (see the individual examples, below) and MultiWell input data are included in the Supplementary Information. A discussion of some details of master equations and the MultiWell code is found in the Appendix.


5


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 29

Example Systems Ethyl + O2: A single-well system Chemical Activation An example of a well-studied system is the reaction of ethyl radicals with molecular oxygen. A potential energy diagram may be found in Figure 2 of Miller and Klippenstein.21 Using MultiWell we have modeled this system as a chemical activation initiated by the interaction of ethyl radicals with oxygen molecules at various temperatures and pressures of helium. Using input parameters that matched as nearly as possible those used by Miller and Klippenstein we were able to reproduce the results in Figure 3 in their paper with good fidelity [MultiWell input data are tabulated in the Supplementary Information]. The processes that follow the initial chemical activation process are: 1. 2. 3. 4.

C2H5O2** C2H5 + O2 C2H5O2** C2H4 + HO2 C2H5O2** + M C2H5O2 + M C2H5O2 + M C2H5O2** + M

Process 4 leads to thermal decomposition of stabilized C2H5O2. (Since process 4 does not produce C2H5O2 with the same energy distribution as the chemical activation process, the product of 4 might be denoted as C2H5O2* to differentiate it from C2H5O2**. This complicates some algebra, but is in fact easily handled in the MultiWell calculation.) Modeling this system by using MultiWell can be accomplished by treating it as a chemical activation process at a fixed temperature and pressure, in which the C2H5 and O2 come together to form an energized C2H5O2. The fractions of C2H5O2, C2H5 +O2 and C2H4+HO2 can be followed as functions of time, along with the reaction flux coefficients and the average vibrational energy in C2H5O2. Under these conditions the path to C2H4 and HO2 dominates. The overall rate constant for the loss of C2H5 + O2 as well as the rate constants for the individual channels can be obtained. An example of this type of calculation is shown in Figure 1.


6

Page 7 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 1: Surviving C2H5O2 and its average vibrational energy (red line). The green and blue lines show the fraction of C2H5O2 under NSED and SED conditions, respectively. The broken line shows the least squares fit of the SED portion. Since, as can be seen, the energy has not relaxed to a SED until 2 microseconds: the decay was fitted to an exponential function starting at that time. The total rate constant during the SED portion is 4.61 × 103 s-1. The linearity of the plot shows, however, that the total rate constant is essentially constant after the first few nanoseconds, even though NSED conditions persist until ~2 microseconds. Thus one might be led to believe mistakenly that NSED conditions have no effect. Figure 2 shows the branching ratio of product channels computed in the same simulations. The fact that the branching ratio depends on time reflects the existence of NSED conditions during the approach to SED. Although the total rate constant for C2H5O2 decomposition is not sensitive to the NSED in this case, the product branching ratios are sensitive to it. (Note that the stochastic noise could be reduced by running more stochastic trials, but that is not necessary for present purposes.)


7


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 29

Figure 2: Product channel branching ratio and average vibrational energy as functions of time for the same calculation illustrated in Figure 1. Shock Initiated Another way that this system can be examined is by simulating “shock” heating. The translational temperature is set at 800 K and the initial internal temperature is set at 300 K. Figure 3 shows the vibrational energy rising to the steady state value and the time delay (the incubation time) for the decomposition. If a rate constant is extracted after the 2 microsecond time period, a value of 4.53 × 103 s-1 is obtained, matching within statistical errors the value obtained from simple thermal decomposition. It is of interest to note that this value combined with the value for the association path computed in the chemical activation simulations yields the equilibrium constant in very good agreement with the value directly computed using the “thermo” code in the MultiWell suite. [Input data are tabulated in the Supplementary Information.] Yet another diagnostic for the approach to SED is the time dependence of the reaction flux coefficients, which are routinely available from MultiWell. An example is shown in Figure 4.

Figure 3: Shock heating from 300 K to 800 K. The rate constant for thermal decomposition of the C2H5O2 was obtained from the exponential fit to the blue line, starting at 2 microseconds. The green curve shows the fraction of C2H5O2 from t = 0 to 2 microseconds. The red curve shows the time dependence of the average vibrational energy.


8

Page 9 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: Reaction flux coefficients as a function of time for the same type of calculation as Figure 3. The lesson from this example is that even for this single well system, where the total rate constant reaches essentially SED values well before steady state is actually achieved, the product branching ratios are affected significantly by NSEDs. cis-Butene-2 Isomerization and Decomposition: A two well system This example, which was adapted from the work of Bedanov, Tsang, and Zachariah,8 was chosen because of its interesting behavior and because it provides a good test of NSED diagnostics in stochastic simulations. The system consists of cis- and trans-butene-2 (i.e. the cis- and trans- wells with index numbers #1 and #2, respectively, in this model), the fragmentation products (1,3-butadiene + H2, with index #3), and transition states for cis/trans isomerization and for fragmentation of cis-butene-2, and, as shown in Figure 5.

Figure 5. Schematic Potential Energy Surface for cis-butene-2 isomerization and decomposition (relative energies with zero point energy corrections expressed in kcal/mol). Adapted from Fig. 1 of Bedanov et al. 8


9


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 29

For the parameters needed to model this system, we used most of the same principal sources as Bedanov et al.8 Harmonic vibration frequencies for cis- and trans-butene-2 were taken from Richards and Nielsen22 and rotational constants computed at the B3LYP/aug-cc-pVTZ level of theory were obtained from the NIST Computational Chemistry Comparison and Benchmark Database.23 The harmonic frequencies for the cis/trans isomerization transition state were taken from Lin and Laidler24 and the rotational constants were assumed to the be the same as those for cis-butene-2. The harmonic vibrations for the decomposition transition state were taken from Alfassi et al.25 Two of the rotational constants were assumed to be the same as the smallest rotational constants of cis-butene-2 and the third was taken as the largest rotational constant in cis-butene-2 divided by a factor of 3.25 (extracted from Alfassi et al.). Relative energies were obtained from Bedanov et al. Simulations were carried out with initial conditions of pure cis-butene-2 with a 600 K thermal internal energy distribution immersed in a bath gas at a translational temperature of 1260 K (i.e. a shock simulation like the one described above). The energy transfer parameter was assumed8 to be α = 2000 cm–1 for collisions with the walls of the Knudsen cell used for Very Low Pressure Pyrolysis (VLPP). The pressure of the butene-2 collider gas was adjusted to give an average wall collision frequency of ~104.3 s–1, which corresponds to the VLPP experiments.25 The computed reaction flux coefficients for the three reactions are shown as a function of time in Figure 6. The figure reflects the fact that several simulations over various periods of time were carried out in order to cover the extended simulated time period with the desired resolution. At the end of each simulation, stochastic noise increases, because relatively few butene molecules remain unreacted at the end. The MultiWell input data are given in the Supplementary Information.] In their work, Bedanov et al.8 reported "time dependent rate constants" (now usually termed4, 26 "reaction flux coefficients") only for the decomposition channel to produce butadiene + H2: r1,3(t). For an irreversible reaction such as this one, the reaction flux coefficient can be identified with the rate constant when an SED exists. Interestingly, Bedanov et al. found that r1,3(t) exhibits two plateaus as a function of time. This behavior is confirmed by the current results shown in Figure 6. The current results are very similar to the topmost curve reported by Bedanov et al. in their Fig. 7, confirming that present stochastic simulations and eigenvalue calculations employed by Bedanov et al. are in agreement. Bedanov et al. explained the first plateau in r1,3(t) as due to a quasi-steady-state that exists in the cis-Well while the trans-Well is still essentially empty and thus acts as a nearly irreversible sink for reactive flux. This condition had been previously discussed by Green et al. 15 and the present results confirm this interpretation. Reaction flux coefficients r1,2(t) and r1,3(t) show that the energy distribution in the cis-Well is almost steady, while r2,1(t) shows that the energy distribution in the trans-Well is far from steady. About 0.1% of the butene has fragmented (to form butadiene + H2) prior to the first plateau, which persists until ~30% has


10

Page 11 of 29

fragmented. Eventually, the populations in the cis- and trans- wells equilibrate with each other, resulting in the second plateau as true SEDs are approached. The second plateau is reached when ~75% has fragmented and persists for the rest of the simulation (97% fragmented at 1.06 s). In the absence of stochastic noise, the plateaus are probably not truly constant, but are expected to exhibit small trends with time. The reaction flux coefficient for fragmentation is ~4 s-1 on the first plateau and ~15 s-1 on the second plateau. In contrast, the "rate constant" obtained by assuming a first order process for the entire simulation is ~3.3 s-1, which is in error by factor of ~4.5. The existence of multiple stages in the reaction is reflected in Figure 7, which shows for the two wells. When a near-SED exists (i.e. during the second plateau), the reaction flux coefficient r1,3(t), can be identified with the rate constant for reaction from the cis-Well to the bimolecular products, because it is an irreversible reaction. For a reversible isomerization reaction, the individual reaction flux coefficients (like r1,2(t) or r2,1(t)) can never be identified directly with the forward or reverse rate constant, even when a true SED exists in both wells. This is because both the forward and reverse elementary rate constants depend algebraically on both reaction flux coefficients.2, 3, 14, 27, 28

Take notice! If the stochastic simulations had not been carried out beyond 0.01 s and if only r1,3(t) had been considered, the existence of the first plateau might have been mistakenly interpreted as a true SED. However, by examining the reaction flux coefficients for reactions originating from the trans-Well it is clear that NSED conditions exist during the first plateau.

Reaction Flux Coefficients (s–1)

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


104 103

VLPP: 600K –> 1260K cis-butene-2 system 0.015 mbar

trans –> cis

102

cis –> trans

101

cis –> butadiene + H2

100 0.1%

10-1 10-4

10-3

% Conversion to Butadiene + H2 3.7% 32%

10-2

81%

10-1

100

Time (s)

Figure 6. Reaction flux coefficients as functions of time.


11


30000 Average Vibrational Energy (cm–1)

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 29

VLPP: 600K –> 1260K, 0.015 mbar cis-butene-2 system 25000

11200 11100

20000

11000

trans-butene-2

10800

15000 10000 5000 10-4

trans-butene-2

10900

10700 10-2

cis-butene-2

10-1 Time (s)

100

cis-butene-2

10-3

10-2

10-1

100

Time (s)

Figure 7. Average vibrational energy in both wells as a function of time.

Alkyl Free Radical Isomerization and Fragmentation: A multi-well system Alkyl free radicals are produced when H-atoms are abstracted from alkanes, such as found in hydrocarbon fuels. Once formed, these radicals rapidly undergo unimolecular isomerization and fragmentation reactions, as well as bimolecular reactions with O2 and other species. These reactions are important in combustion systems at even moderate temperatures because the isomerization and fragmentation barriers are relatively low.7, 8 The system comprised of the alkyl radicals produced by abstracting an H-atom from 2-methylhexane provide a good example of isomerization among multiple wells and fragmentation via multiple channels. Viscolcz et al.29 computed energies, geometries, and vibrational frequencies for the six isomers (wells) and the 15 isomerization transition states in this system. Reversible isomerization takes place via 3-, 4-, 5-, 6-, and 7-member cyclic transition states. Fragmentation, which was not considered by Viscolcz et al., takes place via C–C and C–H bond fission reactions. The reaction paths are illustrated schematically in Figure 8, where the index numbers for the isomers are based on the position of the free radical center and the index numbers for the product sets are arbitrary.30 Barker and Ortiz30 formulated a master equation model for the system, consisting of 6 isomers, 49 reaction paths, and 14 unique sets of fragmentation products by using RRKM theory31 to compute microcanonical rate constants, k(E)s, for the isomerization reactions and the Inverse Laplace Transform (ILT) method31 to compute k(E)s for the fragmentation reactions. The high-pressure limit Arrhenius parameters needed for the ILT calculations were estimated for the fragmentation channels by Barker and Ortiz, based on the work of Yamauchi et al.32 The present work utilizes the Barker and Ortiz model with Argon collider gas. All of the


12

Page 13 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


parameters needed to describe the reaction system have been published previously29, 30, 32 and are collected in the Supporting Information. (Note that because of revisions to the MultiWell Program Suite in the intervening years since 2001, the master equation results reported here are different in some ways from the earlier ones.30)

Figure 8. Reaction paths for isomerization (solid and broken lines) and fragmentation (arrows) in the 2-methylhexyl reaction system (adapted from Figure 2 in Barker and Ortiz.30 Isomerization reactions are shown for cyclic transition states via 6-member rings (thick solid lines), 5- and 7-member rings (thin solid lines), 4-member rings (long broken lines), and 3-member rings (short broken lines). Shock Experiment Simulations As described previously, the simulation was initiated by setting the initial internal (i.e., active energy) temperature Tvib equal to the ambient laboratory temperature (typically 300 K) and the translational temperature Ttrans equal to the final shock temperature (e.g., 1000 K).33-35 The initial energy distribution of the reactant active energy was assumed to be canonical (thermal) at Tvib. All isomerization reactions were treated as reversible and all fragmentation reactions were treated as irreversible. At t=0, all of the population was initially assigned to one of the wells. As the simulation proceeded, population moved from well to well (isomerization) and to fragmentation products. The pressure following the shockwave was assumed to be 1 bar. For the present discussion, the initial reactant was Well-01 (2-methyl-hex-1-yl free radical), but simulations initiated using other wells produced qualitatively similar results. The simulation was carried out using 106 stochastic trials for a simulated time corresponding to 2×104 (20k) collisions with argon bath gas (simulated duration of 3.66 µs). The number of collisions is emphasized, because collisional relaxation depends on the number of collisions.


13


The simulated concentrations (fractions) in the six wells are presented in Figure 9 as a function of the number of collisions. At t=0, all of the concentration resides in Well-01 (i.e. "fraction" f1(t=0) = 1) and the average vibrational energy, = 1792 cm-1, corresponds to Tvib = 300 K. As can be seen in Fig. 10, increases monotonically after the shock, until is reaches a steady value after ~5k collisions. Meanwhile, f1(t) has decreased by almost two powers of ten, but the concentrations in the other wells have increased from zero to finite values. In all of the wells except for Well-01, the concentrations rise to a maximum and then decay due to fragmentation reactions. After long periods of time, the surviving concentrations approach steady ratios, indicating that the surviving concentrations are equilibrating under the conditions of the simulated experiment. The "total" concentration in the wells (∑ fj(t) for the wells) is initially unity and remains unity for a short time (the incubation time) before beginning to decay approximately exponentially, as shown by the solid line in Figure 9. After 5k collisions, the total concentration in the wells has decreased by more than half; by the end of the simulation, it has decreased to 80% of the total C7H15 has already fragmented into smaller species; the yield of Product-08 would have been significantly over-estimated if NSEDs had been neglected. At steady-state, all three of the irreversible reactions shown have similar rate constants (~106 s–1). All three reaction flux coefficients start at zero and show large excursions immediately following the shock. The excursions decay as the NSEDs evolve toward steady state and the flux coefficients approach the values of the respective rate constants. Flux coefficient r1,7 grows from zero and approaches ki,7 as the NSED evolves toward the steady-state SED; r3,8 shows a large positive excursion (r3,8 ≈ ~3×k3,8); and r2,11 approaches k2,11 more quickly than the other two. Yet another way to evaluate this system is to define "effective" rate coefficients ki,e j 0,τ ; see Appendix, Eq. A7. The effective rate coefficients are averages of the

( )

flux coefficients over time. As an illustration, three effective rate coefficients are


16

Page 17 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


given in Figure 12. Consider the effective rate constant corresponding to r3,8: e k3,8 0,τ . If τ is defined as the time needed for 80% of the free radicals to dissociate,

( )

(

)

e the average rate constant during that period is k3,8 0,80% ≈ 1.4 ×106

e k3,8 (0,80% ) ≈ 1.4 ×106 s–1, which is about 40% larger than the true elementary rate

(

)

e constant k3,8. During the same period, k2,11 0,80% is almost the same as the true

(

)

e elementary rate constant, while k1,7 0,80% is about 20% lower than the true value.

Thus, the sensitivities of individual reactions to NSEDs in the same reaction system vary significantly. In calculations initiated in Well-03 (not shown) the sensitivities are different from those shown in Fig. 12, which are for calculations initiated in Well-01, indicating, as is expected, that the NSEDs and individual reaction sensitivities also depend on the way the excitation is initiated.

Figure 12. Individual reaction flux coefficients and effective rate coefficients (Eq. A7, see Appendix) for selected reactions as a function of the fraction of total C7H15 that has fragmented following the shock (the total elapsed simulated time is 3.66 µs). Note that all three reactions involve C–C bond fission and have similar rate constants.

1-Naphthyl+O2: A large aromatic system As an example of a chemical reaction system of a larger complexity we consider here the oxidation of 1-naphthyl by molecular oxygen studied recently.36 The studied system has a total of 18 wells and 5 products. The initial part of the potential energy surface (PES) is shown in Figure 13. The rest of the PES and further details on potential-energy and MultiWell calculations can be found in the reference cited


17


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 29

above. Here we examine this system for shock-initiated decomposition of 1naphthylperoxy at a pressure of 1 atm by performing MultiWell calculations at translational temperature of 1500 K and the initial vibrational temperature of 300 K.

Figure 13: Initial part of the potential energy surface for the 1-naphthyl+O2 system.36

Figure 14 depicts computed species fractions and Figure 15 the corresponding species vibrational energies. These results demonstrate that vibrational equilibration is attained on a finite time scale, the same phenomenon as demonstrated in the previous examples. Also demonstrated, in Figure 16, is that the reaction flux coefficients can serve as clear indicators of the attainment of the equilibrium, as discussed above. In the present example, the equilibrium is attained at about half a microsecond. This time becomes longer at lower temperatures and/or for different reaction channels. Studies indicate that polyaromatic systems exhibit NSED behavior rather often, under practical conditions, and with practical implications. For instance, in a recent study37 of a chemically-activated reaction of phenanthrene radical with OH, the overall reaction rate appeared to be increasing with time up to about 1 ms at 1500 K, i.e., on the same time scale as the oxidation of soot particles by OH measured in hydrocarbon flames.


18

Page 19 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 14: Population fractions versus time for the 1-naphthyl+O2 system

Figure 15: Average vibrational energy versus time for the 1-naphthyl+O2 system

Figure 16: Sum of reaction flux coefficients versus time for the wells in the 1naphthyl+O2 system


19


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 29

Conclusions In models of chemical systems, NSEDs are not always important. When collisional energy equilibration and chemical transformations are well separated in time, and consequently the IEREs and CSEs are clearly separable and identifiable, then all of the popular SED-based treatments suffice and all produce acceptable phenomenological rate coefficients that are essentially indistinguishable from elementary rate constants. However, all of the examples presented in the present work support the seminal work of Tsang and coworkers, who pointed out the ubiquitous existence of NSED conditions. In some cases NSED conditions can be neglected, but in other cases the presence of NSEDs can quantitatively affect rate constants and branching ratios predicted using SED-based methods. In addition, the examples presented here suggest that the effects of NSED conditions may be more significant in larger species and in reaction systems that are more complex. The presence of NSEDs presents two problems. The first problem is the principal issue addressed in the present work: the convenient detection of NSEDs. The second problem posed by the presence of NSEDs is the development and introduction of practical methods accounting for NSEDs in process simulations codes, such as ChemKin. Two measures of the energy distribution for each species are particularly useful for detecting the presence of NSEDs: the average vibrational energy E vib ( t ) , the first moment of the energy distribution; and ri(t), the reaction flux coefficient, which is the average of the energy distribution, weighted by ki(E). The sum of reaction flux coefficients originating in a single well is also useful. All of these are good diagnostic indicators for NSEDs. The reaction flux coefficients are particularly useful, however, because for irreversible reactions in the limit where SED conditions exist they can be identified with reaction rate constants. In those cases, they can be used to quantitatively assess the errors associated with the presence of NSED conditions. The second problem introduced by NSEDs is developing practical methods for chemical modeling. Currently, the most accurate method for modeling a system affected by NSED conditions is to use a master equation to calculate time dependent concentration profiles during the period when NSEDs exist and an elementary reaction mechanism with time independent rate coefficients when SED conditions have been attained to a sufficient degree. Rate constants by themselves are not enough. All of the existing methods for solving the ME can produce the same time dependent concentration profiles: the stochastic method simulates them algorithmically and the matrix solution does the same by including all eigenmodes. The power of the ME simulation approach is that it can accurately model any manifestation of NSED conditions and can make predictions, not just provide descriptions of unusual phenomena.


20

Page 21 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Although feasible and very powerful, master-equation modeling is too cumbersome and computer-intensive at present for practical applications. Phenomenological models derived from master equation solutions go a long way toward pragmatic implementation in detailed reaction models, even when NSEDs exist.1, 4, 5, 38 Because real chemical systems of practical importance are often topologically highly complex, the accuracy of phenomenological models probably cannot be assessed just on the basis of the magnitudes of their eigenvalues, especially when one wishes to use them under non-ideal conditions. The best way to assess their accuracy is to compare their predictions to those of the master equation on which they are based. Just as it is good practice to test other kinetics approximations in individual applications (e.g. the pseudo steady state approximation), the approximate phenomenological models derived from master equations should also be tested. In our opinion, methods that somehow add NSED effects to existing elementary reaction models may be useful in the interim for describing NSED phenomena, but ultimately they are likely to lack the flexibility and rigor needed for predictive theory. Instead of focusing just on interim solutions, the focus should be shifted to inventing new methods and writing new codes that properly account for this important physical reality.

Acknowledgements JRB thanks the National Science Foundation (Atmospheric and Geospace Sciences) and NASA (Upper Atmospheric Research Program) for financial support. MF acknowledges financial support by the US Army Corps of Engineers, Humphreys Engineering Center Support Activity (Contract No. W912HQ-11-C-0035 to UCB). DMG acknowledges support from the Department of Mechanical Engineering, Stanford University. We thank Struan Robertson for valuable discussions about the cis-butene-2 simulations. We also are very happy to acknowledge our valuable associations with Joe Michael, Al Wagner, Larry Harding, the shock tube group, and other researchers at Argonne National Laboratory, who have been at the forefront of scientific developments in the combustion field for many years. Supporting Information Available: Input files for all MultiWell calculations referred to in the text. This includes vibrational frequencies, moments of inertia, heats of formation and energy transfer parameters. This material is available free of charge via the Internet at http://pubs.acs.org.

Appendix Master Equation A master equation (ME) describes chemical reaction systems in which physical and chemical transformations occur on similar timescales and interact with each other.


21


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 29

For unimolecular reaction systems consisting of isomerization and fragmentation reactions, collisional energy transfer is usually the most important physical process, but light absorption and spontaneous photon emission are sometimes important.19 Three key approximations are usually made in formulating a master equation: (a) all translational velocity distributions are thermal (Maxwell-Boltzmann); (b) collisionfree intramolecular vibrational energy redistribution (IVR) is so fast that the "active energy" in each isolated molecule is distributed statistically (i.e., the internal energy distribution is microcanonical); (c) the microcanonical unimolecular reaction rate constant k(E,J) for each reaction (for a specified angular momentum) is described by statistical rate theory. (It is worth mentioning that master equation methods have also been applied to conditions of slow IVR.39 For molecules with very limited numbers of internal states (i.e., diatomics and some triatomics), the master equation can be set up and solved as a set of coupled ordinary differential equations (ODEs), one equation for each quantum state. For larger molecules at high internal energies, there are too many quantum states and this approach is intractable. Instead, an “energy grained” master equation is formulated by dividing the energy scale into multiple small energy bins (i.e., energy grains) and treating the states in each energy grain as a group. A closely related variation is to treat the quantum states as a quasi-continuum characterized by the density of states ρ(E,J). Because the ME depends on both E and J, it is very cumbersome to solve explicitly31, and common practice is to reduce this two-dimensional (2D) ME to one dimension (1D) by invoking various approximations.11, 31 The resulting 1D ME can be solved either as a coupled set of ODEs (typically by eigenvalue methods42-44 or as stochastic simulations19, 42, 44, 45 Gillespie proved that properly formulated stochastic simulation and the corresponding differential equation theoretically converge to the same solution as the number of stochastic trials increases.12, 13, 18 The stochastic approach, as implemented in previous stochastic simulation codes19, 46-48 and in the MultiWell Program Suite,20, 49, 50 was used for all calculations reported here, but the results in principle can be replicated using eigenvalue methods. Thus the present results and conclusions do not depend on how the calculations were carried out. 40, 41

The MultiWell ME is described in the MultiWell User Manual and elsewhere. The reduction of the 2D ME is based on centrifugal corrections51, 52 resulting in a 1D ME that depends on the “active energy”, ε. In the following expressions, primed ε' refers to the current active energy and the unprimed quantity refers to a different active energy: ∞

channels dN (ε ′,t) = F ( ε ′,t ) + ω ∫  P ( ε ′; ε ) N ( ε ,t ) − P ( ε ; ε ′ ) N ( ε ′,t )  d ε − N ( ε ′,t ) ∑ ki ( ε ′ ) dt i=1 0

J

(A1)


22

Page 23 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


where N(ε',t)dε ' is the concentration of a chemical species with active energy in the range ε' to ε'+dε'; ω is the inelastic collision frequency and P ε , ε ′ P ε , ε ′ is the normalized step-size distribution for collisional energy transfer from initial energy ε' to energy ε; F ε ′,t F ε ′,t is a source term (e.g., thermal, chemical, photo activation, or isomerization); and ki(ε') is the unimolecular reaction rate constant for molecules at energy ε' (averaged over the rotational distribution52)52) reacting via the ith reaction channel. For simplicity, an index number denoting the chemical species has been omitted.

(

) (

)

( ) ( )

The inelastic collision frequency ω is governed the temperature, collider gas concentration, and by a Lennard-Jones intermolecular potential between the reactant and the collider gas. The conventional exponential-down step size distribution with energy transfer parameter α was assumed when computing P ε ,ε ′ P ε ,ε ′ .

(

) (

)

All of the parameters for each system are tabulated in Supporting Information. See the original papers for discussion of the selection of the parameters. Different choices of parameters would result in quantitative differences in the simulations, but we believe that all realistic choices will give results that are qualitatively similar and will lead to the same qualitative conclusions. A plethora of terms have been used to describe the time and energy dependent concentrations and quantities derived from them. To avoid ambiguous language, we define the following terms: (a) Population or fraction:

f (t ) =

∫

∞

0

N ( ε ′, t ) d ε ′

∫

∞

0

N ( ε ′, 0 ) d ε ′

(A2)

(b) Energy distribution; when varying with time, a non-Steady Energy Distribution (NSED)

D ( ε ′, t ) = N ( ε ′, t )

∫

∞

0

N ( ε ′, t ) d ε ′

(A3)

After equilibration with the bath gas, this becomes a Steady Energy Distribution

(

) (

)

(SED): D ε ′,∞ D ε ′,∞ . (d) Average active energy (also referred to as the “average vibrational energy”):

Evib ( t ) =

∫

∞

0

ε ′D ( ε ′, t ) d ε ′

(A4)

(e) Reaction flux coefficient for the ith reaction channel: ∞

ri ( t ) = ∫ D ( ε ′, t ) ki ( ε ′ ) d ε ′

(A5)

0

(f) Incubation time: the induction period during which energy relaxation takes place; associated with loss of a reactant, or production of a product.


23


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 29

Effective Rate Constants for Irreversible Reactions The reaction flux coefficient ri,j(t) for an irreversible reaction at steady-state (i.e. with SED) is identified with the corresponding rate constant ki,j. When not at steadystate (i.e. when NSEDs exist), the reaction flux coefficients vary with time. As the NSED evolves with time, it asymptotically approaches the steady-state SED. Thus for an irreversible reaction, one can write

ki, j = lim ri, j ( t )

(A6)

t→∞

an "effective" rate constant for a specified time period can be obtained by averaging ri,j(t)over time. For the time period t0 → τ, an effective rate constant ki,e j t0 ,τ can be

( )

defined as the average over time:

ki,e j ( t0 ,τ ) = (τ − t0 )

−1

∫

τ

r

t0 i, j

( t ) dt

(A7)

( )

In the limit as τ→∞, ki,e j t0 ,τ becomes equal to the rate constant:

lim ki,e j ( t0 ,τ ) = ki, j

(A8)

τ →∞

In the above, t0=0 is the natural choice for shock experiments and other phenomena initiated at t=0.


24

Page 25 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


References 1. Bartis, J. T.; Widom, B., Stochastic Models of the Interconversion of Three or More Chemical Species. J. Chem. Phys. 1974, 60, 3474-3482. 2. Widom, B., Molecular Transitions and Chemical Reaction Rates. Science 1965, 148, 1555-1560. 3. Quack, M., On the Mechanism of Reversible Unimolecular Reactions and the Canonical (“High Pressure”) Limit of the Rate Coefficient at Low Pressures. Ber. Bunsenges. Phys. Chem. 1984, 88, 94-100. 4. Miller, J. A.; Klippenstein, S. J., Master Equation Methods in Gas Phase Chemical Kinetics. J. Phys. Chem. A 2006, 110 (36), 10528-10544. 5. Glowacki, D. R.; Liang, C.-H.; Morley, C.; Pilling, M. J.; Robertson, S. H., Mesmer: An Open-Source Master Equation Solver for Multi-Energy Well Reactions. J. Phys. Chem. A 2012, 116 (38), 9545-9560. 6. Miller, J. A.; Klippenstein, S. J., Determining Phenomenological Rate Coefficients from a Time-Dependent, Multiple-Well Master Equation: ‘‘Species Reduction’’ at High Temperatures. Phys. Chem. Chem. Phys. 2013, 15, 4744-4753. 7. Tsang, W.; Bedanov, V.; Zachariah, M. R., Master Equation Analysis of Thermal Activation Reactions: Energy-Transfer Constraints on Falloff Behavior in the Decomposition of Reactive Intermediates with Low Thresholds. J. Phys. Chem. 1996, 100 (10), 4011-4018. 8. Bedanov, V. M.; Tsang, W.; Zachariah, M. R., Master Equation Analysis of Thermal Activation Reactions: Reversible Isomerization and Decomposition. J. Phys. Chem. 1995, 99, 11452- 11457. 9. Kee, R. J., Miller, J.A., Jefferson, T.H. Chemkin: A General Purpose, Problem Independent,Transportable Fortran Chemical Kinetics Code Package, Sandia National Laboratories: SAND80-8003, 1989. 10. Oppenheim, I.; Shuler, K. E.; Weiss, G. H., Stochastic Processes in Chemical Physics: The Master Equation. MIT Press: Cambridge, MA, 1977. 11. Gilbert, R. G.; Smith, S. C., Theory of Unimolecular and Recombination Reactions. Blackwell Scientific: Oxford, 1990. 12. Gillespie, D. T., A General Method for Numerically Simulating the Stochastic Time Evolution of Coupled Chemical Reactions. J. Comp. Phys. 1976, 22 (4), 403-434. 13. Gillespie, D. T., Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys. Chem. 1977, 81 (25), 2340-2361. 14. Rice, O. K., On the Relation between an Equilibrium Constant and the NonEquilibrium Rate Constants of Direct and Reverse Reactions. J. Phys. Chem. 1961, 65, 1972-1976. 15. Green, N. J. B.; Marchant, P. J.; Perona, M. J.; Pilling, M. J.; Robertson, S. H., Forward and Reverse Rate Coefficients in Equilibrating Isomerization Reactions. J. Chem. Phys. 1992, 96, 5896-5907. 16. Miller, J. A.; Klippenstein, S. J.; Robertson, S. H.; Pilling, M. J.; Green, N. J. B., Detailed Balance in Multiple-Well Chemical Reactions. Phys. Chem. Chem. Phys. 2009, 11, 1128 - 1137.


25


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 29

17. Gillespie, D. T., A Rigorous Derivation of the Chemical Master Equation. Physica A 1992, 188, 404-425. 18. Gillespie, D. T., Markov Processes: An Introduction for Physical Scientists. Academic Press: San Diego, CA, 1992. 19. Barker, J. R., Monte-Carlo Calculations on Unimolecular Reactions, EnergyTransfer, and Ir-Multiphoton Decomposition. Chem. Phys. 1983, 77 (2), 301-318. 20. Barker, J. R.; Ortiz, N. F.; Preses, J. M.; Lohr, L. L.; Maranzana, A.; Stimac, P. J.; Nguyen, T. L.; Kumar, T. J. D. Multiwell-2014.1b Software, J. R. Barker, University of Michigan: Ann Arbor, Michigan, USA (http://aossresearch.engin.umich.edu/multiwell/), 2014. (accessed June 25, 2014). 21. Miller, J. A.; Klippenstein, S. J., The Reaction between Ethyl and Molecular Oxygen. 2: Further Analysis. Int. J. Chem. Kinet. 2001, 33 (11), 654-668. 22. Richards, C. M.; Nielsen, J. R., Raman Spectra of Cis- and Trans-2-Butene in the Gaseous and Liquid States. J. Opt. Soc. Am. 1950, 40, 442-445. 23. NIST Computational Chemistry Comparison and Benchmark Database, Nist Standard Reference Database Number 101. http://srdata.nist.gov/cccbdb (accessed November 12, 2014). 24. Lin, M. C.; Laidler, K. J., Theoretical Aspects of Gas-Phase Thermal Isomerizations. Trans. Faraday Soc. 1968, 64, 94-102. 25. Alfassi, Z. B.; Golden, D. M.; Benson, S. W., Very Low-Pressure Dehydrogenation of Cis-2-Butene. Activation Energy for 1,4-Hydrogen Elimination. Int. J. Chem. Kinet. 1973, 5 (6), 991-1000. 26. Aguda, B. D.; Pritchard, H., Reversible and Irreversible Formulation of Unimolecular Reactions. J. Chem. Phys. 1992, 96, 5908-5914. 27. Rice, O. K., Further Remarks on the "Rate Quotient Law”’. J. Phys. Chem. 1963, 67, 1733-1735. 28. Snider, N. S., Nonequilibrium Effects in the Kinetics of Isomerization Reactions and in the Kinetics of Dissociation and Recombination of Diatomic Molecules. J. Chem. Phys. 1965, 42, 548-555. 29. Viskolcz, B.; Lendvay, G.; Seres, L., Ab Initio Barrier Heights and Branching Ratios of Isomerization Reactions of a Branched Alkyl Radical. J. Phys. Chem. A 1997, 101, 7119-7127. 30. Barker, J. R.; Ortiz, N. F., Multiple-Well, Multiple-Reaction-Path Unimolecular Reaction Systems. Ii. 2-Methylhexyl Free Radicals. Int. J. Chem. Kinet. 2001, 33 (4), 246-261. 31. Forst, W., Unimolecular Reactions. A Concise Introduction. Cambridge University Press: Cambridge, 2003. 32. Yamauchi, N.; Miyoshi, A.; Kosaka, K., Thermal Decomposition and Isomerization Processes of Alkyl Radicals. J. Phys. Chem. A 1999, 103, 2723-2733. 33. Shi, J.; Barker, J. R., Incubation in Cyclohexene Decomposition at High Temperatures. Int. J. Chem. Kinet. 1990, 22, 187-206. 34. Barker, J. R.; King, K. D., Vibrational Energy Transfer in Shock-Heated Norbornene. J. Chem. Phys. 1995, 103, 4953-4966. 35. Barker, J. R.; Stimac, P. J.; King, K. D.; Leitner, D. M., CF3CH3 → HF + CF2CH2: A Non-RRKM Reaction? J. Phys. Chem. A 2006, 110, 2944-2954.


26

Page 27 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


36. Kislov, V. V.; Singh, R. I.; Edwards, D. E.; Mebel, A. M.; Frenklach, M., Rate Coefficients and Product Branching Ratios for the Oxidation of Phenyl and Naphthyl Radicals: A Theoretical RRKM - Master Equation Study. Proc. Combust. Inst. 2015, 35, 1861-1869. 37. Edwards, D. E.; Zubarev, D. Y.; Lester, W. A.; Frenklach, M., Pathways to Soot Oxidation: Reaction of OH with Phenanthrene Radicals. J. Phys. Chem. A 2014, 118 (37), 8606-8613. 38. Georgievskii, Y.; Miller, J. A.; Burke, M. P.; Klippenstein, S. J., Reformulation and Solution of the Master Equation for Multiple-Well Chemical Reactions. J. Phys. Chem. A 2013, 117, 12146−12154. 39. Weston, R. E., Jr.; Barker, J. R., On Modeling the Pressure-Dependent Photoisomerization of Trans-Stilbene by Including Slow Intramolecular Vibrational Redistribution. J. Phys. Chem. A. 2006, 110, 7888-7897. 40. Venkatesh, P. K.; Dean, A. M.; Cohen, M. H.; Carr, R. W., Master Equation Analysis of Intermolecular Energy Transfer in Multiple-Well, Multiple-Channel Unimolecular Reactions. I. Basic Theory. J. Chem. Phys. 1997, 107 (21), 8904-8916. 41. Venkatesh, P. K.; Dean, A. M.; Cohen, M. H.; Carr, R. W., Master Equation Analysis of Intermolecular Energy Transfer in Multiple-Well, Multiple-Channel Unimolecular Reactions. Ii. Numerical Methods and Application to the Mechanism of the C2H5+O2 Reaction. J. Chem. Phys. 1999, 111 (18), 8313-8329. 42. Widom, B., Collision Theory of Chemical Reaction Rates. Advan. Chem. Phys. 1963, 5, 353-386. 43. Hoare, M., Steady-State Unimolecular Processes in Multilevel Systems. J. Chem. Phys. 1963, 38 (7), 1630-1635. 44. Tardy, D. C.; Rabinovitch, B. S., Collisional Energy Transfer. Thermal Unimolecular Systems in the Low Pressure Region. J. Chem. Phys. 1966, 45 (10), 3720-3730. 45. Montroll, E. W., The Application of the Theory of Stochastic Processes to Chemical Kinetics. Adv. Chem. Phys. 1958, 1, 361-399. 46. Barker, J. R., Infrared Multi-Photon Decomposition - a Comparison of Approximate Models and Exact-Solutions of the Energy-Grained Master Equation. J. Chem. Phys. 1980, 72 (6), 3686-3694. 47. Barker, J. R., Direct Measurements of Energy Transfer Involving Large Molecules in the Electronic Ground State. J. Phys. Chem. 1984, 88, 11-18. 48. Vereecken, L.; Huyberechts, G.; Peeters, J., Stochastic Simulation of Chemically Activated Unimolecular Reactions. J. Chem. Phys. 1997, 106 (16), 65646573. 49. Barker, J. R., Multiple-Well, Multiple-Reaction-Path Unimolecular Reaction Systems. I. Multiwell Computer Program Suite. Int. J. Chem. Kinet. 2001, 33 (4), 232245. 50. Barker, J. R., Energy Transfer in Master Equation Simulations: A New Approach. Int. J. Chem. Kinet. 2009, 41, 748-763. 51. Marcus, R. A., Dissociation and Isomerization of Vibrationally Excited Species. Iii. J. Chem. Phys. 1965, 43 (8), 2658-2661.


27


1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 29

52. Weston, R. E., Jr.; Nguyen, T. L.; Stanton, J. F.; Barker, J. R., HO + CO Reaction Rates and H/D Kinetic Isotope Effects: Master Equation Models with Ab Initio SCTST Rate Constants. J. Phys. Chem. A 2013, 117, 821-835.


28

Page 29 of 29

1 2 3 4 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 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table of Contents Image


29

Cosmetic doping--when anabolic-androgenic steroids are not enough.

When normal is not enough.

When novelty is not enough.

When is enough enough? We are only beginning.

Augmentation genioplasty: when bone is not enough.

When is enough, enough?

When There Is Not Enough Evidence and When Evidence Is Not Enough: An Australian Indigenous Smoking Policy Study.

Habitat and Vegetation Variables Are Not Enough When Predicting Tick Populations in the Southeastern United States.

Differences in metacognitive regulation in introductory biology students: when prompts are not enough.

Venereal disease education: facts are not enough.

Cancer cachexia-when proteasomal inhibition is not enough.

Noninferior Antibiotics: When Is "Not Bad" "Good Enough"?

Gout treatment: when a tablet is not enough.

Communicating about communicating: when innate is not enough.

Treatment of gonococcal infection: when one drug is not enough.

Floral Organogenesis: When Knowing Your ABCs Is Not Enough.

When enough is not enough: Information overload and metacognitive decisions to stop studying information.

Evidence about inotropes: when is enough, enough?

EU air quality limits are not stringent enough, say researchers.

Understanding the poor iontophoretic transport of lysozyme across the skin: when high charge and high electrophoretic mobility are not enough.

Correction: When Optimal Feedback Control Is Not Enough: Feedforward Strategies Are Required for Optimal Control with Active Sensing.

When is enough enough? Ethical dilemmas on the burn unit.

The clinical psychologist as program consultant: When is enough enough?

When good isn't good enough.