Subscriber access provided by ORTA DOGU TEKNIK UNIVERSITESI KUTUPHANESI
Comment
Comment on “When Rate Constants Are Not Enough” by John R. Barker, Michael Frenklach, and David M. Golden James A. Miller, Stephen J Klippenstein, Struan H. Robertson, Michael J. Pilling, Robin J. Shannon, Judit Zádor, Ahren W Jasper, C. Franklin Goldsmith, and Michael P. Burke J. Phys. Chem. A, Just Accepted Manuscript • DOI: 10.1021/acs.jpca.5b06025 • Publication Date (Web): 23 Dec 2015 Downloaded from http://pubs.acs.org on December 23, 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 14
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
Comment on “When Rate Constants Are Not Enough” by John R. Barker, Michael Frenklach, and David M. Golden James A. Miller,1* Stephen J. Klippenstein,1 Struan H. Robertson,2 Michael J. Pilling,3 Robin Shannon,3 Judit Zádor,4 Ahren W. Jasper,4 C. Franklin Goldsmith,5 and Michael P. Burke6 1
Chemical Sciences and Engineering Division, Argonne National Laboratory, Argonne, IL 60439, USA 2 Dassault Systèmes, BIOVIA, Science Park, Cambridge, CB4 0WN, UK 3 School of Chemistry, University of Leeds, Leeds, LS2 9JT, UK 4 Combustion Research Facility, MS 9055, Sandia National Laboratories, Livermore, CA 94551-0969, USA 5 School of Engineering, Brown University, Providence, RI 02912, USA 6 Department of Mechanical Engineering, Department of Chemical Engineering, and Data Sciences Institute, Columbia University, New York, NY 10027, USA
*
[email protected] Abstract We discuss “phenomenological” reactions, rate constants, and flux coefficients in the context of their relationship to master-equation methods for obtaining rate constants. The isomerizationdissociation of 2-butene is used as an example for illustrating some of the pitfalls that can exist in trying to obtain rate constants from population-time histories alone (with associated, derivative quantities such as flux coefficients). We attach a considerable amount of Supporting Information that elaborates on some of the key points made in the text. Included in the Supporting Information is a description of two methods for obtaining rate constants from the master equation when its transition matrix is not readily available.
Comment Introduction In their paper titled “When Rate Constants are Not Enough”, Barker, Frenklach, and Golden (BFG)1 raise an important question. A key aspect of this question is the identification of conditions under which a rate-constant description in chemical kinetics is not sufficient. Their discussion centers around reactions that involve a single or multiple potential energy wells, single or multiple products, and an interaction between collisional energy transfer with bath gas molecules and chemical reaction. This type of problem is usually treated theoretically using a master equation (ME) approach, in which the time evolution of an energy-resolved molecular system is described by a set of coupled, linear differential equations (or integro-differential equations) that incorporate both energy transfer and chemical reaction. BFG argue that rate constants for the macroscopic reactions modeled cannot be time independent if the energy distribution of reactants and intermediates in the wells is not in a steady state, thus rejecting one of the foundational concepts 1 ACS Paragon Plus Environment
The Journal of Physical Chemistry
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
of chemical kinetics. They use the acronym NSED (non-steady-state energy distribution) to describe the populations under these conditions. The purpose of this Comment is to demonstrate that it is not necessary to abandon the concept of time-independent rate coefficients and that no revision of this foundational concept is required to describe these systems accurately. We discuss an established methodology for calculating time-independent rate constants that evolves naturally from an alternative method of solution of the ME based on a matrix approach. Consider the transition matrix of an energy- and species-resolved, multiple-well master equation describing the collisionally- and chemically-induced transitions in molecules or complexes of the system on a microscopic scale. The eigenvectors of this transition matrix generally separate into a set with small (in magnitude) eigenvalues describing the chemical transformations and a set with larger-magnitude eigenvalues describing the energy relaxation process. A rate constant description, which is formulated to describe only the chemical transformations on a macroscopic scale, is then appropriate for all times significantly beyond that characterizing the slowest decaying energy relaxation process (normally a small multiple of that time). This is an absolute correlation that needs no lengthy discussion. Nevertheless, we discuss this point for an example problem in the Supporting Information. Unfortunately, BFG confuse the situation by first addressing the issue of collisional and chemical time scale separation in very general terms, with which we basically agree. But then they present examples, and discussion of those examples, that instead focus primarily on effects that occur well beyond the energy-relaxation timescale. These cases can clearly be described by rate constants, as we demonstrate in this Comment. Their confusion is largely related to an unwarranted emphasis on steady and non-steady population distributions over energy levels (SEDs and NSEDs). Their intent apparently is to exploit SED conditions to calculate rate constants, as flux coefficients in some cases. We show that this can lead to rate coefficients whose numerical values are incorrect, and thus are unsuitable to be used in chemical kinetic models. In contrast, we present three distinct methods to obtain rate coefficients appropriate for kinetic modeling. The BFG article cries out for a clear, unambiguous definition of a rate constant (or rate coefficient). To this end it is convenient to refer to the IUPAC definition of an elementary reaction. Such a reaction is defined by IUPAC as one “for which no reaction intermediates have been detected or need to be postulated in order to describe the chemical reaction on a molecular scale. An elementary reaction is assumed to occur in a single step and to pass through a single transition state.”2 This definition, relying on a single transition state, is unnecessarily restrictive for our current purposes; many reactions involve a number of transition states. We wish to define a phenomenological reaction as an elementary reaction with the single-transition-state condition replaced by the condition that its rate constant be time independent and not depend on initial conditions over a very wide range of thermodynamic conditions (temperature, pressure, and possibly compositions). Other restrictions (related to molecularity, reaction order, etc.) that must be satisfied by phenomenological reactions are the same as those for any elementary reaction. We must also clarify what is meant by an intermediate - in this context an intermediate must be a readily identifiable molecule or radical, not just a transient complex that (perhaps) suffers a few collisions in a well and goes on to react. For example, the addition-elimination reaction of ethyl with O2,
2 ACS Paragon Plus Environment
Page 2 of 14
Page 3 of 14
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
C2H5 + O2 → C2H4 + HO2, is treated as a single reaction, even though the reactants pass through the C2H5O2 well and two transition states to get to products; C2H5O2 is not recognized in the phenomenological description of this reaction, because it is not actually formed as a thermodynamic species. This type of reaction is called “well-skipping” (perhaps the term “well-mediated” is less confusing, but “well-skipping” is in common usage, along with “formally direct”). Whether one uses well-skipping or wellmediated, the term applies only to the phenomenological description; at the microscopic level the molecular configuration of C2H5O2 in not skipped en route to products. We should point out that the phenomenological two-step sequence C2H5 + O2 → C2H5O2 C2H5O2 → C2H4 + HO2 occurs in parallel with the “well-mediated” reaction in chemical kinetics models. In this sequence, the peroxy radical is collisionally stabilized (thermalized), followed by collisional excitation and reaction; it is favored by higher pressures and the well-mediated reaction by lower pressures. It is important to understand what is meant by the term phenomenological in describing a rate coefficient. A reasonably good definition is that a phenomenological rate constant (or rate coefficient) is the one that goes into the phenomenological rate law. Critical to the development and application of these long-standing rate laws is the idea that the phenomenological rate coefficient for a reaction cannot be explicitly a function of time. Furthermore, a phenomenological rate coefficient is explicitly a function of macroscopic, thermodynamic variables, not microscopic mechanical variables (e.g. molecular constants of the motion). As a consequence, phenomenological rate coefficients in the present context are not necessarily local. For a reaction involving several potential energy wells, for example, the phenomenological reactions and their associated rate constants can be viewed as properties of the entire system under consideration. They can describe a reaction between configurations (wells or sets of bimolecular fragments) that are not adjacent to each other (well-skipping). They cannot, in general, be interpreted as the probability per unit time that a reactant (or reactants) will react to form products, i.e. they are not necessarily flux coefficients (defined explicitly in Eq. 4 below). These rate coefficients, k(T,p) (suppressing the possible dependence on composition), are to be contrasted with the microscopic rate coefficients, k(E,J) or k(E), where E is the vibrational-rotational energy in the reactant molecule and J is its total angular momentum quantum number. The rate coefficients k(E,J) differ from k(T,p) in that the former depend on microscopic, molecular constants of the motion, not thermodynamic variables. The k(E,J) are intrinsically local (no well-skipping), because they represent a coarse-grained, statistical description of the exact microscopic dynamics of the motion. They can be interpreted unambiguously as probabilities per unit time of reaction occurring – they are necessarily flux coefficients. In his seminal paper on the subject, Widom3 refers to k(T,p) as phenomenological because it encapsulates the consequences of a complex sequence of microscopic, molecular events, in some cases masking the details. This treatment of phenomenological reactions and their associated rate constants is exactly what is needed in macroscopic chemical kinetics models that are used, e.g., in atmospheric and combustion chemistry.
3 ACS Paragon Plus Environment
The Journal of Physical Chemistry
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 14
This brings us back to the BFG paper and its unwarranted emphasis on SEDs and NSEDs. Let us be clear. There is no essential connection between SEDs and rate constants (as we have defined them). In describing the CSE (chemically significant eigenmode, eigenvalue, eigenvector) method of obtaining rate constants, BFG indicate that in this method “an SED is considered to exist and time independent rate constants can be defined where the magnitudes of the CSEs are much smaller than those of the IEREs” (IEREs are internal energy relaxation eigenmodes, eigenvalues, or eigenvectors). SEDs play no role whatsoever (explicit or implicit) in the CSE methodology. The classic example illustrating this point is equilibration in the two-well system studied by Quack4 and others5. In Quack’s problem the single CSE (in addition to the equilibrium eigenmode) is separated by more than an order of magnitude from the IEREs. In this case neither the forward nor reverse flux coefficients are ever constant during the reaction, indicating that NSED conditions always exist during this time. Nevertheless, rate constants can be determined easily for the forward and reverse reactions, and their sum is equal to the modulus of the non-zero CSE. CSE methodology The CSE methodology6 is described in more detail in SI-I and in even more detail in the references. A rate-constant treatment of chemical kinetics (in the macroscopic sense) is possible only when chemical time scales can be separated from vibrational-rotational energy relaxation time scales. In a matrix formulation of the ME the reciprocals of these scales correspond to the negatives of the eigenvalues of the transition matrix. They generally separate into CSEs and IEREs; there are many IEREs and relatively few CSEs. As long as the magnitudes of the IEREs substantially exceed those of the CSEs, little or no chemical reaction occurs during the rotational-vibrational energy relaxation process. A discussion of how to handle cases when this condition is not satisfied is well beyond the scope of this Comment, but it can be found elsewhere.6f,7 It is, however, not at issue in the present discussion. Determining the rate constants rests on identifying the solution for
the macroscopic populations, X i xi ( E )dE for wells ( xi (E)dE is a suitable measure of the 0
population at energies between E and E + dE in well i), as
X i (t )
Nchem
ae j 0
jt
i=1,…,S
ij
(1)
where Nchem = S – 1 is the number of CSEs (excluding the λ0 = 0 equilibrium eigenmode). The coefficients aij contain information about the initial conditions and the jth eigenvector, corresponding to the λj eigenvalue. Each term in Eq.(1), other than the 0 0 term, describes how a different CSE contributes to bringing a particular species population to its equilibrium value from an arbitrary initial condition. Provided the CSEs are well separated in magnitude, each CSE describes the equilibration of one group of previously-equilibrated species with another (a group may be a single species). Eq. (1) is readily identified as the solution to a system of first-order (or pseudo first-order) macroscopic, phenomenological rate equations,
dX i (t ) kil X l X i kli dt l i l i
i 1,..., S
4 ACS Paragon Plus Environment
(2)
Page 5 of 14
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
The details are given elsewhere,6f, 7 but it is a straight-forward task to invert Eq. (1) and obtain the rate constants that go into Eq. (2) (see SI-I for one derivation). Note that Eq. (2) inherently contains all phenomenological reactions that skip wells (well-mediated reactions). They are essential to describe the correct time evolution of the system (see SI-III for an example of what happens when they are excluded) and evolve naturally from the master equation analysis. There are no approximations in this method other than the separation of scales (CSEs from IEREs), and this really is not an approximation. It is a real, general, physical constraint on the existence of rate constants and rate laws as we know them. The rate constants derived reproduce virtually the exact population histories, Xi(t), that come from the ME (after a short time delay) as long as no significant reaction occurs during the vibrational-rotational relaxation period. This point is illustrated below and in SI-II. We cannot overemphasize that rate constants are characteristic of the transition matrix (i.e. the entire system under consideration) and not of the population vectors (distributions) on which they act. To determine rate constants one never needs to integrate the ME in time – rate constants can always be determined directly from the eigenvalues and eigenvectors of the transition matrix6d, 6g (see also SI-I). Time integration is performed for one of 3 reasons: to simulate directly some laboratory experiment, to add physical insight, or because one does not have easy access to the transition matrix. Of course, using the stochastic ME approach, BFG do not have the eigenvalues and eigenvectors of the transition matrix at their disposal, clearly a weakness of the approach. One can, however, still extract the correct phenomenological rate coefficients in many cases, as demonstrated in SI-III and SI-IV.
Demonstration of the CSE methodology We shall use as an example the cis-butene-2 isomerization/dissociation problem discussed by BFG (see Fig. 5 in BFG for the PES, a variant of which is also given at the beginning of our Supporting Information) and studied previously by Bedanov, Tsang, and Zachariah.8 Our primary purpose in this example is to apply the CSE method to a simple problem and interpret the BFG results in light of the CSE description. The reactions involved are as follows: cis-butene-2 ⇄ trans-butene-2 cis-butene-2 → butadiene + H2 trans-butene-2 → butadiene + H2 The rate constants for forward and reverse isomerization reactions are labeled k21 and k12. The rate constants for the subsequent reactions are k31 and k32; the latter is a well-skipping reaction. The kinetics of this system of phenomenological reactions conforms to a set of three coupled differential rate equations (c.f. Eq. 2). This set of rate equations can be expressed in matrix form; there are two dependent variables (because the fragment concentration is treated as a sink, the differential equation for its concentration can always be decoupled from the rest), dX Kr X dt
(3)
5 ACS Paragon Plus Environment
The Journal of Physical Chemistry
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 X is the vector of species populations and 𝐊 𝐫 is a 2x2 matrix of rate constants. These rate equations are exactly the rate equations required to describe the kinetics of this reaction system in a chemical kinetics model. SI-I shows that the eigenvalues of 𝐊 𝐫 are identical to those of the CSEs, and that the solution of the master equation, and generation of the rate coefficients through a CSE analysis, provides exactly the quantities required for a chemical kinetics model. The exact relationship between the rate coefficients and the eigenvalues is discussed in the SI-I. The particular case we consider is almost the same as the one considered by BFG. The only differences are that we use a temperature of 1200 K (compared with BFG's 1260 K) and a collision rate with the wall of 1.38 x 104/s (compared with BFG's 2.0 x 104/s.) The essential features of the two cases are the same. Some numerical results of significance from the CSE analysis are given in Table I. In this problem there are 2 CSEs and a multitude of IEREs. The CSEs are well separated from the IEREs (by a factor of about 300, see Table I), so that the fundamental assumption of the method is easily satisfied. The forward and reverse rate constants for the cis-trans isomerization, k21 and k12, satisfy detailed balance (k21/k12 = Keq) to several significant digits. The method also gives dissociation rate constants for both the cis and trans isomers, k31 and k32; again the latter is a well-skipping (or well-mediated) reaction. Examination of the eigenvectors shows that the faster-relaxing of the 2 CSEs primarily relates to equilibration of the populations of the 2 wells, while the slower of the 2 CSEs equilibrates the equilibrated pair of isomers with the fragments. Because the fragments are assumed to be an infinite sink (both here and by BFG), the final equilibrium condition corresponds to all the initial population having reacted to form fragments. In Fig. 1 we compare the populations, Xi(t), as a function of time that come directly from the solution of the master equation with the solution to a system of rate equations that use the rate constants derived above. The two results are virtually identical. That is not an accident. The rate constants are derived to be an “exact” replication of the long-time behavior of the solution to the ME. By “long-time” we mean the behavior after vibrational-rotational energy relaxation has died out. A log-log plot of the populations, shown in Fig. 1b, reveals that the rate-constant description becomes accurate after roughly 2 ms. This is about 4τ, where τ is the time constant corresponding to the slowest-relaxing IERE. This delay time is weakly sensitive to the initial conditions and can be even shorter than 2 ms (see SI-II). SI-II also shows that the conclusion stays the same when starting in the trans well. In the butene-2 reaction system, under these conditions, the extent of reaction under the energy relaxation timescale is negligible (
