984
J. Phys. Chem. 1996, 100, 984-992
Recombination of Methyl Radicals. 2. Global Fits of the Rate Coefficient† Jan P. Hessler* Gas Phase Chemical Dynamics Group, Chemistry DiVision, Argonne National Laboratory, 9700 South Cass AVenue, Argonne, Illinois 60439
Paul J. Ogren Department of Chemistry, Earlham College, Richmond, Indiana 47374 ReceiVed: May 2, 1995; In Final Form: August 21, 1995X
The temperature- and pressure-dependent behavior of the cross section for optical absorption by the methyl radical is carefully considered, so we may define a criterion for selecting and correcting measurements of the rate coefficient for the recombination of methyl radicals, CH3 + CH3 f C2H6. The low-temperature data of Slagle et al., Hippler et al., and Walter et al. and the high-temperature data of Gla¨nzer et al., Hwang et al., and our latest results (previous paper in this issue) are used to define a data set which contains 217 points. Subsets of isothermal data show that the temperature dependence of the high-pressure rate coefficient may be described by the simple exponential function A∞ exp{-T(K)/T∞}. Four different formulations for the pressuredependent behavior in the falloff region are used for the global fits: (1) the asymmetric Lorentzian broadening function of Gilbert et al.; (2) the Gaussian broadening function of Wang and Frenklach; (3) the empirical “a equation” introduced by Gardiner; (4) the extension of Lindemann’s expression suggested by Oref. All formulations reproduce the data, but Oref’s “J equation” produces the least correlation between the best-fit parameters, the least uncertainty in these parameters, and the smallest uncertainty in the predictions. These results are k∞(T) ) 8.78 × 10-11 exp{-T(K)/723} cm3 s-1, k0(T) ) 9.04 × 10-27 cm6 s-1, and Jexp(T) ) {exp[T(K)/268] - 1}2.
I. Introduction The recombination of methyl radicals proceeds on a barrierless potential energy surface and represents a prime example of a simple fusion reaction. In a previous paper1 we reported measurements of the rate coefficient for recombination of the methyl radical to form ethane over the temperature range 11751750 K and the pressure range 114-230 kPa (1.13-2.27 atm). In this paper, we demonstrate that the temperature and pressure dependence of the cross section for optical absorption must be determined for each different experimental apparatus. We use the detailed behavior of the cross section as a criterion to select a set of data. With this set, four different empirical formulas for the pressure-dependent rate coefficient are used to extract the temperature-dependent behavior of the limiting rate coefficients and produce global fits. An important new aspect of this work is the use of Monte Carlo simulations to estimate the confidence envelopes of the predictions and the confidence limits of the best-fit parameters. The empirical formulations for the behavior in the falloff region suggested by Gilbert et al.,2 Wang and Frenklach,3 and Gardiner4 are discussed. In addition, the extension of Lindemann’s model in terms of a “J equation” introduced by Oref5 is discussed. All of these formulations adequately reproduce the experimental data and show that the temperature dependence of the high-pressure rate coefficient follows the simple exponential equation A∞ exp{T(K)/T∞}. Unfortunately, different formulas for the behavior in the falloff region give statistically different values for the parameters A∞ and T∞. Oref’s “J equation” produces the least † Work performed under the auspices of the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences, under Contract No. W-31-109-Eng-38. * To whom correspondence should be addressed. X Abstract published in AdVance ACS Abstracts, December 15, 1995.
0022-3654/96/20100-0984$12.00/0
correlation between the best-fit parameters, the least uncertainty in these parameters, and the smallest uncertainty in the predictions. Because of the complexity of the potential energy surface for methyl-methyl recombination, a global potential energy surface based upon ab initio calculations does not exist. The most recent and complete ab initio calculations were reported by Robertson et al.,6 who produced a “semiglobal” surface which is based upon “ab initio calculations for an extensive, but judiciously restricted, set of geometries”. At the present time, dynamical calculations based upon this improved surface are being performed.7 The earliest calculation of the temperature dependence of the high-pressure rate coefficient was performed by Gorin.8 His ideas have been extended and improved upon by Gilbert and co-workers.9 An alternate approach was suggested in 1977 when Quack and Troe10 extended the statistical adiabatic channel model11 to this reaction. The latest calculations of these two approaches predict there is a slight positive temperature dependence of the high-pressure rate coefficient. Wagner and Wardlaw12 have used variational transition state theory and obtain a significant negative temperature dependence in the high-pressure rate coefficient. The most recently reported dynamical calculations are by Walter et al.,13 who used a parametrized potential energy surface. The dynamics were calculated by a “weak-coupling corrected variational activated complex RRKM calculation based upon a FTST description of the transition state”. In these calculations two parameters, R which controls the presumed exponential switching of structures and conserved frequencies along the reactions path, and 〈∆E〉tot which is the average energy transferred in an Ar-C2H6* collision, were adjusted to obtain agreement with experimental data below 1000 K. This adjustment produced a substantial positive temperature dependence of 〈∆E〉tot and, more importantly, gave a significant negative temperature dependence to © 1996 American Chemical Society
Recombination of Methyl Radicals. 2
J. Phys. Chem., Vol. 100, No. 3, 1996 985
the high-pressure rate coefficient. The functional form of this temperature dependence agrees with the global fits presented here. However, its magnitude is sufficiently different from theoretical calculations that a new evaluation of the dynamics is warranted. In section II of this paper, we discuss the cross section for optical absorption, review the experimental situation, and use the cross section measurements to define a set of acceptable experimental data. In section III, we describe the fitting procedure and Monte Carlo simulations used to produce global fits of the data. In section IV, we use selected subsets of isothermal data to determine the temperature-dependent behavior of the limiting rate coefficients and present the results of four different global fits. Section V is devoted to a discussion of the implications of these results, and the conclusions are presented in section VI. II. Cross Section for Optical Absorption Most of the recent experiments designed to measure the recombinational kinetics of the methyl radical use optical absorption to measure the time-dependent absorption coefficient, R(t). If the rate of formation of the methyl radical is much faster than the rate of recombination, this coefficient may be written
R(t) )
R0 1 + 2krec[CH3]0t
(1a)
where the inital absorption is R0, the rate coefficient for recombination is krec, and the maximum concentration of the methyl radical is [CH3]0. The cross section for optical absorption, σ, is given by R0/[CH3]0. In high-temperature experiments performed with shock tubes, pyrolysis of a precursor with a known yield is generally used to produce the methyl radicals. Therefore, since the maximum concentration of the methyl radical may be calculated from the initial conditions of the experiment, both the absorption cross section and the rate coeficient for recombination may be determined in a single experiment. If the rate of dissociation of the precursor is comparable to the rate of recombination eq 1a becomes more complex,14,15 but the above conclusion is still valid. If an alternate approach is used to produce the methyl radicals, for example a precursor with an unknown yield or simple photolysis, then it may not be possible to calculate the initial methyl concentration and eq 1a must be rewritten to give
krec 1 1 ) +2 t σ R(t) R0
(1b)
In experiments of this type, the cross section must be obtained from other measurements. The correlation between the rate coefficient for radical recombination and the cross section is apparent in eq 1b. Perhaps less apparent, but just as important, is the correlation in eq 1a.16 Since most of the recombinational experiments have used the second approach, we must carefully consider the behavior of the cross section for optical absorption. Unfortunately, it is not possible to predict how the cross section will behave as the pressure is increased.17 Line narrowing has been observed in the hydroxyl radical,18 and significant frequency shifts have been observed in NO.19 Collisional broadening of hydrogen has also been observed in a flat-flame H2-air burner.20 Optical absorption by the methyl radical near 216 nm has an additional complication because the absorption band is composed of many overlapping rovibronic transitions. Therefore, the cross section should be measured
as a function of temperature, pressure, and gas composition for each set of experimental conditions. Furthermore, one must exercise caution in using cross sections determined in one laboratory to quantitatively correct or infer something about measurements performed in another laboratory. The preceding statements are supported by the following experimental observations. Macpherson et al.21 have measured the absorption cross section of methyl at 216.36 nm with a bandwidth of 0.6 nm between 3.5 and 53.3 kPa (0.035 and 0.53 atm) at 298 K. They obtain a pressure-independent cross section of 4010 ( 260 pm2 which is in agreement with that determined by end-product analysis, 4120 ( 240 pm2. On the other hand, Hippler et al.22 performed similar experiments (216.2 nm with a spectral resolution of 0.1 nm) at room temperature over the pressure range 0.10-21 MPa (1-207 atm). From 0.1 to 3.1 MPa (130 atm) the cross section is independent of pressure at 4200 pm2 but decreases to 1836 pm2 at 21 MPa (207 atm). Clearly, above 3.1 MPa (30 atm) a pressure-dependent cross section is observed. At much higher temperatures Hwang et al.23 measured methyl absorption at 216.5 nm with a spectral resolution of 0.8 nm and, also, observed a pressure-dependent cross section. At 1500 K as the pressure was increased from 1.03 MPa (10.2 atm) the cross section increased from 333 to 815 pm2 near 10 MPa (99 atm) and then decreased to 640 pm2 at 14 MPa (138 atm). Now that we have established the general behavior and role of cross section measurements, we proceed to select data for the global fits. The low temperature data of Slagle et al.24 may be used because they have experimentally verified the pressureindependent behavior21 of the cross section over the range of their 296-577 K experiments. We include the low-pressure data of Walter et al.13 These data were obtained with a mass spectrometer and significantly extend the low-pressure measurements at room temperature and 407 K. To these results we add the high-pressure room-temperature data of Hippler et al.22 from 0.1 to 1.0 MPa (1-10 atm). At higher pressures diffusioncontrolled kinetics begins to decrease the recombinational rate. We do not use the kinetic data of Macpherson et al.21 because they only report limiting high- and low-pressure rate coefficients, not original data. We include the pioneering work of Gla¨nzer et al.25,26 Although they do not present measurements for the cross section over the pressure range of their experiments, they did not detect any pressure-dependent changes in the cross section.27 Also, we combine our data with the high-pressure data of Hwang et al.23 Initially, we included the data of Zaslonko and Smirnov,28 but their results above atmospheric pressure deviate by much more than 3 standard deviations from the other results. We assign a weight to each point of 1/σ2i where σi is the standard deviation of the ith point. Standard deviations were not reported by Hippler et al. and Hwang et al.; on the basis of the scatter of their data, we assign standard deviations of 10 and 20%, respectively, to each of their points. To utilize experimental results obtained with optical techniques at 700, 810, and 905 K reported by Slagle et al.,24 we must correct them. Recall in Figure 4 of our previous paper1 the cross sections at these temperatures are anomalously high. These cross sections were not determined experimentally but were based upon a quadratic extrapolation of the lower temperature cross sections. This behavior was based on unpublished calculations.24 Now that we have several different sets of high-temperature data, we feel a better approximation to the temperature-dependent behavior of the cross section is σ(T) ) A exp{-T(K)/B}. To correct their data, we use this expression to extrapolate their lower-temperature cross section measurements. We obtain σ(pm2) ) 7683 exp{-T(K)/515.5},
986 J. Phys. Chem., Vol. 100, No. 3, 1996
Hessler and Ogren
which gives σ(700 K) ≈ 1976 ( 119 pm2, σ(810 K) ≈ 1596 ( 144 pm2, and σ(906 K) ≈ 1325 ( 155 pm2 where the uncertainties are at the 1σ (68.3%) confidence level. Since the rate coefficients determined by the optical techniques were extracted by measuring the slope, k/σ, of inverse absorption profiles vs time, eq 1b, we may correct the published rate coefficients by multiplying by the ratio of the cross sections. This implies that to correct the 700, 810, and 906 K data we multiply by 0.919, 0.818, and 0.694, respectively. With these corrections we have a data set which contains 217 points and spans the ranges 296 e T(K) e 1742 and 0.020 e P(kPa) e 23 790. This set will be used to perform the global fits. III. Fitting Procedure and Monte Carlo Simulations Recently, Pilling29 has reviewed the various methods of fitting kinetic data to obtain limiting rate coefficients. In this work, we consider different empirical formulations for the pressuredependent behavior in the falloff region and determine the limiting rate coefficients, k∞(T) and k0(T). The LevenbergMarquardt method30 of nonlinear least-squares analysis is used to determine the best-fit values of the independent parameters. Once a best-fit set of parameters has been determined, the best way to estimate confidence limits for multiparameter models, such as those used to describe kinetic rate coefficients, is to perform a Monte Carlo simulation.31 In our case, the Monte Carlo simulation is performed by creating a “synthetic” set of data {kiMC, i ) 1, N} from the original data. Each synthetic datum is generated by calculating a random number R from a Gaussian distribution of random numbers32 and then applying the equation
kiMC ) kie + Rσi
(2)
The ith measured rate coefficient is kie, and its uncertainty is σi. A different random number is calculated for each datum of the set. Note, in some discussions the best-fit values of the independent parameters are used in the analytic functions to generate the synthetic set of data. Since the functions we will be using are only approximations, this approach could introduce significant systematic differences and could destroy or distort important information. Therefore, it is crucial that the experimental data be used to generate the “synthetic” set of data. After each synthetic set of data has been calculated, we analyze it to determine another set of best-fit independent parameters. With these parameters we then calculate the increase in χ2, ∆χ2, with respect to the original, experimental, set of data and save both ∆χ2 and the new set of independent parameters. Depending upon the number of fit parameters, this simulation is repeated up to 4000 times. Once these simulations have been completed, the values of ∆χ2 are arranged in ascending order and the parameters associated with each value of ∆χ2 are retained. The confidence limits of the independent parameters are given by a volume element in the multidimensional space of the independent parameters that contains a given fraction of the distributed results.33,34 For example, the 68.3% confidence limit of a parameter is given by the extreme values of the parameter associated with the lower 68.3% of the values of ∆χ2. The confidence limits of the predictions are generally described in terms of a confidence enVelope35 which is calculated by determining the upper and lower bounds of the prediction from the set of independent parameters within a given volume element. To illustrate these concepts, we fit the 18 points of the isothermal recombinational data, krec(P), at 577 K to the expression
Pr krec(P) ) k∞ F(Pr,T) 1 + Pr
(3)
where the high-pressure rate coefficient is k∞; the reduced pressure is Pr, which is defined as k0P/k∞RT where k0 is the low-pressure rate coefficient. P, R, and T have their usual meaning, and the broadening factor is F(Pr,T). It may be written
F(Pr,T) ) [Fcent(T)]X(Pr,T)
(4)
where Fcent(T) is the centering parameter and X(Pr,T) is the broadening function. There are several empirical expressions for the broadening function.36 Perhaps the most popular is the asymmetric Lorentzian function suggested by Gilbert et al.:2
{ [
X(Pr,T) ) 1 +
log{Pr} + c
]}
2 -1
N - d(log{Pr} + c)
(5)
where c ) -0.4 - 0.67 log{Fcent}, N ) 0.75 - 1.27 log{Fcent}, d ) 0.14. With this broadening function we fit the set of isothermal data to determine the best-fit values of three independent parameters: k∞, k0, and Fcent. For 18 points χ2 is 6.67. The reduced χ2, which is defined as χ2/(number of degrees of freedom), is 0.44, and the root-mean-squared deviation is 7.0 × 10-13 cm3 s-1. The fact that the reduced χ2 is much lower than 1.0 indicates that the model reproduces the data to within the experimental uncertainties. The best-fit values are 4.58 × 10-11 cm3 s-1, 6.78 × 10-27 cm6 s-1, and 0.226, respectively. The 68.3% confidence limits are 3.78 × 10-11 e k∞(cm3 s-1) e 6.03 × 10-11, 0.85 × 10-27 e k0(cm6 s-1) e 70.4 × 10-27, and 0.051 e Fcent e 0.68. The ratios of the maximum-to-minimum values of these confidence limits are 1.60, 82.8, and 13.3, respectively. To determine the confidence envelope which describes the uncertainty in our calculation of krec(P), we use each set of the independent parameters {k∞,k0,Fc} associated with the smaller 68.3% of the ∆χ2 values and calculate krec(P) over a pressure range of interest. This range may extend beyond the range of the data. At each pressure, Pj, used in the calculation we determine a kmin(Pj) and kmax(Pj). The curves defined by the locus of these points represent the confidence envelope for the calculation of the rate coefficient. An example of such a confidence envelope is shown in Figure 1. The dotted lines represent the 68.3% confidence envelope. In simple terms, 68.3% of the curves of krec(P) determined in the Monte Carlo simulations fall within this confidence envelope. The confidence envelope shown in Figure 1 tends to diverge as the pressure is increased or decreased. To better understand this behavior we compare the asymmetric Lorentzian broadening function, eq 5, with the data as a function of the reduced pressure, Pr. For this comparison we use eq 4 in eq 3 and take the decadic logarithm to obtain
(
Xexp(Pr) ) log
)
Pr kexp - log /log Fcent k∞ 1 + Pr
(6)
where Pr is given by k0Pexp/k∞RT; k0, k∞, and Fcent are the bestfit values determined by the least-squares analysis and kexp is the experimental rate coefficient at pressure Pexp. The comparison is shown in Figure 2. All of the experimental points have a reduced pressure greater than 1. Therefore, the data provide a significant amount of information about the highpressure rate coefficient, some information about the centering parameter, and very little information about the low-pressure rate coefficient. These trends are reflected in the ratios of the
Recombination of Methyl Radicals. 2
Figure 1. Results of least-squares analysis and Monte Carlo simulation of isothermal data at 577 K. The error bars represent 1σ experimental uncertainties, the solid line is the best fit, and the dashed lines represent the 68.3% confidence envelope.
J. Phys. Chem., Vol. 100, No. 3, 1996 987
Figure 3. Low-pressure rate coefficients determined by least-squares analysis of the isothermal data sets. Asymmetric Lorentzian broadening function, Fcent ≈ exp{-T(K)/300}. The error bars represent 68.3% confidence limits determined by Monte Carlo simulations.
IV. Temperature-Dependence of k0, k∞, and Results of Global Fits
Figure 2. Asymmetric Lorentzian broadening function from a leastsquares analysis of the 577 K isothermal data set. The error bars represent 1σ uncertainties of the experimental data.
maximum-to-minimum confidence limits given above. We may also note that since these data appear on only one side of the broadening function, it will be very difficult to distinguish between different broadening functions. In addition, although the functional form of the pressure-dependent behavior in the falloff region will have some influence on the value of the highpressure rate coefficient, the low-pressure rate coefficient will be determined completely by the choice of this function.
Before proceeding with the global fits, we must determine the temperature-dependent behavior of the limiting rate expressions k0(T) and k∞(T). To accomplish this, we choose an analytic expression for the pressure-dependent behavior in the falloff region, as we did in the previous section, and analyze several sets of isothermal data. Plots of the temperature-dependent rate coefficients will provide the functional form. These functions and their parameters will form the starting points for the global fit. Unfortunately, as we showed above, the form of the analytic expression will influence the behavior of the limiting rate coefficients. Therefore, we shall consider several different analytic expressions. First, we consider the suggestion of Gilbert et al., eqs 3-5. We use isothermal data sets at 296, 407, 513, 577, 700, 810, and 906 K. These contain data over a large enough pressure range that we may simultaneously determine the best-fit values of k0, k∞, and Fcent and their 68.3% confidence limits. When the best-fit values of the four highest temperatures are plotted, they do not show a uniform trend with temperature. Furthermore, their confidence limits are very large; several orders of magnitude. Since we wish to determine the temperaturedependent behavior of the low- and high-pressure rate coefficients, we used the three lowest temperature data sets to establish the temperature-dependent behavior of Fcent; Fcent(T) ≈ exp{-T(K)/300}. With this expression for the centering parameter, we then refit the above seven sets of data to determine k0(T) and k∞(T). The best-fit values are shown in Figures 3 and 4 along with their 68.3% confidence limits. Both rate coefficients may be described by Berthelot’s expression37
kj(T) ) Aj exp{-T/Tj}
(7)
where j ) 0 or ∞. This expression is a special case of van’t Hoff’s empirical four-parameter equation.38 This observation
988 J. Phys. Chem., Vol. 100, No. 3, 1996
Hessler and Ogren
Figure 4. High-pressure rate coefficients determined by least-squares analysis of the isothermal data sets. Asymmetric Lorentzian broadening function, Fcent ≈ exp{-T(K)/300}. The error bars represent 68.3% confidence limits determined by Monte Carlo simulations of the isothermal data. The solid line is the best-fit line of the calculated highpressure rate coefficient to A∞ exp{-T(K)/T∞} and the dashed lines represent the 68.3% confidence envelope from the global fit.
is not surprising since there is no barrier to methyl-methyl recombination and the temperature dependence of the limiting rate coefficients is expected to be weak.39 Now that we have established the temperature dependence of the limiting rate coefficients we proceed with the global fits. We assume Fcent(T) is given by its standard form:
Fcent(T) ) (1 - a)exp{-T(K)/T***}
(8)
and determine both a and T*** from the global fit. A global fit of the 217-point data set to determine six parameters gives
k∞(T) ) 7.71 × 10-11 exp{-T(K)/1298} cm3 s-1 (9a) k0(T) ) 8.06 × 10-26 exp{-T(K)/280} cm6 s-1 (9b) Fc(T) ) (1 - a)exp{-T(K)/333}
(9c)
The best-fit value of the parameter a is 0 to within the confidence limits of the fit. Therefore, we fixed a at 0 and recalculated the Monte Carlo simulations. The value of χ2 is 170.4, with five varied parameters the reduced χ2 is 0.804, and the rootmean-squared deviation is 8.90 × 10-13 cm3 s-1. The confidence envelopes and data for isothermal sets at 296, 577, 906, 1350, and 1525 K are shown in Figure 6. Scatter plots from the Monte Carlo simulations indicate there is significant correlation between the parameters and the distributions of the best-fit values from the Monte Carlo simulations are not symmetric. The 68.3% confidence limits are 7.14 × 10-11 e A∞(cm3 s-1) e 8.43 × 10-11 , 884 e T∞(K) e 2129, 4.9 × 10-26 e A0(cm6 s-1) e 11.5 × 10-26, 202 e T0(K) e 411, and 251 e T***(K) e 569. The maximum-to-minimum ratios of these limits are 1.2, 2.4, 2.3, 2.0, and 2.3, respectively.
Figure 5. Results of global fit and Monte Carlo simulation for the asymmetric Lorentzian broadening function, eqs 9. The error bars represent 1σ experimental uncertainties, the solid lines are the best fits at 296, 577, 906, 1350, and 1525 K, and the dashed lines represent the 68.3% confidence envelopes.
Wang and Frenklach3 have suggested that the above broadening function tends to predict a broadening factor which is too small at low and high pressures. They suggest the Gaussian blending function
{ [
X(Pr,T) ) exp -
]}
1 log{Pr} - R(T) 2 σ(T)
2
(10)
where a temperature-dependent shift is R(T) and a temperaturedependent width is σ(T). We have introduced the factor of 1/2 so the width agrees with the standard definition. Unfortunately, this broadening function contains two temperature-dependent parameters. To extract these from the data, we note the Lorentzian and Gaussian functions are approximately the same near the origin; their major differences occur in the long range tails of the functions. Therefore, we equate eq 10 to eq 5 at log{Pr} ) 0 and solve for R(T). This gives
R(T) ) R - 0.3T(K)/T***
(11)
where R0 ) 0.40 and T*** is defined in eq 8. When we equate eq 10 to eq 5 at log{Pr} - R(T) ) σ(T) and solve for σ(T), we obtain
σ(T) ) σ0 + 0.4T(K)/T***
(12)
where σ0 ) 0.54. Equations 11 and 12 may be viewed as the first two terms of power series expansions for the temperaturedependent shift and width functions. When we fix the value of T*** at 333 K, from eq 9c above, and determine k0(T) and k∞(T) from the isothermal data sets, we find the low-pressure rate coefficient is independent of temperature and the best-fit values of the high-pressure rate coefficient are within the confidence envelope shown in Figure 4. To perform a global fit, we characterize the broadening function by three independent parameters: T***, R0, and σ0. For the first fit we fix R0 at 0.40 and σ0 at 0.54 and vary four parameters: A∞, T∞, A0, and
Recombination of Methyl Radicals. 2
J. Phys. Chem., Vol. 100, No. 3, 1996 989 0.84. The distribution of the best-fit values from the Monte Carlo simulations is not Gaussian. The maximum-to-minimum ratios of the above confidence limits are 1.10, 1.43, 4.2, 1.10, and 1.6, respectively. When we performed a third fit by allowing R0 to vary the scatter plots showed the parameters are significantly correlated. In particular, the correlations between A0 and R0 and between σ0 and R0 are maximal. Therefore, when we tried to vary R0 to determine its best-fit value the Monte Carlo simulations showed that A0 could be as large as 3.7 × 10-24 cm6 s-1 and vary over 8 orders of magnitude. Clearly, it is very difficult to extract a low-pressure rate coefficient if both the shift parameter, R0, and the low-pressure rate coefficient are to be extracted from the current set of data. This observation is consistent with Figure 2. We simply do not have data in the pressure region needed to locate the shift parameter. Another empirical equation has been introduced by Gardiner.4 This is the a equation given by
ka ) (k0P/RT)a + k∞a
(15)
where a is a metric exponent which is temperature dependent. If a ) -1, eq 15 reduces to the Lindemann-Hinshelwood expression. The temperature dependence of the metric exponent may be approximated by Figure 6. Results for xJexp determined by least-squares analysis of the isothermal data sets. The low-pressure rate coefficient was held fixed at 3.4 × 10-26 cm6 s-1. The error bars represent the 68.3% confidence limit. The solid line is the result of the global fit given in eq 20c.
T***. For the 217-point data set the value of χ2 is 175.6 with a reduced χ2 of 0.824 and a root-mean-squared deviation of 9.03 × 10-13 cm3 s-1. The best-fit values of the parameters are
k∞(T) ) 7.78 × 10-11 exp{-T(K)/1007} cm3 s-1 -25
k0(T) ) 1.10 × 10
6 -1
(13a)
cm s
(13b)
Fc(T) ) exp{-T(K)/180}
(13c)
The 68.3% confidence limits are 7.48 × 10-11 e A∞(cm3 s-1) e 8.10 × 10-11 , 911 e T∞(K) e 1121, 0.91 × 10-25 e A0(cm6 s-1) e 1.34 × 10-25, and 174 e T***(K) e 187. The maximum-to-minimum ratios of these limits are 1.08, 1.23, 1.47, and 1.08, respectively. We use the results of this fit as the starting point for a second fit in which we extend the list of variable parameters to include σ0. This five-parameter fit reduces χ2 by 2% to 172.0, the reduced χ2 to 0.811, and the root-mean-squared deviation to 8.93 × 10-13 cm3 s-1. The bestfit parameters are
k∞(T) ) 7.68 × 10-11 exp{-T(K)/1095} cm3 s-1
(14a)
k0(T) ) 1.92 × 10-25 cm6 s-1
(14b)
Fc(T) ) exp{-T(K)/177}
(14c)
σ(T) ) 0.67 + 0.4T(K)/177
(14d)
The confidence envelopes for isothermal data at 296, 577, 906, 1350, and 1525 K are almost indistinguishable from the results shown in Figure 5. From the Monte Carlo simulations we obtain the 68.3% confidence limits: 7.30 × 10-11 e A∞(cm3 s-1) e 8.04 × 10-11, 935 e T∞(K) e 1340, 0.99 × 10-25 e A0(cm6 s-1) e 4.19 × 10-25, 169 e T***(K) e 186, and 0.52 e σ0 e
(
a(T) ) -1/ a* +
T T***
)
(16)
where a* and T*** are constants; T*** should be numerically equal to the same parameter in eq 8. From an analysis of several isothermal sets we conclude k0 is independent of temperature and k∞, again, follows eq 7. A least-squares analysis of the data set gives χ2 ) 206.9, a reduced χ2 of 0.971, and a rootmean-squared deviation of 9.8 × 10-13 cm3 s-1. The best-fit parameters are
k∞(T) ) 6.44 × 10-11 exp{-T(K)/3709} cm3 s-1
(17a)
k0 ) 8.27 × 10-27 cm6 s-1
(17b)
(
(17c)
a(T) ) -1/ 0.0 +
)
T(K) 189
The fact that a* is 0 to within experimental error is somewhat surprising since we had expected a* to be approximately 1.0. The 68.3% confidence limits of the other parameters are (6% for A∞, -40% and +200% for T∞, (25% for k0, and (55% for T***. Scatter plots of the Monte Carlo simulations indicate that A∞, k0, and T*** follow the parabolic approximation but that T∞ is asymmetric and correlated with the other three parameters. Confidence envelopes of the isothermal data are not as good as the other results, but they are reasonable. In 1989, Oref5 generalized Lindemann's expression for unimolecular rates and derived a “J equation”. We write this equation as
k(T,P) ) [-(k∞ + k0P/RT) + [(k∞ + k0P/RT)2 + 4(Jexp - 1)k∞k0P/RT]1/2/2(Jexp - 1)] (18) where Jexp is a temperature-dependent parameter to be determined by the data. When the three parameters k∞, k0, and Jexp are determined by the isothermal data sets, k0(T) is, again, independent of temperature near 10-26 cm6 s-1. We use this value for all temperatures to determine the temperaturedependent behavior of Jexp(T) and k∞(T). The high-pressure rate coefficient follows Berthelot’s expression and is only a few
990 J. Phys. Chem., Vol. 100, No. 3, 1996
Hessler and Ogren
Figure 7. Results of global fit and Monte Carlo simulation for the “J equation” of Oref, eq 18. The error bars represent 1σ experimental uncertainties, the solid lines are the best fits at 296, 577, 906, 1350, and 1525 K, and the dashed lines represent the 68.3% confidence envelopes.
percent different from the values determined with empirical broadening functions. The results for Jexp(T) are shown in Figure 6. Since the high-pressure rate coefficient follows Berthelot’s expression, we use eq 10 of Oref’s original paper5 to calculate the temperature-dependence of Jexp(T). The result is
Jexp(T) ) {exp[T(K)/TJ] - 1}2
(19)
where TJ is to be determined by an analysis of the data. The global fit of the 217-point data set to determine four parameters gives
k∞(T) ) 8.78 × 10-11 exp{-T(K)/723} cm3 s-1
(20a)
k0(T) ) 9.04 × 10-27 cm6 s-1
(20b)
Jexp(T) ) {exp[T(K)/268] - 1}2
(20c)
The value of χ2 is 171.8, the reduced χ2 is 0.807, and the rootmean-squared deviation is 8.93 × 10-13 cm3 s-1. The confidence envelopes for isothermal data at 296, 577, 906, 1350, and 1525 K are shown in Figure 7. Note that the confidence envelopes are narrower than in Figure 5 and do not diverge as the pressure is increased. Scatter plots from the Monte Carlo simulation indicate that the empirical parameters of this formulation are the least correlated and the most precisely determined. In fact, of the four this is the only situation where the χ2 hypersurface in the vicinity of the best-fit parameters is described by the parabolic approximation. The 68.3% (1σ) confidence limits are 3.6%. 5.8%, 20%, and 3.8% for the parameters A∞, T∞, A0, and TJ, respectively. V. Discussion To facilitate comparison between the various formulations we have summarized the parameters for the limiting rate
coefficients in Table 1. One might be tempted to compare the values of χ2 and conclude that the asymmetric Lorentzian function suggested by Gilbert et al. provides the best fit of the data and, therefore, should be used to extrapolate the behavior of the rate coefficient into regions where measurements have not been performed. Although this may be true, such a choice can not be based upon the values of χ2 alone. For example, all of the values of the reduced χ2 are less than unity, which indicates that all of the empirical formulas adequately reproduce the data. Therefore, one could argue that they are all valid and that more and better data are needed to provide a definitive test of these formulas. This is indeed true, and as we have shown in Figure 5, obtaining data at much lower pressures will provide important new information. Finally, when we consider the correlation between the parameters, the behavior of the χ2 function in the vicinity of the best-fit parameters, and the fractional values of the confidence limits it appears the “J equation” of Oref is better than the other formulations. We know of many instances where experimental data has been analyzed by the asymmetric Lorentzian formulation suggested by Gilbert et al. Clearly, the results presented here indicate that the Gaussian approach of Wang and Frenklach, the “J equation” of Oref, and Gardiner’s “a equation” should be tested with additional experimental results. Despite the fact that we cannot choose one formulation over the others, all of them indicate that the mathematical form of the temperature-dependent behavior of the high-pressure rate coefficient is best described by Berthelot’s formula, eq 7, which is a special case of van’t Hoff’s four-parameter equation:37
k(T) ) ATβ exp{-∆E/kBT}exp{DT}
(21)
The general use of the three-parameter temperature-dependent preexponential function suggested by Kooji (eq 21 with D ) 0), which has physical meaning for reactions which proceed over a barrier, has no physical foundation when the reaction has no barrier.39 When we used Kooji’s three-parameter equation and varied the parameters A, β, and ∆E/kB, the strong correlations between these parameters produce a very large confidence envelope for the prediction of the high-pressure rate coefficient. By using the four-parameter equation of van’t Hoff and setting ∆E/kB ) 0 (there is no barrier), the temperature dependence is described with the parameters β and D. In this case, we found that β was also equal to 0. Perhaps, a physical interpretation can be given for the parameter D. In Figure 8, we show the 68.3% confidence envelopes for the temperaturedependent behavior of the high-pressure rate coefficient. The relatively large area of the confidence envelope for the asymmetric Lorentzian broadening function is due to the significant correlation between the high-pressure rate coefficient and the broadening function. These correlations exist because the Lorentzian broadening function extends to much higher reduced pressure than the Gaussian function. An unfortunate aspect of Figure 8 is the clear indication that the temperature dependence of the high-pressure rate coefficient depends upon the formula used to describe the behavior in the falloff region. Surely, more work is needed in this area. There have been several calculations of the high-pressure rate coefficient for methyl-methyl recombination. The basic Gorin approach, developed in 1939, has recently been improved and extended by Pitt, Gilbert, and Ryan.9 The results of these latest calculations are shown as the open circles in Figure 8. An alternate approach is based upon the statistical adiabatic channel model.10,11 More recent calculations by Troe40 and Cobos and Troe41 use a simplified version of this model and suggest a slight positive temperature dependence, e.g., as the temperature is
Recombination of Methyl Radicals. 2
J. Phys. Chem., Vol. 100, No. 3, 1996 991
TABLE 1: Comparison of the Results of Various Formulationsa method asym Lorentz (5 variables) Gauss (4 variables) Gauss (5 variables) “a equation” (4 variables) “J equation” (4 variables) a
χ2 170 176 172 207 172
1011A∞ (cm3 s-1) +0.72 7.71-0.57 +0.32 7.78-0.30 +0.36 7.68-0.38
6.44 ( 0.39 8.78 ( 0.32
T∞ (K) +831 1298-414 +114 1007-96 +245 1095-160 +7418 3709-1485
723 ( 42
1026A0 (cm6 s-1) +3.4 8.1-3.2 +2.4 11.0-1.9 +22.7 19.2-9.3
T0 (K) +131 280-78
0.83 ( 0.21 0.90 ( 0.18
All confidence limits are given by Monte Carlo simulations and are at the 68.3% level.
confidence limits for the parameters of a model and confidence envelopes for predictions of the limiting rate coefficients. Monte Carlo simulations, or for that matter any other technique of data reduction, cannot answer a question if the answer is not contained within the data. For example, to distinguish between the various formulations of the behavior in the falloff region, data at lower pressures are needed. On the other hand, the Monte Carlo simulations provide realistic boundaries for the predictions. In certain situations these boundaries may be more important than the predictions.
Figure 8. Confidence envelopes (68.3%) for the high-pressure rate coefficient determined by the asymmetric Lorentzian (slanted lines), Gaussian with five variables, (vertical lines), “J equation” (cross-hatched lines), and “a equation” (uniform dots) formulations. The open circles are the latest calculations for a Gorin model,9 the solid line is the latest SACM calculation,41 and the solid dots are the latest RRKM calculations.13,42
increased from 300 to 1300 K the high-pressure rate coefficient should increase by 15%. These calculations are shown as the solid line in Figure 8. Most recently, Wagner and Wardlaw12 and Walter et al.13 have performed calculations based upon variational transition-state theory and recommended expressions for both the high- and low-pressure rate coefficients as well as the behavior in the falloff region. These calculations are particularly pertinent because they used the data of Slagle et al. to calibrate the calculations. Unfortunately, they did not correct the data at 700, 810, and 906 K; therefore, their calibration contains a systematic error. Of these, the latter13 is considered more reliable because it includes a substantial positive temperature dependence to 〈∆E〉tot, the average energy transferred in an Ar-C2H6* collision. The results of this calculation42 are shown as the solid circles in Figure 8. Note that although these results were initially presented in terms of Kooji’s equation they are described much better by Berthelot’s equation. VI. Conclusions The problems associated with extracting limiting rate coefficients from experimental data are well documented.29,43 To our knowledge, this is the first report of the use of Monte Carlo simulations to analyze kinetic data and, thereby, produce
Acknowledgment. We wish to thank Julia Willis, an undergraduate student from Knox College, for help with the data reduction. She participated in the Spring 1995 Science and Engineering Research Semester programs. This program is supported by the Division of Educational Programs. We thank Albert F. Wagner for providing us with the results of his calculations, helpful discussions, and critical comments on the manuscript. Other helpful discussions with Joe V. Michael, Michael J. Pilling, Struan H. Robertson, and Ju¨rgen Troe are acknowledged. We thank William C. Gardiner Jr. and Michael Frenklach for kindly providing us with references to their latest work. This work was supported by the U.S. Department of Energy, Office of Basic Energy Sciences, Division of Chemical Sciences. References and Notes (1) Du, H.; Hessler, J. P.; Ogren, P. J. J. Phys. Chem., previous paper in this issue. (2) Gilbert, R. G.; Luther, K.; Troe, J. Ber. Bunsen-Ges. Phys. Chem. 1983, 87, 169. (3) Wang, H.; Frenklach, M. Chem. Phys. Lett. 1993, 205, 271. (4) Gardiner Jr., W. C. In 12th IMACS World Congress on Scientific Computations; Paris, France, 1988; p 582. (5) Oref, I. J. Phys. Chem. 1989, 93, 3465. (6) Robertson, S. H.; Wardlaw, D. M.; Hirst, D. M. J. Chem. Phys. 1993, 99, 7748. (7) Wardlaw, D. M., private communication. (8) Gorin, E. J. Chem. Phys. 1939, 7, 256. (9) Pitt, I. G.; Gilbert, R. G.; Ryan, K. R. J. Phys. Chem. 1995, 99, 239. (10) Quack, M.; Troe, J. Ber. Bunsen-Ges. Phys. Chem. 1977, 81, 329. (11) Quack, M.; Troe, J. Ber. Bunsen-Ges. Phys. Chem. 1974, 78, 240. (12) Wagner, A. F.; Wardlaw, D. M. J. Phys. Chem. 1988, 92, 2462. (13) Walter, D.; Grotheer, H.-H.; Davies, J. W.; Pilling, M. J.; Wagner, A. F. in Twenty-Third Symp. (Int.) on Combustion; The Combustion Institute: Seattle, 1990; p 107. (14) Chien, J.-Y. J. Am. Chem. Soc. 1948, 70, 2256. (15) Talat-Erben, M. J. Chem. Phys. 1957, 26, 75. (16) Ogren, P. J.; Hessler, J. P. Int. J. Chem. Kinet. 1995, 27, 719. (17) Berman, P. R. in AdVances in Atomic and Molecular Physics; Bates, D. R., Benderson, B., Eds.; Academic Press: New York, 1977; Vol. 13, p 57. (18) Rea, E. C.; Chang, A. Y.; Hanson, R. K. In Western States Combustion Meeting; unpublished, 1989. (19) Chang, A. Y.; Di Rosa, M. D.; Hanson, R. K. J. Quantum Spectrosc. Radiat. Transfer 1992, 47, 375. (20) Kessler, W. J.; Allen, M. G.; Davis, S. J. J. Quantum. Spectrosc. Radiat. Transfer 1993, 49, 107. (21) Macpherson, M. T.; Pilling, M. J.; Smith, M. J. C. J. Phys. Chem. 1985, 89, 2268. (22) Hippler, H.; Luther, K.; Ravishankara, A. R.; Troe, J. Z. Phys. Chem. NF (Frankfurt) 1984, 142, 1.
992 J. Phys. Chem., Vol. 100, No. 3, 1996 (23) Hwang, S. M.; Wagner, H. G.; Wolff, T. In Twenty-Third Symp. (Int.) Combustion; The Combustion Institute: Seattle, 1990; p 99. (24) Slagle, I. R.; Gutman, D.; Davies, J. W.; Pilling, M. J. J. Phys. Chem. 1988, 92, 2455. (25) Gla¨nzer, K.; Quack, M.; Troe, J. Chem. Phys. Lett. 1976, 39, 304. (26) Gla¨nzer, K.; Quack, M.; Troe, J. In Sixteenth Symp. (Int.) Combustion; The Combustion Institute: Seattle, 1977; p 949. (27) Troe, J., private communication, 1994. (28) Zaslonko, I. S.; Smirnov, V. N. Kinet. Catal. 1979, 20, 575. (29) Pilling, M. J. Int. J. Chem. Kinet. 1989, 21, 267. (30) Bevington, P. R.; Robinson, D. K. Data Reduction and Error Analysis for the Physical Sciences, 2nd ed.; McGraw-Hill Book Co., Inc.: New York, 1992. (31) Press, W. H.; Flannery, B. P.; Teukolsky, S. A.; Vetterling, W. T. In Numerical Recipes: The Art of Scientific Computing; Cambridge University Press: New York, 1989; p 498. (32) Box, G. E. P.; Muller, M. E. Ann. Math. Stat. 1958, 19, 610. (33) Lampton, M.; Margon, B.; Bowyer, S. Astrophys. J. 1976, 208, 177.
Hessler and Ogren (34) Avni, Y. Astrophys. J. 1976, 210, 642. (35) Meyer, S. L. In Data Analysis for Scientists and Engineers John Wiley & Sons, Inc.: New York, 1975; p 359. (36) Pawlowska, Z.; Gardiner Jr., W. C.; Oref, I. J. Phys. Chem. 1993, 97, 5024. (37) Laidler, K. J. In Chemical Kinetics; Harper & Row: New York, 1987; p 39-48. (38) van’t Hoff, J. H. Lectures on Theoretical and Physical Chemistry; Edward Arnold: London, 1898; Part 1. (39) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell Scientific Publications: Oxford, 1990. (40) Troe, J. J. Chem. Phys. 1981, 75, 226. (41) Cobos, C. J.; Troe, J. J. Chem. Phys. 1985, 83, 1010. (42) Wagner, A. F., private communication. (43) Troe, J. In Twenty-Second Symp. (Int.) Combustion; The Combustion Institute: Seattle, 1988; p 843.
JP951218O
