J. Phys. Chem. 1995, 99, 13452-13460
13452
Fitting of Pressure-Dependent Kinetic Rate Data by Master Equatioflnverse Laplace Transform Analysis S. H. Robertson, M. J. Pilling," and D. L. Baulch School of Chemistry, University of Leeds, Leeds LS2 9JT, U.K.
N. J. B. Green Department of Chemistry, King's College, London, Strand, London WC2R 2LS, U.K. Received: March 10, 1995; In Final Form: June 5, 199.5@
Techniques for the analysis of pressure-dependent kinetic rate data based on the master equatiodinverse Laplace transform (ME/ILT) method are presented. It is shown that the efficiency of the solution of the ME is substantially enhanced if the ME is replaced by a diffusion equation (DE) and if a reduced matrix, based on a restricted energy range, is used. The speed of the ILT calculation is enhanced by using the fast Fourier transform technique. Case studies for the reactions, K3H7 C3H6 -t H (1) and CH3 CH3 C2H6 (2), are presented. An improved of 87.7 kJ mol-' is obtained from the enthalpy change for reaction 1 calculated from published experimental data (cf. 90.0 obtained by Seakins et al. J Phys. Chem. 1992, 96, 9841). A similar analysis of reaction 2 produces a x2 surface that is very flat and has no well defined minimum because of strong correlation between fitted parameters. Thus, it is difficult to extract a unique parameter set for ky. It is found that the modified Arrhenius form 5.4 x 10-8(T)-'.'o exp((-160 K)/T) cm3 molecule-' s-' gives a satisfactory fit to the available experimental data.
-
1. Introduction The analysis, interpretation, and extraction of physically meaningful parameters from experimental kinetic rate data remains a fundamental problem facing kineticists. The complex dependence of unimolecular rate coefficients on pressure and temperature is particularly challenging. Considerable progress in data analysis has been made in recent years. One of the most widely used techniques is that due to Troe and c o - ~ o r k e r s , l - ~ based on the Lindemann-Hinshelwood model of unimolecular fall-off. The basic Lindemann-Hinshelwood model is retained, and factors correcting for the energy dependence of the microcanonical rate coefficient and weak collision effects are appended. The great advantage of these correcting factors is that they can be represented by relatively simple analytic formulas, an important consideration in the modeling of combustion or atmospheric processes where rate coefficient calculation has to be rapid. A disadvantage of this method is that it is not always possible to extract meaningful parameters from rate data or to assess the validity of theoretical models. Furthermore, in modeling it is often required to extrapolate rate coefficients beyond the region in which they are directly experimentally accessible-to do this requires a sound theoretical foundation; there is no reason for any simple parametrization to be accurate outside of the fitting range. Consequently, there remains a need for less empirical techniques which can be used to analyze data and produce a realistic numerical representation that can then be parametrized using, for example, the T r ~ e ' method -~ to produce an analytical representation for application in realistic reaction schemes. The purpose of this paper is to describe a series of techniques based on the master equation4 (ME). This method has often been seen as complicated and computationally expensive; however, with the dramatic increase in the power @Abstract published in Advance ACS Abstracts, August 15, 1995.
0022-365419512099- 13452$09.00/0
+
-
and abundance of computers and improvements in the technique, ME calculations can become a routine method in data analysis. In section 2 the theoretical background is outlined. The ME is formulated, and diffusion equation (DE) approximations to it are introduced. The advantage of DE approximations is that they can be computed much faster than MEs, which is an important consideration in data analysis. The speed of both the ME and DE calculations of rate coefficients can be greatly enhanced by the use of the reduced matrix approximation (RMA), which is also discussed. Section 3 discusses the inverse Laplace t r a n ~ f o r m -(ET) ~ method for relating microcanonical rate coefficients to experimental limiting high pressure canonical rate coefficients and improvements in the speed of this technique. In sections 4 and 5 two specific examples are considered. The first of these is the unimolecular decomposition of the isopropyl radical,
This reaction plays an important role in combustion and pyrolysis and has been studied by Seakins et a1.* An improved analysis of the data is presented here. The second system is the recombination reaction, CH,
+ CH, - C,H,
(2)
which has received intense ~ c r u t i n y , ~ not - ' ~ only because of its importance in combustion but also because it is ideally suited to theoretical a n a l y s i ~ I ~ and - ' ~ has become a benchmark by which various methods can be assessed.
2. Master Equation There is now a large volume of literature on the application of the ME to gas phase reactions, so only a brief summary is given here. For a unimolecular decomposition reaction the ME
0 1995 American Chemical Society
Fitting of Pressure-Dependent Kinetic Rate Data
J. Phys. Chem., Vol. 99, No. 36, 1995 13453
is written as4*18319
coefficients. The solution of eq 8 is
where @(E,t)is the time-dependent probability density that a given molecule will have rovibrational energy E at time t, w is the collision frequency, P(E(E’) is the transition probability density of the molecule going from energy E’ to energy E on collision, and k(E) is the microcanonical rate coefficient for decomposition. There are a number of restrictions on P(E(E’): it must be normalized because the m o u n t of reactant present is not altered by collisional excitatiodde-excitation processes:
LffiP(EIE’)dE = 1
(4)
Furthermore, in the absence of reaction the system must approach the equilibrium Boltzmann distribution, AE);that is, the system must obey detailed balance,
P(E1E’)flE’) = P(F1E)flE)= W I E ) ME) exp(-PE)lQCa)
(5)
where /3 = l/kT, N(E) is the rovibrational density of states function, and QV)is the corresponding rovibrational partition function of the molecule. Clearly, e(E,t) is dependent on the form of P(E1E’); in common with many previous treatments the exponential down mode120-21 will be adopted here as a physically reasonable model for P(EIE’), P(EIE’) = A@) exp(-a(E’ - E ) )
E‘ L E
(6)
where A is a diagonal matrix of the eigenvalues of M and U is the corresponding eigenvector matrix. The eigenvalues of M are all negative, and for reasonably large threshold energies ( 2 10 000 cm-I) there will be one value, 11, whose magnitude will be orders less than the others. At long times, which correspond to experimental times scales, the decay will be dominated by 11and will be exponential. Equation 10 can then, to a very good approximation, be written as
e@)= UI exp(2,t)
(11)
and the unimolecular rate coefficient kd can be identified with -11. Furthermore, the vector u1 is the steady state normalized population density for the reacting system. Rate coefficient calculation thus becomes a question of determining the eigenvalue of the smallest magnitude for the “spatial operator” in eq 8, Le., the matrix M. For data analysis purposes the calculation of the eigenvalue has to be as efficient as possible. If the grain size is 100 cm-I and the energy range to be considered is 50 000 cm-’, then M will be of order 500. The first thing to note is that it is inefficient to calculate all of the eigenvalues of such a matrix. There are available a number of routines that will calculate only certain eigenvalues. Gilbert and c o - w o r k e r ~ ~have ~ , ~made ~-~~ much use of the Nesbet algorithm for single eigenvalue calculation. Efficiency can be further enhanced by two approximations: the diffusion equation approximation (DE) and the reduced matrix approximation ( M A ) . Each of these is considered in turn. 2.1. Diffusion Equation. A thorough analysis of the use of use of DESto approximate MEs used in unimolecular systems has been given by Green et al.*’ Here we summarize the main results with particular attention being paid to the drift-determined DE. The Kramers-Moya128 expansion of eq 3 gives
where a-l is the mean energy transferred in a deactivating collision, ( A E ) d , and A(E)is an energy-dependent normalization coefficient. This model is very widely used, but currently there is considerable debate as to its suitability, particularly in view of the increasing evidence that so-called supercollisions have significant effects on rate coefficient^.^^,^^ The exponential down model will be retained here because of its ease of use, but it is emphasized that the techniques described below are general and apply to any transition density model. Equation 3 can be solved analytically only for a few simple where systems; for most other cases numerical solution is required. Solution is usually effected by using the graining t e c h n i q ~ e , ‘ ~ , * ~ whereby the energy is partitioned into a series of contiguous intervals or grains, each of which is characterized by a number of states, a mean energy, and a mean microcanonical rate Though exact, eq 12 is clearly difficult to solve, so it is often coefficient. Equation 3 is thus approximated by truncated at the second term in the expansion to give the familiar diffusion equation:
a where m is the number of grains and is chosen such that the population of the mth grain makes a negligible contribution to the rate coefficient. Equation 7 can be written more concisely using matrix notation as Q l d t = %(t)
(8)
where M is given by
M = o(P - I) - K
(9)
the elements of P being the Pv of eq 7, I is the identity matrix, and K is a diagonal matrix containing the microcanonical rate
=
where D(E) (=a2(E)/2) is the infinitesimal variance (or diffusion coefficient) and p(E) (=-al(E)) is the infinitesimal mean (or drift coefficient). The effect of this truncation is that the equilibrium distribution-i.e., the steady state distribution of the system at infinite time in the absence of reaction-is not, in general, the Boltzmann distribution. This can be rectified by imposing the condition that the equilibrium solution is the Boltzmann distributionAE), with the effect that there now exists a direct relationship between p(E) and D(E),
13454 J. Phys. Chem., Vol. 99, No. 36, 1995
Robertson et al. 1, I
1
0
7,
1
Starting with either the definition of p(E) or D Q ,one can then obtain a complementary definition of the other function that would yield the Boltzmann distribution as its equilibrium distribution. In the past it has been common to use eq 15 and the definition of D(E) to determine p(E). Green et al.27 demonstrated that integration of eq 15 to give D(E) in terms of P(E)
-2 -
r.
5
-4
-
-5
-
-3
i
-7 -6 -8 0
'
resulted in a better DE approximation. But what is the advantage of such an analysis? DES are on the whole more amenable to solution-it is occasionally possible to obtain analytic solutions. Even if this is not the case, numerical solution is often simpler. In the present context the numerical analysis of the diffusion operator leads to a tridiagonal matrix representation. The determination of the eigenvalues of such a matrix is considerably faster. How much faster depends on the size of the matrix, but for a matrix of the order of 500 the increase in speed is a factor of 3. 2.2. Reduced Matrix Approximation. The basis of this approximation is the observation that for systems with relatively high thresholds for reaction the distribution far below the threshold tends to a Boltzmann distribution. This behavior is illustrated in Figure 1, which compares a log plot of the probability density of the isopropyl system at equilibrium, the Boltzmann distribution, and a plot of the probability density for a lower pressure system, i.e., the eigenvector U I of eq 11. Figure 2 shows an enlarged plot of the data in Figure 1 for the threshold region showing depletion below the threshold. Figure 1 clearly shows that at lower energies the two distributions are identical, which is physically plausible because it is expected that detailed balance will be the main determinant of relative population densities away from the reaction threshold. This well-established observation can be used to reduce the size of the matrix required for a rate coefficient calculation. Consider the eigenvalue expression for the first eigenvalue 111. Written out in full, this equation is m
m C ~ u ( u l )j m(uI)j - kj(ul)j = ~ l ( u l ) j
(17)
j= I
where (ul)i is the ith component of the eigenvector u1. Summing over i and invoking the normalization condition (eq 4)yield the expected result,
This result can be rewritten as
\
'i 1
m
so00
12000
2
lMM0
m
Ucni'
Figure 1. Ratio of steady state fall-off distribution to the Boltzmann distribution. The reaction threshold is marked by the vertical broken line. 0.05
r
1
\
-0.15 -0.2
\
t
-0.25 9000
11ooo
Elcm"
13W
15WO
Figure 2. As in Figure 1 but showing detail close to the threshold.
than or equal to that of the cth grain, F,. This can be approximated by the corresponding Boltzmann distribution ratio.
Combining eqs 18-20 gives
The accuracy of this approximation improves as c is reduced; however, the size of the corresponding matrix increases. A suitable choice of c is needed to reconcile these features. It depends largely on the system of interest and the width of the graining, but a value of lOkT below the threshold is generally a reasonable estimate. Also, because the distribution declines rapidly in magnitude above the threshold, there is usually little error in truncating the matrix at lOkT above the threshold. In any application these cutoff values must be tested carefully. The size of the matrix required now depends on the temperature, but for lower temperature the matrix will be quite small. As will be shown below, for certain high threshold systems-where the calculation of the rate coefficient is susceptible to numerical e r r o r e q 21 is a convenient way of calculating the forward rate coefficient. 3. Inverse Laplace Transform
where c refers to a grain that is below the threshold. Since kj is zero below the threshold, the numerator of the second quotient in eq 19 is the same as the numerator of eq 18. Comparison of eq 19 with eq 18 suggests that the second quotient in eq 19 is the eigenvalue of a reduced matrix, which has dimension m c. This eigenvalue will be written as 11;. The first quotient in eq 19 is the probability of the molecule having an energy greater
The inverse Laplace transform (ILT) technique may be used to calculate microcanonical rate coefficients from experimental Arrhenius expressions of the high pressure canonical rate coefficient. The idea was fist proposed by Slates who observed that the standard expression that relates the microcanonical and canonical rate coefficients can be written in terms of a Laplace transform:
Fitting of Pressure-Dependent Kinetic Rate Data
J. Phys. Chem., Vol. 99, No. 36, 1995 13455
where ky@) is the high pressure limiting rate coefficient. Equation 23 shows directly how the microcanonical rate coefficient may be obtained if we have k y @ ) . Moreover, the uniqueness theorem of Laplace transforms guarantees that there The is a one-to-one correspondence between k(E) and difficulty with eq 23 is that /cy@) needs to be determined over a wide range of temperatures (strictly an infinite range), which presents experimental difficulties because of the strong dependence of k:@) on p. Errors in k;@) can have dramatic effects on the inferred functional form of k(E). Such difficulties are less pronounced for the association rate coefficient, k y @ ) , which in general is a much weaker function of p. k y @ ) and ky@) are connected via the equilibrium constant K@), so eq 23 can be rewritten as
WE) = Y ' [ QK@) @ky@)I/N(E) )
E >0
Davies et al.' have shown that if the general form of be written as
(24)
rate coefficients for various temperatures and pressures in the fall-off region. This is an example of a system for which the parameters A", n", and E" are already quite well-known, Harris and P i t t ~having ~ ~ measured the association reaction rate coefficient as a function of temperature and close to the high pressure limit. The objectives of Seakins et al. were to obtain estimates of the enthalpy of reaction, @, and (A,?& for He from the kinetic data. The general strategy for this and other parameter determinations is as follows: (i) Use eq 26 to obtain @E). For the present example this involves making an initial guess for do. (ii) Insert the k(E) obtained into the ME/DE and calculate the rate coefficients for the experimental pressure and temperature points for a specific value of ( A E ) d . (iii) Compare the calculated and experimental rate coefficients and adjust the parameters as indicated by the data in accordance with a prescribed criterion. (iv) Repeat this cycle until the best fit criterion is met. As with most fitting procedures the criterion used was that of least squares; that is the best values of the parameters were those that minimized
x2,
kif$) can
for n" > - 1.5, then eq 24 becomes
A@ - x]""+~'~ dx (26) where N e ) is the density of states for the unimolecular reactant N,C) is the convoluted densities of states of the product species, Al$ is the zero-point energy difference between the reactants and products, and r(.)is the gamma function. C' is given by
where Ma and MB are the masses of the dissociation products. The feature that makes the use of eq 26 particularly attractive in data analysis is that it is a convolution and as such can be calculated very rapidly by the use of the fast Fourier transform (FIT)method.30 The use of the technique clearly depends on the parameters A", n", and E", which are not always known and which may be the very parameters that are required from a least squares fit to the fall-off data. In the examples below the use of the ILT in this manner is demonstrated. The main shortcoming of the technique is that it yields no information about the angular momentum dependence of the reaction. Work in progress in this laboratory on two dimensional MEs, which incorporate the (E,J) dependence of the microcanonical rate coefficient, indicates that the errors introduced by J-averaging are not important for the systems studied here.3i Finally, it is noted that this form of ILT produces association and dissociation data that are thermodynamically consistent, based as it is on the equilibrium constant. 4. Isopropyl Radical Decomposition
Seakins et a1.8 studied reaction 1 using laser flash photolysis coupled with photoionization mass spectroscopy and obtained
x2 is thus a function of M0and (A@*, and the best fit corresponds to its global minimum. The nature of the calculation of the rate coefficients means that x2 is a nonlinear function of the parameters, so a nonlinear least squares method is required, for example, the Levenberg-Marquadt algorithm,30 in which the gradient vector is followed to the minimum. The difficulty with this method for the present system is that analytic derivatives are not available and have to be calculated using finite difference which significantly increases the processor time. Further difficulties with this method arise if the nature of the surface is such that either it has local minima in which a search might get trapped or the surface is flat so that the derivatives obtained from finite difference have large errors. As will be shown shortly, this latter difficulty was encountered for the present system and the method was abandoned in favor a simple grid search of parameter space. In performing a grid search it is necessary to calculate x2 for each grid point. The number of grid points for such a search must, inevitably, be large, so fast calculation of x2 is essential. For this reason both the DE and RMA approximations were used. Before the search was undertaken, a check on the accuracy of the above approximations was performed. The results are shown in Table 1. (The values of molecular parameters used are given in Appendix 1.) Values for @ and ( A @ d were fixed at 12 250 and 210 cm-I, respectively. The energy range considered was 0-30 000 cm-I, which was split into 300 grains of lOO-cm-' width. Full ME solutions for 18 of the experimental conditions reported by Seakins et a1.* were calculated, and these are shown in the third column of Table 1 . The corresponding RMA approximation (with a lower boundary set at lOkT below the threshold and the upper boundary set at 12kT above) is shown in the fourth column. The respective CPU times for each set of 18 rate coefficients were 600 and 48 s on a Sun SPARCstation 10. The DE values and the corresponding RMA values are shown in the fifth and sixth columns, respectively, of Table 1 . The respective CPU times for these results are 24 and 5 s. The agreement between the various approximations and the full ME solution is satisfactory and well within experimental error.
Robertson et al.
13456 J. Phys. Chem., Vol. 99, No. 36, 1995 900
YnaK I
600
s
500
-
400
-
810 K
8W K
300 -
-93K-80 K
100 -
200
770 K 760 K 750
0 ' 0
"i I
5
IO
20 [He]/ld6molec/cm 15
25
30
35
Figure 4. Plot of isopropyl dissociation data and calculated fall-off curves.
12000
12500
12250
A% /cm" Figure 3.
x2 surface for the isopropyl decomposition.
Contours are
at intervals of 1.5 x lo4 s - ~ ,
[He] x molecule/cm'
ME
760.0 760.0 760.0 760.0 760.0 760.0 800.0 800.0 800.0 800.0 800.0 800.0 800.0 830.0 830.0 830.0 830.0 830.0
6.02 9.02 12.00 18.10 24.00 30.00 3.00 6.01 9.00 12.10 18.00 24.00 30.00 3.00 6.00 9.00 12.00 18.00
61.63 75.98 87.46 105.99 120.06 131.92 105.98 159.80 199.98 234.07 286.33 328.99 364.86 195.52 299.75 379.59 445.91 554.24
5. Methyl Radical Recombination The recombination of methyl radicals occurs via the scheme
TABLE 1: Comparison of ME, RMA/ME, DE, and RMA/DE Approximations for Isopropyl Decomposition temdK
rigorous analysis of the errors of these estimates is currently being investigated.
WS-' ME/RMA DE 61.98 76.34 87.82 106.35 120.40 132.25 107.16 161.19 201.46 235.59 287.90 330.55 366.41 198.38 303.22 383.36 449.88 558.40
60.39 14.49 85.80 104.06 117.94 129.67 103.85 156.68 196.17 229.71 281.17 323.23 358.63 191.73 294.07 372.54 437.78 544.41
CH3 + CH3 === C 2 H t
(29)
DE/RMA 60.74 74.86 86.18 104.45 118.34 130.05 105.03 158.07 197.66 231.26 282.79 324.89 360.31 194.59 297.53 376.32 441.76 548.65
For the grid search, the DE/RMA method was used and x2 values for 500 grid points were calculated. The full experimental data set of 71 points was used. Initially, grid points close to the optimal values obtained by Seakins et aL8 were explored, but it became apparent that this minimum was located in a region of parameter space which consists of a long trench. The x2 surface obtained from a wider search is shown in Figure 3; the contour lines are at intervals of 15 000 s - ~ . This surface shows a trench with a flat bottom, and, as stated above, this is the sort of minimum that Levenberg-Marquadt algorithms have difficulties in searching, because the minimization path tends to bounce between the two steep walls of the trench. It is also apparent that the best value of (AQd is very strongly correlated to the value of r\H",. Even so, good definition of both parameters, and particularly of @, is obtained. Finally, it is clear that although previous estimates for these parameters were satisfactory, better estimates are do = 12 275 cm-' and ( a d = 230 cm-'. Figure 4 shows a plot of the experimental data together with calculated fall-off curves. From this value of the reaction enthalpy change an estimate of the enthalpy of formation of the isopropyl radical of 87.7 kJ/mol is obtained (cf. 90.0 kJ/mol obtained by Slagle et al.). A method for the
where M is the bath gas. This is a specific example of the more general recombination,
considered by Smith et a1.26who demonstrated that the overall recombination rate coefficient was related to the dissociation rate coefficient by the equilibrium constant regardless of the pressure of the system. Calculation of the dissociation rate coefficient represents one possible route to the rate coefficient for recombination. Another route is to modify the rate ME given in eq 3 to take account of the reactive input:
The term Ri(t) represents reactive gain due to bimolecular association and is constrained by detailed balance at equilibrium. The equilibrium limit is found by substitutingj for gi(t) in eq 31:
Summing over i gives (33) where k;P is the bimolecular recombination rate coefficient in the high pressure limit. Combining eqs 32 and 33 yields
J. Phys. Chem., Vol. 99, No. 36, 1995 13457
Fitting of Pressure-Dependent Kinetic Rate Data where 4; partitions the rate of formation between the grains of the adduct. If A and B remain in thermal equilibrium with the bath gas, the same relationship holds throughout the course of reaction. (35) where R'(t) represents the overall rate of association. R'(t) is nonlinear in t, so relaxation is nonexponential. When an association has taken place, there are two possible fates for the molecule; either it will be deactivated or it will redissociate before deactivation. The observed association rate coefficient measures the overall flux from reactants through to deactivated product. This rate can be modeled using a ME method given a suitable operational definition of deactivation. In the calculations presented here, the molecule is considered to be deactivated if the energy reaches at least lOkT below threshold. The deactivation rate is therefore modeled as the flux into a cemetery state placed at this energy. The solution for this irreversible system is characterized by a fast transient period followed by a steady state regime in which the shape of the population distribution remains constant. The steady state solution can be obtained from the ME, eq 31, by setting the time derivative and all reactivation rates from the cemetery state to zero. This effectively reduces the size of the matrix by removing all states below 10kT. In this form matrix inversion can be regarded as the recombination equivalent of the RMA. The steady state population distribution g can then be obtained by inversion of the truncated matrix.
This is a specific example of the more general solution to eq 36 obtained by Smith et a1.26 The observed rate of association, Rob, is given by the difference between the input rate R'(t) and the unimolecular dissociation rate:
= R'(t) - C k g ,
(37)
where k, is the observed rate coefficient for association. Substitution for gi from eq 36 and setting [A] = [B] = 1 yields
k, = kr(1
+ Xki(M-'I$)J i
This approximation was tested for the methyl recombination system, by comparing rate coefficients obtained from using eq 38 with those from the eigenvalue analysis for the reverse dissociation reaction coupled with the equilibrium constant. Microcanonical rate coefficients, used in both the inversion and eigenvalue analyses, were obtained using ILT. (The molecular parameters used in this study are gathered together in Appendix 2.) The dissociation rates were found using the M A approximation described in section 2. This system illustrates one of the benefits of this method, because if the full matrix is used, unrealistic values of the rate coefficients are obtained for low temperatures. This is a consequence of the high threshold energy for ethane dissociation, which means that the population of states at and above the threshold is very small leading to large numerical error. The results are shown in the third column of Table 2. In general the agreement is satisfactory except at low pressures,
TABLE 2: Comparison of Dissociation and Association Rate Coefficients for CH3 CH3 at T = 300 K
+
Lg([Hel/ (molecule ~ m - ~ ) ) 4.0 5 .O 6.0 7 .O 8 .O 9.0 10.0 11.0 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23.0
Lg(10"k,l(cm3 molecule-' s-l)) kdKe, ea 38 ea 39 -7
-9.8842 -8.8842 -7.8842 -6.8842 -5.8842 -4.8845 -3.8870 -2.9091 -2.0044 - 1.1348 -0.3583 0.2433 0.6075 0.7542 0.7875 0.7919 0.7923 0.7924 0.7923 0.7926
-2.0152 -2.0152 -2.0152 -2.0152 -2.0151 -2.0 146 -2.0095 - 1.9635 - 1.7105 - 1.0834 -0.3506 0.2450 0.6083 0.7550 0.7882 0.7925 0.7930 0.7930 0.7930 0.7930
-9.8886 -8.8886 -7.8886 -6.8886 -5.8887 -4.8889 -3.8914 -2.9134 -2.0079 - 1.1375 -0.3601 0.2426 0.6073 0.7542 0.7875 0.7919 0.7923 0.7924 0.7924 0.7924
where the rate coefficient appears to converge to an asymptote, rather than continue to fall-off as expected. The origin of this error is the truncation of the matrix: back-reaction is slightly underestimated, and the difference appears at low pressures where back-reaction becomes increasingly important. This error can be corrected by increasing the size of the matrix used, at the expense of CPU time. A more satisfactory approach is to recognize that in the steady state the overall association rate is the same as the total rate into the cemetery state; hence, the rate coefficient is given by
where Pc,i is the transition probability from state i to the cemetery state. The advantage of this approach is that as i increases, Pc,i decreases rapidly so that the sum in eq 39 converges rapidly. This approach was also tested, and the results are shown in the fourth column of Table 2. Note the improved agreement at lower pressures with the rate coefficient calculated from the dissociation rate coefficient. The relative amounts of CPU time by the eigenvalue and inversion methods depends on the size of the matrix involved, which in turn depends on the temperature. For 1000 K the inversion method is the faster by 33%. At lower temperatures the methods tend to converge in CPU time requirement. Having obtained satisfactory agreement between forward and reverse rate coefficients, the method described by eq 39 was then used to analyze the experimental data for methyl recombination. Again a grid search approach was used. The data to be fitted were taken from ref 10-14, a total of 210 experimental points. The experimental methods used are summarized in Table 3. Recent work by Hessler and co-workersI4 suggests that the high temperature rate coefficients measured by Slagle et al.'I are in error because of an incorrect estimate of the absorption cross section, u, of the methyl radical. The directly measured quantity in the absorption experiment reported by Slagle et al. is klu. The u values needed in the 537-906 K range were obtained by extrapolation of direct measurements for the range 296 IT < 537 K and were based on a model of the absorption
Robertson et al.
13458 J. Phys. Chem., Vol. 99, No. 36, 1995
TABLE 3: Experimental Methods Used To Investigate Methyl Recombination ref
10 11 12 13 14
method laser flash photolysis/ UV absorption laser flash photolysis/ photoionization mass spectroscopy and laser flash photolysis/UV absorption shock tube/UV absorption laser flash photolysis/ UV absorption and discharge flow/ mass spectroscopy shock tubenaser absorption
process.I0 Direct measurements at higher temperatures using shock tubes indicate that the extrapolations were in error. These rate constants were corrected before fitting, using the interpolated values of u obtained by Hessler and c o - ~ o r k e r s . ' ~ The parameters of key interest in the present case are the Arrhenius parameters, A", n", and E". A fit of all three parameters simultaneously is possible but computationally expensive, and because of this a constraint on the parameters was introduced. At 300 K, experimental data are available close to the high pressure limit and a value for the high pressure rate coefficient at 300 K has been evaluated by Baulch et al.33to be 6.0 x lo-'' cm3 molecule-' s-l. The Arrhenius parameters n" and E" were floated, and the parameter A" was adjusted so as to give this value. The energy transfer parameter (AE)d was also fixed at 270 cm-l. The results are shown in Figure 5,
TABLE 4: Arrhenius Curves for Methyl Recombination ~~
~
ref
A"/(cm3 molecule-' s-')
n-
E"/K
17 33 34
1.5 x 10-7 6.0 x 10-7 3.53 x 10-8
-1.18 0.0 -0.973
329 0 312
,',
Figure 5.
x2 surface for the methyl recombination:
x2
which shows a plot of the surface of the fit. This surface was obtained from 63 x2 points which took approximately 2 h to generate. It is clear from Figure 5 that the surface is characterized by a broad flat minimum similar to that found for isopropyl, but in this case the minimum encompasses a great deal more of the parameter space. As a consequence it is not possible to assign with confidence a given set of parameters for the high pressure rate coefficient because no one point is statistically more significant than another. There are two possible explanations for this: the data may be inadequate to fit these parameters, or there may be inconsistencies in the data that prevent a definitive fit from being obtained. Arrhenius parameters determined by Wagner and Wardlaw,I7 Baulch et al.,33and Golden et al.34are collected in Table 4 and plotted in Figure 6. Also shown are two curves obtained from the above fitting exercise. The parameter sets (A" = 4.48 x lo-' cm3 molecule-' s-l, n" = - 1.40, F I R = 280 K) and (A" = 0.54 x cm3 molecule-' s-', n- = -1.10, E"/R = 160 K) are taken form opposite ends of the trench in Figure 5. As can be seen in both cases, the temperature dependence is stronger than those from previous estimates.
units of lo2*cm6 molecule-2 s - ~ .
x2
Fitting of Pressure-Dependent Kinetic Rate Data
J. Phys. Chem., Vol. 99, No. 36, 1995 13459 Acknowledgment. An EPSRC Research Fellowship (S.H.R.) is gratefully acknowledged. Thanks are also owed to Dr. Jan Hessler for helpful discussions and for access to experimental results prior to publication.
-i. h
Appendix 1 Data for isopropyl decomposition. (All values given in cm-I.) Vibrational frequencies of C3H7. Ref. 34 Model1 A Model I1 x
3052.0 2968.2 2968.2 2967.7 2920.0 2887.0 2850.0 1468.0 1464.0 1462.0 1440.0 1378.0 1338.0 1292.0 1200.0 1165.0 1053.8 921.7 879.0 748.1 369.2 364.0 External rotational constants of C3H7, 1.528 0.2645 0.2645 0
0.6 -
O
O
0
0
O
@ T
A
-
+
+
y
o
Internal rotational constants of C3H7.
T
h
0.4
1
0.2
1
+ +
5.36053
i
Barrier height to intemal rotation of C3H7. 255.0
0 1
-0.4 -o.2 -0.6
I '
0
Vibrational frequencies of C3H6. 296 K 474 K 810 K 906 K 0.5
1
1.5
2 2.5 Log ([AT]10'6/moleccm-3)
3
0
t 0
1
X
I 3.5
Figure 7. Plot of methyl recombination data and calculated fall-off
3091.0 3022.0 2991.0 2973.0 2952.8 2931.0 1652.8 1458.5 1442.5 1414.01378.0 1298.0 1177.5 1044.7 990.0 934.5 919.0 912.0 575.2 428.0 External rotational constants of C3H6.
curves. A selection of rate coefficients and the calculated fall-off curves are plotted in Figure 7 for the parameter set (A" = 5.4 x lo-* cm3 molecules-' s-l, rim = -1.10, E"/R = 160 K). The other parameter set above gives a similar fit. 6. Summary
It has been shown that the analysis of experimental kinetic rate data can be addressed using ME techniques directly, thus avoiding the use of parametrized fall-off models. Two approximations are employed to facilitate the rapid calculation of canonical rate coefficients: the diffusion equation approximation replaces the master equation collision operator with a diffusion operator which can be represented in numerical form by a tridiagonal matrix. The reduced matrix approximation originates from the observation that only states in the immediate vicinity of the threshold are significantly depleted from the Boltzmann distribution. For recombination reactions, a similar reduction in the size of the matrix can be achieved by employing a cemetery state. It was shown that the DE approximation reproduces the results of the ME results to within 2% and that the RMA further accelerates computation with only a minor reduction in accuracy. A further requirement for the ME is a set of microcanonical rate coefficients k(E) which can be obtained rapidly by the ILT method in conjunction with the fast Fourier transform technique. The combination of these techniques allows the rapid calculation of points on the x2 surface which can then be analyzed to locate the best fit parameters. In future work it is intended to combine these techniques with automatic methods for determining minima and a rigorous analysis of errors.
1.544 0.2905 0.2905 Internal rotational constant of C3H6. 5.36053 Barrier height to internal rotation of C3H6. 699.0 Appendix 2 Data for the methyl radical recombination. (All values given in cm-I.) Vibrational frequencies for C&. 2985 2985 2969 2969 2954 2896 1469 1469 1468 14681388 137911901190995 822822 Rotational constants for C2H6. 2.671 0.6625 0.6625 Internal rotational constant for C2H6. 5.342 Barrier height to internal rotation for C2H6. 1012.32 Vibrational frequencies for CH3. 31623162304413961396606
13460 J. Phys. Chem., Vol. 99, No. 36, 1995 Rotational constants for CH4. 4.789 9.578 9.578 References and Notes (1) Troe, J. J . Phys. Chem. 1979, 83, 114. (2) Troe, J. Ber. Bunsen-Ges. Phys. Chem. 1983, 87, 161. (3) Gilbert, R. G.; Luther, K.; Troe, J. Ber. Bunsen-Ges. Phys. Chem. 1983, 87, 169. (4) Oppenheim, I.; Shuler, K. E.; Weiss, G. H. Stochastic Processes in Chemical Physics: The Master Equation; MIT: Cambridge, MA, 1977. (5) Slater, N. B. Proc. Leeds Philos. Lit. SOC.1955, 6, 259. (6) Pritchard, H. 0. Quantum Theory of Unimolecular Reactions; Cambridge University Press, Cambridge, 1984. (7) Davies, J. W.; Green, N. J. B.; Pilling, M. J. Chem. Phys. Lett. 1986, 126, 373. (8) Seakins, P. W.; Robertson, S. H.; Pilling, M. J.; Slagle, I. R.; Gmurczyk, G. W.; Bencsura, A.; Gutman, D.; Tsang, W. J . Phys. Chem. 1993, 97, 4450. (9) Glaenzer, K.; Quack, M.; Troe, J. Chem. Phys. Lett. 1979,39, 304. (10) Macpherson, M. T.; Pilling, M. J.; Smith, M. J. C. J. Phys. Chem. 1985,89, 2268. (11) Slagle, I. R.; Gutman, D.; Davies, J. W.; Pilling, M. J. J . Phys. Chem. 1988, 92, 2455. (12) Hwang, S. M.; Wagner, H. G.;Wolff, T. Twenty-Third Symposium (International)on Combustion; The Combustion Institute: Pittsburgh, PA, 1990; p 99. (13) Walter, D.; Grotheer, H.-H.; Davies, J. W.; Pilling, M. J.; Wagner, A. F. Twenty-Third Symposium (International) on Combustion; The Combustion Institute: Pittsburgh, PA, 1990, p 107. (14) Hessler, J. P. Private communication. (15) Wardlaw, D. M.; Marcus, R. A. J . Chem. Phys. 1985, 83, 3462.
Robertson et al. (16) Wardlaw, D. M.; Marcus, R. A. J . Phys. Chem. 1986, 90, 5383. (17) Wagner, A. F.; Wardlaw, D. M. J . Phys. Chem. 1988, 92, 2462. (18) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell: Oxford, U.K., 1990. (19) Forst, W. Theory of Unimolecular Reactions; Academic: New York, 1973. (20) Penner, A. P.; Forst, W. J. J . Chem. Phys. 1977, 67, 5296. (21) Murrell, J. N.; Bosanac, S. D. Introduction to the Theory ofAtomic and Molecular Collisions; Wiley: New York, 1989. (22) Morgulis, I. N.; Sapers, S. S.; Steele, C.; Oref, I. J . Chem. Phys. 1989, 90, 923. (23) Bemshtein, V.; Oref, I. J . Phys. Chem. 1993, 97, 12811. (24) Hanning-Lee, M. A.; Green, N. J. B.; Pilling, M. J.; Robertson, S. H. J. Phys. Chem. 1993, 97, 860. (25) Smith, S. C.; McEwan, M. J.; Gilbert, R. G. J . Chem. Phys. 1989, 90, 1630. (26) Smith, S. C.; McEwan, M. J.; Gilbert, R. G. J. Chem. Phys. 1989, 90, 4265. (27) Green, N. J. B.; Robertson, S. H.; Pilling, M. J. J . Chem. Phys. 1994, 100, 5259. (28) Moyal, J. E. J. R. Stat. SOC.B 1949, 11, 145. (29) Widder, D. V. The Laplace Transform; Princeton University Press: Princeton, NJ, 1946. Chandler, D. An. Introduction to Modern Statistical Mechanics; Oxford, University Press: Oxford, U.K., 1987. (30) Press, W. H.; Flannery, B. P.; Teukolsky, S. A,; Vetterling, W. T. Numerical Recipes; Cambridge University Press: Cambridge, U.K., 1988. (31) Robertson, S. H.; Green, N. J. B.; Pilling, M. J.; Baulch, D. L. To be published. (32) Harris, G.W.; Pitts, J. N., Jr. J . Chem. Phys. 1982, 77, 3994. (33) Baulch, D. L.; Cobos, C. J.; Cox, R. A,; Esser, C.; Franck, P.; Just, Th.; Kerr, J. A.; Pilling, M. J.; Troe, J.; Walker, R. W.; Wamatz, J. J . Phys. Chem. Re5 Data 1992, 21, 411. (34) Stewart, P. H.; Larson, C. W.; Golden, D. M. Combust. Flame 1989, 75, 25. JP9506954
