J. Phys. Chem. 1986, 90, 2616-2634
2616
LASER CHEMISTRY, MOLECULAR DYNAMICS, AND ENERGY TRANSFER The Effect of VlbratlonaI-Rotational Dlsequlllbrlum on the Rate Constant for an Atom-Transfer Reaction Carmay Lim and Donald G. Truhlar* Department of Chemistry, University of Minnesota, Minneapolis, Minnesota 55455 (Received: August 8, 1985; In Final Form: January 13, 1986)
We report the first realistic quantitative calculation of the competition between reaction and vibrational-rotational energy transfer in respectively perturbing and reestablishing the internal-state equilibrium of reactants and products in a typical, fast, second-order bimolecular reaction. The reaction considered is C1 HBr HCl + Br at 300 K, including V-V,R,T; V-V,T; V-R,T; V-T; R-R,T; and R-T energy-transfer processes. The 110-state nonlinear master equation is solved by a variety of techniques, including numerical integration, iteration to a quasi-steady state, eigenvalue-eigenvector analysis of the symmetric matrix obtained by linearization in the deviations of the state concentrations from their values at chemical equilibrium, and eigenvalue-eigenvector analysis of the nonsymmetric matrix obtained by linearization in the deviations of the state concentrations from their values for typical experimental initial conditions. We show that vibrational-rotational nonequilibriumeffects are less than 1%under typical experimental conditions (i.e., initial rate studies in the absence of significant back reaction), but that the conventional picture of reactive time scale > vibrational relaxation time scale > rotational time scale is not valid. The eigenvalue of the rate matrix corresponding to the eigenvector that carries most of the reactive flux is the fourth or fifth smallest nonzero eigenvalue, in contrast to the conventional picture in which it is the smallest nonzero eigenvalue. As a consequence, after the early-time quasi-steady state in which the local equilibrium rate coefficient is observable (within l%), the reaction rate coefficient decreases and phenomenological analysis of the net rate of change of the concentration of reactants may show up to four additional quasi-steady states involving successively slower reaction rates before chemical equilibrium is finally achieved.
+
I. Introduction Since chemical reaction rates are not measured under conditions of complete chemical equilibrium, and since chemical reactions perturb the internal-state population distributions of both reactants and products, even a small perturbation from chemical equilibrium puts one in the regime of nonequilibrium kinetics. Gas-phase collision theory and transition-state theory, however, are most readily used to calculate the local-equilibrium rate constant which is an average of the initial-stateselected, fmed initial translational energy reaction rate probabilities over an equilibrium distribution of reactant In contrast the measurable phenomenological rate constant depends on the nonreactive energy-transfer probabilities in the system as well as on the reaction pr~babilities.~” It is of great interest for applications of collision theory to gasphase rate processes to estimate the deviation of phenomenological rate constants from local equilibrium ones under various conditions and for various classes of reaction. It is especially important to consider fast reactions, for which the nonequilibrium effects are largest. Nonequilibrium effects on dissociation reactions, both unimolecular and bimolecular, have been widely studied.622 Most (1) M. A. Eliason and J. 0. Hirschfelder, J. Chem. Phys., 30,1426 (1959). (2) J. Ross, J. C. Light, and K. E. Shuler, in Kinetic Processes in Gases and Plasmas, A. R. Hochstim, Ed., Academic Press, New York, 1969, p 281. (3) M. M. Kreevoy and D. G. Truhlar, in Investigation of Rates and Mechanisms of Reactiom, Part 1,4th ed.,C. F. Bernasconi, Ed., Wiley, New York, 1986, p 13. (4) B. Widom, Science, 148, 1555 (1965). ( 5 ) R. K. Boyd, Chem. Rev., 77, 93 (1977). (6) H. 0. Pritchard, The Quantum Theory of Unimolecular Reactions, Cambridge University Press, Cambridge, 1984. (7) P. J. Robinson and K. A. Holbrook, Unimolecular Reactions, Academic Press, New York, 1972.
0022-3654 /86/2090-26 16%01SO10 I
,
-+
work on nonequilibrium effects on nondissociation reactions has centered on translational d i s e q u i l i b r i ~ m , ~usually ~ - ~ ~ using approximate solutions of the Boltzmann equation. Internal-state disequilibrium effects may be addressed in terms of the discrete-index master equation. Eyring and co-workers presented (8) W. Forst, Theory of Unimolecular Reactions, Academic Press, New York, 1973. (9) J. Troe, in Physical Chemistry, An Advanced Treatise, Vol. VIB, W. Jost, Ed., Academic Press, New York, 1975, p 835. (10) D. C. Tardy and B. S. Rabinovich, Chem. Rev., 77, 369 (1977). (1 1) B. J. Gaynor, R. G. Gilbert, and K. D. King, Chem. Phys. Lett., 55, 40 (1978). (12) D. C. Tardy and B. J. Malius, J. Phys. Chem., 83, 93 (1979). (13) W. Forst, J. Phys. Chem. 83, 100 (1979). (14) J. Troe, J. Phys. Chem. 83, 114 (1979). (15) A. J. Stace and P. V. Sellers, Chem. Phys., 50, 147 (1980). (16) H. 0. Pritchard and S. R. Vatsya, J . Chem. Phys., 59,2575 (1981). (17) A. Lifschitz, A. Bar-Nun, A. Burcat, A. Ofir, and R. D. Levine, J . Phys. Chem., 86, 791 (1982). (18) N. Snider, J . Chem. Phys., 77, 789 (1982). (19) H. W. Schranz and S. Nordholm, Chem. Phys., 74, 365 (1983). (20) H. 0. Pritchard, Spec. Period. Rep.: React. Kine?., 1, 243 (1975). (21) D. G. Truhlar, N. C. Blais, J.-C. J. Hajduk, and J. H. Kiefer, Chem. Phys. Lett., 63, 337 (1979). (22) J. E. Dove and S. Raynor, J . Phys. Chem., 83, 127 (1979). (23) C. F. Curtiss, Ph.D. Thesis, Wisconsin University Report CM-476 (1948). (24) I. Prigogine and E. Xhrouet, Physica (Amsterdam), 15, 913 (1949). (25) I. Prigogine and M. Mathiew, Physica (Amsterdam), 16, 51 (1950). (26) R. D. Present, J . Chem. Phys., 31, 747 (1959). (27) R. D. Present and B. M. Morris, J . Chem. Phys., 50, 151 (1969). (28) B. Shizgal and M. Karplus, J. Chem. Phys., 52,4262 (1970); 54,4345 (1971); 54, 4357 (1971). (29) B. Shizgal, J . Chem. Phys., 55, 76 (1971). (30) R. Kaprai, S . Hudson, and J. Ross, J. Chem. Phys., 53,4387 (1970). (31) K. Koura, J . Chem. Phys., 59, 691 (1973). (32) N. Xystris and J. S . Dahler, J . Chem. Phys., 68, 354 (1978).
b 1986 American Chemical Societv
Effect of Vibrational-Rotational Disequilibrium early examples of such treatments for irreversible few-state first-order or pseudo-first-order models of r e a ~ t i o n . ~Poulsen ~,~~ and Shizgal have presented more realistic examples for slow four-center exchange reaction^.^^^^^ There has been less work on atom-transfer reactions, but Finkelman and Dove studied vibrational-state models in which rotational degrees of freedom are not treated e ~ p l i c i t l y , ~and ' * ~ Gutkowicz-Krusin ~ has studied a model in which the reactants have two states and the products have 01-18.~~In addition to the papers presenting calculations, several authors have provided general discussions of the effect on reaction rates of nonequilibrium internal-state distributions of reactants or p r o d ~ c t s . ~ ~ ~ ~ In the present study we attempt to make as realistic as possible an estimate of the effect of internal-state disequilibrium on the rate constant for a fast, bimolecular, reversible atom-transfer reaction. The reaction chosen is
C1 + HBr
1
7HC1 + Br
AHo,, = -15.7 kcal mol-I
(R1) at 300 K. Various experimental studies have yielded a room temperature rate constant for this reaction of 7 X 10-l2, 8 X 10-l2, or 1 X lo-" cm3 molecule-' s-', re~pectively.~'-~~ (The energy of activation is only 0.848or 0.949kcal mol-'.) Reaction R1 is an interesting case because the forward reaction rate is reasonably large, though not atypical of exothermic reactions, and because sufficient data are available to make reasonable estimates of all the required state-to-state rate coefficients. Throughout the study we assume translational energy is well equilibrated. Some aspects of the present study have been reported elsewhere in two preliminary communication^,^^-^^ and additional details are available in a thesis.52 Section I1 presents the theoretical framework for this study and the notation for reactive and nonreactive rate coefficients. Section I11 summarizes the state-to-state rate coefficients we used for reactive and nonreactive proce'sses. Section IV describes the methods used to calculate the phenomenological rate coefficients. Section V discusses the results, and section VI is a brief summary. 11. Theory We shall calculate the phenomenological rate constants k l and k-] which are defined such that they obey the usual rate law d[HCl]/dt = k,[Cl][HBr] - k-,[HCl][Br] (2.1)
and the rate quotient law where K l e is the equilibrium constant for ( R l ) in the exothermic direction. W e assume isothermal, isobaric conditions at a temperature of 300 K. We then consider three general questions: (i) After an initial transient period, does the phenomenological rate (33) R. J. Zwolinski and H. Evrine. J . Am. Chem. SOC..69.2702 (1947). (34j S. J. Yao and B. J. Zwolin&, Adu. Chem. Phys., 21, 91 (1971j. (35) L. L. Poulsen, J . Chem. Phys., 53, 1987 (1970). (36) B. Shizgal, J . Chem. Phys., 57, 3915 (1972). (37) M. Finkelman, Ph.D. Thesis, University of Toronto, Canada, 1976. (38) T. C. Clark, J. E. Dove, and M. Finkelman, Acta Astronaut., 6,961 (1979). (39) D. Gutkowicz-Krusin, Physicu A, 97A,425 (1979). (40) H.Eyring, T. S. Ree, T. Ree, and F. M. Wanslass, Spec. Pub1.Chem. SOC.,16,25 (1962). (41) B. Widom, J . Chem. Phys., 55,44 (1971). (42) B. Widom, J . Chem. Phys., 61,672 (1974). (43) N.S . Snider and J. Ross, J . Chem. Phys., 44, 1087 (1966). (44) K. H.Lau, S. H. Lin, and H. Eyring, J. Chem. Phys., 58, 1261 (1973). (45) S.A. Reshetnyak and L. A. Shelepin, Tr. Fir. Imt. i m P. N. Lebedeua, Adad. Nauk SSSR, 106.90 (1979). (46) A. I. Osipov, Inzh-Fiz. Zh., 38, 351 (1980). (47) K. Bergmann and C. B. Moore, J . Chem. Phys., 63, 643 (1975). (48) C.-C. Mei and C. B. Moore, J . Chem. Phys., 67,3936 (1977). (49) R.Rubin and A. Persky, J . Chem. Phys., 79,4310 (1983). (50)C. Lim and D. G. Truhlar, J . Phys. Chem., 89, 5 (1985). (51) C. Lim and D. G. Truhlar, Chem. Phys. Lett., 114, 253 (1985). (52) C. Lim, Ph.D. Thesis, University of Minnesota, Minneapolis, 1984, Chapter 4.
The Journal of Physical Chemistry, Vol. 90, No. 12, 1986 2617
law hold with constant rate coefficients (a regime for which the phenomenological rate law holds with constant or nearly constant rate coefficients will be called a quasi-steady state or a steadykinetics regime) or do the rate coefficients change continuously throughout an experiment and never settle down? (ii) If the rate coefficients do not depend on the extent of reaction during a single experiment, to what extent are they independent of the initial composition of the reactive system? (iii) To the extent that the phenomenological rate coefficients are independent of extent of reaction and initial composition, how do they compare to the equilibrium rate constants? The first two questions examine the validity of the phenomenological rate law for bimolecular reactions, and the last question gives a quantitative estimate of the deviations of the phenomenological rate constants from their local equilibrium estimates. In order to obtain quantitative answers to these questions we obtained exact and approximate solutions to the nonlinear master equation for the system. The master equation includes reactive collisions; V-V,T and V-V,R,T collisions (in which vibrational quantum numbers change on two molecules in a single collision, the translational energy changes, and possibly also one or more molecules change their rotational energy); V-T and V-R,T collisions (in which the vibrational quantum number changes on one molecule and the vibrational energy change is balanced by a change in translational energy or changes in both the rotational energy of the same molecule and translational energy); R-R,T collisions (in which rotational quantum numbers changes on two molecules in a single collision and the translational energy changes); and R-T collisions (in which the rotational energy of one molecule and the translational energy change). No electronically excited species are considered. For all the calculations we assume that Ar is present and that it acts as a third body. To simplify the notation in the following equations, we denote C1, HBr, HCl, Br and Ar by the letters A, B, C, D, and M, and state-to-state rate constants are written with initial quantum numbers in subscripts and final ones in superscripts. We denote the concentrations of HBr and HCl in the states with vibrational and rotational quantum numbers u a n d j by nBpJ and nc,uj, respectively. Then the master equation for ( R l ) is given by the coupled set
2618
The Journal of Physical Chemistry, Vol. 90, No. 12, 1986 [Cl,Br] = [Cl]
+ [Br]
(2.7)
Lim and Truhlar where the xe's represent normalized Boltzmann distributions and kIe/k-le = Kle
[HBr,HCl] denotes the total concentration of HBr and HC1 [HBr,HClI =
CC(n~,uj + U J
(2.8)
nc,UJ)
and [N] denotes the total concentration of all species [N] = [Ar]
+ [HBr,HCl] + [Cl,Br]
where K I Cis the equilibrium constant for reaction (Rl). The phenomenological rate constants are defined by requiring that they satisfy (2.1) and (2.2). Thus kl can be calculated from
(2.9a)
We will also use the notation [R,P] to denote the total concentration of reactants and products [R,P] = [Cl,Br] + [HBr,HCI] = [N] - [Ar] (2.9b) To obtain eq 2.3-2.6 we neglected termolecular recombination of C1 and Br to form molecules, wall reactions, and secondary reactions of the products, such as the C1+ HCl reaction (although back reaction is included). Thus, reaction R1 is treated as a single-step, reversible second-order reaction in the absence of competitive or side reactions. Notice that because of these assumptions, [Cl,Br], [HBr,HCl], and [N] are all constants (independent of time) in eq 2.7-2.9. We also made the following additional simplifying approximations in the treatment of energy-transfer processes: For V-V,R,T and V-V,T transitions in which the collision partner is HBr, we neglected V-R,T transitions in which both vibrational and rotational quantum numbers of HBr change in a single collision, except for V-V,R,T collisions with HCl in which the collision partner also changes its vibrational quantum number. For V-T transitions of HBr, we assumed identical rate constants for collisions with either Cl or Br. For R-T transitions of both molecules we assumed identical rate constants for all collision partners. For HC1, we neglected V-R,T transitions except for collisions with C1. For V-T transitions of HCl we assumed that the state-testate rate constants for collisions with Ar, HBr, or another HCl differ only by overall scale factors aAr,aHBr, and aHCI,independent of u j , and v', from those for = aH ]. We also assumed collisions with Br, and further that aHBr that the forward state-testate reaction rates k$$'are independent of the state u, j of HBr. These assumptions are made mainly because of lack of more detailed data; they simplify the master equation relatively little. Furthermore, Keren et aLs3have argued that the fine details of the rotational distributions arising from V-R,T transfer are largely erased by the generally much faster R-R,T and V-V,T processes; thus the lack of detailed information about the initial and final rotational quantum number dependence of the vibrational energy-transfer rate constants is not a serious limitation. Three sets of reaction rate coefficients are of interest:' the one-way flux coefficients kIf and k-lf, the equilibrium rate constants kle and kmIe,and the phenomenological rate constants kl and k-l. We can rewrite the left-hand side of (2.1) or (2.3) in terms of the reactive fluxes to give d[Br]/dt = R I f- R-lf = klf[HBr][C1] - k..lf[HCl][Br] (2.10)
k, =
k$$~J'xB,vj
(2.1 1)
B . u j C,u'J'
(2.12) where the x's represent the internal-state distributions, normalized separately to unity for each molecule, a t time t . In the limit as t tends to infinity the one-way flux coefficients become equal to the equilibrium rate constants, which are given by
d[HCl]/dt [Cl] [HBrIG
(2.16)
and k-l can be calculated from (2.2) and (2.16). In eq 2.16, G is the unitless reaction affinity for ( R l ) and is defined by G= 1-
[HCl][Br]/([HBr][C1]Kle)(2.17) = 1 - R-I/RI
(2.18)
where R1 and R-, are the forward and reverse phenomenological rates, and are given by
R I = k,[HBr][Cl]
(2.19)
R-l = k-, [HCl] [Br]
(2.20)
The deviations of the one-way flux coefficients from the equilibrium rate constants are represented by 71
= klf/kle - 1
7-1 =
k-lf/k-le - 1
(2.21) (2.22)
and the deviations of the forward and reverse phenomenological rate constants from their local equilibrium estimates are measured by RIN = K - 1 (2.23) where we have defined K
K
to be
= k,/kle = k-l/k-le
(2.24)
The second equality in (2.24) results from (2.2). We shall also calculate the relative fractional state populations of HBr and HCl fB,vj
=
XB,vj/XB,vJe
fc,uj = XC,uj/XC,aje
(2.25) (2.26)
111. System
We consider 74 vib-rotational states of HC1, 34 vib-rotational states of HBr, and the ground electronic states of C1(2P3iz)and Br(2P3/2)for a total of 110 concentration variables in the master equation. All states up to an energy of 31 kcal mol-l from the bottom of potential of well of HCl are included in our calculations; for HCl the states included are j = 0-25 in the u = 0 manifold, j = 0-20 in u = 1, j = 0-16 in u = 2 , and j = 0-9 in u = 3, and for HBr those included are j = 0-20 in u = 0 and j = 0-12 in u = 1. The energy levels of HC1 and HBr are calculated from
Equating the right-hand sides of (2.3) and (2.10) we get klf =
(2.15)
EUJ= EZb + EL,J
(3.1)
where
+ X) + X)' Eujroi = BJjCi + 1) + DJ'jCi + =
W,(U
W,X,(U
(3.2) (3.3)
(2.14)
In these equations we is the vibrational frequency, ugeis the first anharmonicity constant, B, is the rotational constant, and D, is the centrifugal distortion constant.54 The equilibrium constant is calculated from the energy levels of eq 3.1-3.3 and the reaction enthalpy (AHD, = -15.7 kcal mol-'). We assume no population of electronically excited states. At 300 K the equilibrium constant K l e is 2.178 X IO".
(53) E.Keren, R. B. Gerber, and A. Ben-Shad, Chem. Phys., 21, 1 (1976).
(54) G. Herzberg, Spectra of Diatomic Molecules, Van Nostrand Reinhold, Princeton, NJ, 2nd ed., 1950.
(2.13)
The Journal of Physical Chemistry, Vol. 90, No. 12, 1986 2619
Effect of Vibrational-Rotational Disequilibrium Most of the rate coefficients needed for the present study were collected by Keren et al?3 in conjunction with a study of rotational disequilibrium under steady photolysis conditions in the C1+ HBr pulsed laser. In some respects, however, the data they used were insufficient for the present study and had to be extended. Here we summarize the data we actually used. State-to-state rate constants for V-V,R,T and V-V,T energy transfer in HC1 collisions with HC1 or HBr, V-R,T energy transfer in HC1-Cl collisions, R-R,T and R-T energy transfer in HCl 0 and 2 1 V-T deacticollisions with any partner, and 1 vations in HC1-Br collisions were taken from Keren et R-R,T and R-T rate constants for HBr with any collision partner were assumed to have the same form as those for HCl with the preexponential factor scaled down by a factor of 4.5, to reproduce the ratio of collision numbers of Agrawal and S a ~ s e n a . The ~~ 3 2 rate constant for V-T energy transfer in HCl-Br collisions 1 value, based on the was assumed to be 1.5 times the 2 Landau-Teller relation.56 The scale factors for V-T rate constants for HC1-Ar and HC1-molecule collisions are taken to be aAr= and aHBr = 0.1, based on the compilation of Leone,57 and the V-T rate constants for HBr are taken from that compilation and from the assumptions that the efficiency of HCl as collision partner is independent of state and that of C1 is the same as that of Br. For the reactive process HBr(uj) C1 s HCl(u'j') Br (R1)
-
-
+
-
+
+
the stateto-state rate constants for the forward direction are given by kc,V'J'= klCpc,dpC,V'J (3.4) B.UJ
where k l c ,the equilibrium rate constant in the exothermic direction, is given the value used by Keren et al.,53 namely 7.6 X cm3 molecule-' s-l. Combining this with our calculated equilibrium constant yields a backward rate constant kleof 3.489 X cm3 molecule-' s-l. pC+"and pc-"'J'are the fractions of reactive collisions which populate HC1 in vibrational level u'and in rotational statej'of that level, respectively, and are taken from Keren et ales3and extended to the full set of states by surprisal Details of the surprisal analysis and a complete table of the state-to-state reaction rate constants are given elsewhere.52 The fraction of reactive collisions that populate HC1 in vibrational levels u = 0, 1, 2, and 3 are 0.18, 0.59, and 0.24, and 6.9 X respectively. The ratio ~ , 2 / ~ = , 12.5 is in reasonable agreement with a more recent determination,61which yields 3. Reverse inelastic and reaction rate constants are obtained by microscopic reversibility. The relations in section I1 are independent of the specific form of kE;$'. If we specialize to the present case in which eq 3.4 assumes that kg;:iJ'is independent of u a n d j of HBr, eq 2.1 1 and 2.13 yield the additional relation k l f= k l e
(3.5)
IV. Methods The master equation for (Rl), eq 2.3-2.6, is a coupled set of ordinary differential equations which are first-order (with respect to the time derivative), nonlinear (in the concentration terms), and stiff, i.e., have a set of inherent time constants that span many orders of magnitude. The nonlinearity precludes the possibility of an analysis along the lines of the analytic studies performed for isomerization62d4or dissociation reactions.6H8 Furthermore, P. M. Agrawal and M. P. Saksena, J . Chem. Phys., 61,848 (1974). L. Landau and E. Teller, Phys. Z . Sowjerunion, 10, 34 (1936). S. R. Leone, J . Phys. Chem. Ref.Data, 11, 953 (1974). R. D. Levine and J. Manz, J . Chem. Phys., 63, 4280 (1975). H. Kaplan, R. D. Levine, and J. M a n , Chem. Phys., 12,447 (1976). (60) E. Pollak and R. D. Levine, Chem. Phys. Lett., 39, 199 (1976). (61) M. A. Wickramaaratchi and D. W. Setser, J . Phys. Chem., 87, 64 (1983). (62) J. Wei and C. D. Prater, Adu. Catal., 13, 203 (1962). (63) N. S.Snider, J . Chem. Phys., 42, 548 (1965). (64) B. Widom, Science, 148, 1555 (1965).
the stiffness of the master equation presents difficulties since these stiff equations are characterized by the presence of transient components that, although negligible relative to the other components of the numerical solution, constrain the step size of traditional numerical methods (e.g., Adams of Runge-Kutta formulas) to be of the order of the shortest time scale in the problem. However, several methods are now available for the solution of stiff equations, and we shall employ a variableorder method based on backward differentiation multistep formulas, as originally analyzed and implemented by Gear69 and later modified and studied by Hindmarsh and B ~ r n e . ' ~We ~ ~use the EPISODE program of Hindmarsh and Byme70-72to solve the master equation for the present study, but we also examine three other methods, which, although less general, offer significant advantages in terms of computational efficiency (the second method, which is based on simple iteration) or interpretability (the third and fourth methods, which involve eigenvalue-eigenvector analyses). In previous we have shown that a nonstiff iterative method provides an efficient and convenient way to obtain phenomenological rate constants for dissociation reactions, even when the associated master equation is stiff. Here we shall use the iterative technique as the second method to calculate quasisteady-state distribution and phenomenological rate constants for (Rl). Although we cannot solve the master equation analytically for (Rl), we can find analytic solutions to the linearized master equation obtained by expanding the original master equation to first order in the deviations of the concentrations from their values at chemical equilibrium. This is the third method used here, and we shall compare the results of eigenvalue-eigenvector analyses of linearized master equations to those obtained by the stiffequation solver and the iterative method. The fourth method used here is based on eigenvalue-eigenvector analysis of the master equation linearized about initial conditions. In the rest of this section we shall outline the procedures employed to implement these four methods. A . Numerical Integration. We used the EPISODE program to integrate the master equation for the time dependences of the 110 concentration variables for several sets of initial conditions. In this subsection we discuss the input data which we use in EPISODE. We mentioned above that [Cl,Br], [HBr,HCl], and [N] are all conserved quantities for the reaction considered here. Two additional conserved quantities are defined by YCl = ([Cll + [HClI)/[Nl (4.1) YBr = (IBr] -k [HBrl)/[Nl (4.2) We study systems with [HCl] = [Br], which would be the case for systems composed initially of only C1 and HBr, but we start the integration with products present. For a given ratio of the forward to reverse phenomenological rates R 1 / R l and a given set of Ycl and YBr,we calculate the total concentrations of Br, C1, HBr, and HC1 by solving a quadratic equation. We then assume that the initial internal-state distributions of HCl and HBr are Boltzmann distributions at 300 K, and this allows us to calculate the initial values of the concentrations nc,uJ and nB,uJ, which, in conjunction with [Cl] and [Br], form a concentration vector of length 110. Given the appropriate input parameters and the initial concentration vector, EPISODE computes the new concentration vector at each desired output value of t . Hence, the various rate coef(65) C. Lim and D. G. Truhlar, J . Phys. Chem., 87, 2683 (1983). (66) C. A. Brau, J . Chem. Phys., 47, 1153 (1967); 47, 3076 (1967). (67) W. L. Hogarth and D. L. S. McElwain, J . Chem. Phys., 63, 2502 (1975). (68) W. L. Hogarth and D. L. S. McElwain, Chem. Phys., 19,429 (1977). (69) C. W. Gear, Numerical Initial Value Problems in Ordinary Differential Equations, Prentice-Hall, Englewood Cliffs, NJ, 1971. (70) G. D. Byrne and A. C. Hindmarsh, ACM Trans. Math. Software, 1, 71 (1975). (71) A. C. Hindmarsh and G. D. Byrne, Lawrence Livermore Laboratory Report UCID-30112 (1975). (72) A. C. Hindmarsh and G. D. Byrne, in Numerical Methods for Differential Systems, L. Lapidus and W. E. Schiesser, Eds., Academic Press, New York, 1976, pp 147-166. (73) C. Lim and D. G. Truhlar, J . Phys. Chem., 88, 778 (1984).
2620 The Journal of Physical Chemistry, Vol. 90, No. 12, 1986
ficients can be calculated a t each output t value from the corresponding concentration vector at that time and its first-order time derivative. B. Iterative Method. The iterative method has been previously described in ref 65. In the present case we use it to converge to the first one or two quasi-steady states (where k , and k-, have become constant) for a system commsed initially of onlv CI and HBr according
Lim and Truhlar finite. Furthermore, the conservation relations imply that there are three linear dependences in the 110 rate equations; thus S has three zero eigenvalues. The eigenvalues A, and eigenvectors 12,) of S satisfy the equations SIR,) = A,lZ,)
(4.15)
(Z,lZ”) = 6,”
(4.16)
The solution of (4.14) is easily shown to be
Iw) = Cbre-hFr12,)
(4.17)
L
where the time-independent coefficients b, are determined from b, = (Z,lw)o
(4.18)
and 1 ~ =) Iw) ~ at time t = 0. Writing (4.17) in terms of components wi and substituting (4.12)-(4.13), we have
ni = nie + Cb,e-AJXi,, where [Cl], denotes initial concentration of C1, the superscript ( i )denotes the ith iteration, and Ar is the time step. The derivatives on the right-hand side of (4.3) and (4.4) are computed according to eq 2.4 and 2.5, respectively. The first equality in (4.5) is a consequence of our choice of initial conditions, and (4.6) results from the conservation of atoms in the system. We first assume a set of n&J, n&. [Br]lol, and [C1](olcorresponding to a given value of (Rl/Rl)$l.These initial concentrations along with the state-to-state rate constants and an approximate time step are substituted into the right-hand side of (4.3) and (4.4) to find the second approximation of the internal-state concentrations of HBr and HC1 which, in turn, give the second approximation to the total concentrations using (4.5)-(4.7). This process is repeated until the phenomenological rate coefficients becomes constant within an acceptable tolerance. The choice of time step is arbitrary, and for the present study s. The iterative method is able we use a time step, At = 5 X to converge to a quasi-steady state very efficiently, but then it takes much longer to follow the time evolution of the concentrations toward their equilibrium values as stiffness arises. Consequently, we do not follow a set of initial concentrations for all the way to chemical equilibrium; instead a given (R,/R-l)lOl we switch to EPISODE when it is desired to follow the concentrations all the way to chemical equilibrium. C. Eigenvalue Analysis in Terms of Deviations from Equilibrium. The procedure for linearizing the master equation for a second-order bimolecular reaction has been discussed by Snider and Ross.43 We now outline the steps taken to linearize the master equation for ( R l ) about equilibrium, and we discuss the eigenvalue-eigenvector solutions to the linearized equations. We define
ni = nc,Uj nj = nB,uj
n: = nC,,jc, i = 1, ..., 74
(4.8)
nie = nB,,;t, i = 75, ..., 108
(4.9)
=
rC11
n109e
~ I I= O
[Brl
~ I I =O [Brie ~
n109
=
ic11e
where the nC’s and [ ]e(s are concentrations at chemical equilibrium. The infinitesimal displacements from equilibrium are defined by (4.12) 5,” = ni - n?, i = 1 , ..., 110 We then substitute (4.8)-(4.12) into (2.3)-(2.6) and neglect all terms. Next we let wi = (i’/(n:)1/2
where we have defined xi,,
The eigenvectors
Ix,)
= Zi,,(ni?”*
(4.13)
The resulting coupled equations in vector notation become (d/dt)Jw) = - S ( W ) (4.14) where Iw) is a vector length 110 and S is a symmetric matrix of order 110. Since S is Hermitean and the displaced system always returns to equilibrium, the eigenvalues of S are positive semide-
(4.20)
are normalized according to = 6,”
(X,I*lX”)
(4.21)
where 9 is a diagonal matrix whose elements are given by \kij = (n,”)-16ij
(4.22)
It is useful to introduce a short notation for sums of eigenvector components over all the rotational states of a given vibrational level. These sums are denoted as follows 26
=
To,‘
Cxi,,
(4.23)
i= 1
41
(4.24) 64
(4.25) l d
(4.26) (4.27)
=
(4.28)
108 Xi,, i=96
Summing (4.19) over all states of HC1 yields 3
[HCl] = [HCl],
+ Cb,,e-hJZ
Tu,,‘
(4.29)
u=o
(4.10) (4.1 1)
(4.19)
P
A similar sum gives [HBr] and (4.10)-(4.11) yield [Cl] and [Br]. We note that all chemical concentrations must deviate from their equilibrium values by the same amount because of the stoichiometry of ( R l ) . Hence we define
E‘ = [HBr] - [HBr], = [CI] - [CI], = -([HCl] - [HCl],) = -([BrI - [Brl,) (4.30) Since these relations must be satisfied for all initial conditions and for all times, we have 3
-k
T1,>
=
xl09,p
=
-cT~,,c = u=o
-XllO,p
(4.31)
for all nonzero eigenvalues, p = 1, ..., 107. Let us order the nonzero eigenvalues, which are positive, so that AI < A2 < ... < AIo7. If A, > 1, every term beyond the pth term in eq 4.19 has become negligible compared to the
Effect of Vibrational-Rotational Disequilibrium
The Journal of Physical Chemistry, Vol. 90, No. 12, 1986 2621
It0)= C(c,/A,)(l -e-*Wr) c
TABLE I. Typical Experimental Conditions (DF and FP) and Initial Conditions Used in This Work (I-V) ~~
~
~
Ptotal,
torr
PAr,
YcI/YBr
torr
DF' 0.3-10 0.3-10 FPb 5-10' 5-10' I I1 111 IV V
(4.45)
13 13 13 0.3 200
12.3 12.3 12.3 0.28
lO-'-2O 10"-0.05 0.11 lo4 lo-"
188.5
If the nonzero positive eigenvalues are ordered such that AI
Ycl
YBr
10-9-0.02 10-8-0.004 5.7 x 10-3 5.7 x 10-4 5.7 x 10-6 5.7 x 10-6 5.7 x 10-6
106-0.3 10-"-0.08 5.2 X 5.7 X 5.7 x 10-2 5.7 x 10-2 5.7 x 10-2
-k YBr
< A2 < ... < A107 and if A, > (A,+l -
10-'-10.3 10-"-0.08 5.8 X 5.8 X 5.7 x 10-2 5.7 x 10-2 5.7 x 10-2
AJI, terms before the pth will contribute negligibly to the rate. If a single eigenmode p dominates eq 4.45, then
YCl
'DF = discharge flow. *FP = flash photolysis.
pth term, and terms before the pth term contribute negligibly to the net rate. If a single eigenmode p dominates the sum, then eq (4.19) and (4.29) reduce to
+ bpe-AJXi,,
ni = nf
+ bpe-AH'CTOJ,C
(4.33)
d[HCl]/dt = -A,([HCl] - [HCl],)
(4.34)
On the other hand if we substitute (4.30) into (2.1) and neglect terms of the order of
we have
- d F / d t = lk~([HBrle+ [Clle) + k-~([HClle+ [Brle)lF (4.35) Equating (4.34) and (4.35), using (4.30), gives
kl = Ap/([HBrIe
+ [ a l e + ([HClIe + [Brle)/K~'l
+ IP)
(4.37)
where I p ) is a time-independent initial rate vector whose components are calculated from (2.3)-(2.6) by using concentrations a t t = 0. L is a nonsymmetric matrix with eigenvalues A, and eigenvectors I$,,) which satisfy the equations
4$,)
= A,I$,)
(4.38)
If we form a matrix of column eigenvectors @
=
I$"))
(l$l)?l$2)9
(4.39)
the inverse may be obtained by solving 9-19 = 1
by Gaussian elimination. If we have
($,,-I1
($u-Il$,)
(4.40)
denotes row u of the inverse =.%v
(4.41)
W e let
19) = Xa,(t)I$,)
(4.42)
P
and (4.43)
where the time-independent constants c,, are calculated from c,
= ($L1lP)
Using (4.38)-(4.41), we obtain the solution to (4.37) as
(d/dt)I€O) = C,e-*,'Wr)
(4.47)
into (2.1) and neglect terms of the order of
(4.44)
we have
+ k,[N] = RIo
(4.49)
Rlo = kl[Cl]o[HBr]o
(4.50)
dt0/dt where
The solution to (4.49) is given by
9 = Rlo/(kl[N])(l - g k l L N l t )
(4.51)
so that d[O/dt = Rloe-kl[N1'
(4.52)
From (4.47) and (4.52), we see that kl
(4.36)
Eigenvalues A, and eigenvectors lx,) which satisfy (4.15) and (4.16) were found by using the subroutine RS (a subprogram of the EISPACK74 System). D. Eigenvalue Analysis in Terms of Deviations about Initial Conditions. W e use eq 4.8-4.12, where the n:, [...I,, and :4 are replaced by ny, [..I0, and E?, to denote concentrations at t = 0. Substituting these into (2.3)-(2.6) and neglecting all terms results in the coupled equations (d/dt)lP) = -LIP)
(4.46)
On the other hand, if we substitute Eo = [HBr], - [HBr] = [Cl], - [Cl] = [HCl] = [Br] (4.48)
US0
so that
= (c,/AP)(1 - e - * P ? l $ , )
so that
(4.32)
3
[HCl] = [HCl],
15')
= A,/[NI
(4.53)
Eigenvalues A,, and eigenvalues I$,,) which satisfy (4.38) and (4.41) were found by using the subroutine RG (a subprogram of the EISPACK74system) for various initial Boltzmann distributions of HCl and HBr. The eigenvectors I$,) are normalized but are not orthogonal to one another.
V. Calculations, Results, and Discussion Calculations are performed for five sets of initial conditions which are typical of discharge-flow-type and/or flash-photolysis-type experiments. This is summarized in Table I which gives, for each of the sets of initial conditions, the total pressure and the partial pressure of Ar as well as two of the conserved quantities defined above. For each set of initial conditions, the phenomenological rate coefficients are obtained in the absence and in the presence of V-V,R,T and V-V,T energy transfer. For convenience we denote the collection of all V-V,R,T and V-V,T rate processes simply as V-V energy transfer. The calculations performed with these various initial conditions and with and without V-V energy transfer should serve to illustrate (i) the effects of the total pressure, (ii) the role of V-V energy transfer in equilibrating the vibrational distributions on the time scale of chemical reaction, (iii) the dependence, if any, of the phenomenological rate constants on chemical composition as determined by the ratio of initial atom concentration to initial molecule concentration, and (iv) the constancy or variation of phenomenological rate constants with time and thus the nature of the steady-kinetics regime for (Rl). We shall first present the numerical results obtained via the iterative technique and EPISODE. We shall also present eigenvalues and eigenvectors obtained by linearizing the master equation about chemical equilibrium, and these analytical results will be compared to the numerical results. Finally the 1 10-component system will be lumped into a 8-component one corresponding to one state each of C1 and Br and six vibrational states of HC1 and HBr. The master equation for the lumped system will be linearized about chemical equilibrium, and the resulting eigenvalues and eigenvectors will be used to predict the time evolution of the pheno(74) B. T. Smith, J. M. Boyle, J. J. Dongarra, B. S. Garbow, Y.Ikebe,N. C. Klema, and C. B. Moler, Matrix Eigensystem Routines-EISPACK Guide, Spinger-Verlag, Berlin, 1976.
2622
The Journal of Physical Chemistry, Vol. 90, No. 12, 1986
Lim and Truhlar
menological rate coefficients and the internal-state distributions of HC1 and HBr. The master equation for the lumped system will also be linearized about initial conditions, and the resulting eigenvalues and eigenvectors obtained by this procedure will also be discussed. A . Results Obtained by Numerical Integration and the Iterative Method. The results for initial conditions 111 (see Table I) obtained by using the iterative method and EPISODE are presented in Table 11. The different runs are arbitrarily labeled by numbers in order to facilitate reference to the results of a particular run. The second column of the table indicates whether the results were obtained with the iterative method or with EPISODE. We investigated whether the quasi-steady state that is reached for some final value of R1/R-,is independent of the initial distribution by using initial distributions corresponding to Boltzmann distributions at various temperatures. These different starting distributions (xHC''and xHB:) are given in columns 3 and 4. The initial ratio is of the forward to reverse phenomenological rates (Rlo/R-lo) given in column 5. For each quasi-steady state that was observed, the range of R l / R I for which the phenomenological rate constant holds steady is given in column 7. Columns 8 and 9 give the ranges of time and C1 mole fractions corresponding to these ranges of R , / R , . [Note: for Tables I1 and 111, mole fractions are computed for reactants and products only, Le., XCI
= [ClI/ [RJI
= fHC1l / [R,P1 In all cases in Table 11, XHBr, defined similarly, equals 9.998 X lo-' and XHcl= X,, = 1.000 X lo4.] The forward and reverse phenomenological rate coefficients k , and Ll,their deviations K from the local equilibrium values, as defined in eq 2.24, the reverse one-way flux coefficient kif, and its deviation from kleare tabulated in columns 10-14. The intervals tabulated for a given quasi-steady state correspond to periods during which kl and k-, are constant to two or more significant figures. Time intervals between the end of the quasi-steady state and the beginning of another correspond to periods during which k , and k-, are varying more than this. We did not tabulate klr since, according to eq 3.5, it is always equal to kle for the sets of state-to-state rate constants employed in this study. For convenience we have associated an index A, Bl, B,, ..., or C with each quasi-steady state. This index is given in column 6 . The rationale for the way these indices are assigned is given below. All runs using EPISODE proceed to equilibrium. For runs using the iterative method, we will state explicitly in the discussion below where each run is stopped. Run 1 uses the iterative method and 300 K Boltzmann initial distributions of the molecules. The results of run 1 show that the last transients decay very rapidly, and, while the back reaction is still negligible, the system reaches a first quasi-steady state where the phenomenological rate coefficients deviate negligibly from their local equilibrium values. As R,/R_,decreases beyond 4 X lo6, the system departs from the initial quasi-steady state and kl and k-, decrease. A second quasi-steady state where K = 0.3 is reached when R I / R Iis approximately 5 X lo2. The progress of this later steady-kinetic regime is followed in run 2 by using EPISODE starting with the final R , / R Ivalue of run 1. Run 2 shows that the system leaves the second quasi-steady state when Rl/R-,is equal to 12. Although it is not recorded in Table 11, the system proceeds to chemical equilibrium without reaching a third quasi-steady state. Run 3 is like run 1 except that we changed xHCyto correspond to a 350 K Boltzmann distribution. Run 3 shows basically the same results as runs 1 and 2; Le., we find two quasi-steady states with K = 1.0 and 0.3, respectively, but the second one does not persist all the way to equilibrium. Run 4 utilizes the iterative method with initial Boltzmann distributions at 1500 and 300 K for HC1 and HBr, respectively, and run 5 continues the integration with EPISODE starting with the final distributions of run 4 corresponding to R I / R 1of 2.1 X lo5.The first quasi-steady state of run 1 is missed in run 4, and the first quasi-steady state reached corresponds to the second quasi-steady state of runs 1 and 2 where K = 0.3 for 3.3 X lo5 3 Rl/R-I3 3.8 X lo4. As the system approaches equilibrium, XHCl
w C i
Effect of Vibrational-Rotational Disequilibrium a second quasi-steady state appears where the nonequilibrium factor K is 3.9 x In run 6,we start with Boltzmann distributions of HC1 and HBr at 300 and 350 K, respectively, using the iterative method, and run 7 continues the time evolution of the system up to equilibrium with EPISODE. We find the same three quasi-steady states as obtained in runs 1-5: a first quasi-steady state where K is essentially unity, a second where k l and k-l are a factor of 3 less than their equilibrium values, and a last one which is characterized by a large nonequilibrium effect ( K = 3.9 X lov3). In all cases, the reverse flux coefficient starts off at a value corresponding to an equilibrium distribution of HCl molecules a t a given initial temperature which is not necessarily the temperature of the study, 300 K. When the initial temperature is 300 K, it first increases then decreases until it attains its equilibrium value when R I / R - ,equals unity. We made studies for conditions I and I1 of Table I for only one set of initial distributions, namely, those for 300 K initial Boltzmann distributions of molecules. The results of these calculations are in Table I11 where column 2 gives the initial conditions and columns 1,3-8 and 11-15 correspond to columns 1, 2, and 5-14 of Table 11. Columns 9 and 10 of Table I11 give the range of mole fractions of C1 and HCl for each quasi-steady state. [In all cases 0.8 6 XHBrG 0.9998 and XB, = XHc1.1 In all the runs of Table I11 the iterative method is used initially to obtain the first quasi-steady state after which EPISODE continues the integration up to equilibrium. Runs 1 and 2 are duplicated in Table I11 to facilitate comparison of the different cases. For runs 8, 1, and 2, we observe two quasi-steady states, the first of which is characterized by K = 1.00-0.99 for both conditions I and 111, and the second of which is characterized by K = 0.044 for conditions I and by K = 0.31 for conditions 111. For run 9, we observe three quasi-steady states characterized by K = 1.00-0.99, 0.32, respectively. Comparing run 8 with runs 9, 1, and 7.2 X and 2, we find that the period during which the earliest quasisteady state exists increases with decreasing values of Ycl/ YB~. The runs in Table I11 as well as additional runs without V-V energy transfer also show that range of R I / R Ior t corresponding to quasi-steady states for which K 1 and the earliest quasi-steady states are often characterized by K = 1.0 to 2 significant figures), that the eigenvalue dominated by reaction is changing slightly during the course of the reaction. The observation of K = 1.000 is also consistent with the general argument (given in subsection C) that K must be less than or equal to unity. Thus linearization about initial conditions provides a significant qualitative improvement in satisfying this constraint. We also note that the smallest nonzero eigenvalue for conditions I1 in Table X gives rise to a nonequilibrium factor K = 0.017 which is close to the K value obtained by quasi-steady state B2 of run 11 in Table 111. Another interesting feature that is shown in Table X is that the eigenvalues A, ( p C 5) do not have any component of reaction and are associated only with vibrational relaxation since &, = 0 (g = 1, ..., 4). Column 8 in Table X shows that the relative percentage differences between A, and for p = 1-5 for conditions I and I1 are large except for p = 3 and p = 5 whereas for conditions I11 the relative percentage differences between A, and X, (p = 1, ..., 5) are all small. We note that irrespective of whether V-V energy transfer is included, the relative percentage difference for a given p decreases as Ycl/YBrdecreases. This is in accord with our earlier prediction in subsection C that vibrational re-
x,
Lim and Truhlar
2630 The Journal of Physical Chemistry, Vol. 90, No. 12, 1986
TABLE IX: Results Obtained from Reaction Simulation Based on Lumped Eigenvalues and Eigenvectors for Expansion about Equilibrium
condition
2,
s
I1
I11
K
.I
6.0"-1 .r5 7.04-6.0-' 3.O"-9.0" 8. lW5-2.O4 7~ ~ - 2 . 0 - ~ 4.0"-4.1-' 8.0-4-2 .O-*
I
RilR-i. 2.3'2-5.911 5. 12-1 .O 7.2I3-6.7l2 7.18-6.85
0.928-0.932 4.441-2 0.993-0.994 0.3 32-0.320 7.229-3 1.OO-0.996 3.8 57-'
I .03-1 .o 3.915-1.2'2 1.5'-1.0
k , , cm3 molecule-' s-'
k-', cm3 molecule-' s-'
( 7.05-7.08)-'2
(3.24-3.25)-23 1.549-24 3.47-23 ( I . 16-1.1 2)-23 2.522-25 (3.49-3.47)-23 1.346-25
3.375-13 7.55-12 (2.53-2.43)-12 5.49444 (7.60-7.57)-12 2.932-14
TABLE X: Results of Eigenvalue Analysis of Lumped System for the Master Equation Linearized about Initial Conditions
condition
kl,y, cm3 molecule-I
I
1 2 3 4 5
2.496-4 4.4244 8.4924 1.2325 1,8315
s-' 1.036-12 1.836-12 3.525-12 5.112-12 7.600-l2
I1
1 2 3 4 5
3.127-' 2.8934 5.9454 8.8764 1.8315
111
1 2 3 4 5
7.2442 2.7244 5.6654 8.4984 1.8315
s-'
k-,,p cm3 molecule-'
P - X,l/
4.757-24 8.430-24 1.618-23 2.347-23 3.489-23
1.363-' 2.416-I 4.638-1 6.726-' 1.000
k,lkLy 1.737* 9.802' 5.107' 3.521' 2.368'
1.298-" 1.201-'2 2.468-12 3.684-12 7.600-'2
5.957-25 5.512-24 1.133-23 1.691-23 3.489-23
1.708-2 1.580-I 3.247-' 4.847-1 1.000
1.3873 58 1.4992 5.4 7.294' 0.61 4.886' 158' 2.368' 0.22'
3.007-14 1.131-'2 2.351-12 3.527-12 7.600-12
1.380-25 5.191-24 1.079-23 1.619-23 3.489-23
3.956-3 1.488-' 3.094-I 4.641-1 1.000
5.9863 1.5922 7.655' 5.104' 2.368'
S-I
Ky
9 % +a+ 0 67 0 34 0 0.60 1.43a 0 3.8b -5.636-'
$1,"
4
$2,~
$3,~
$4,~
$5,~
$6,~
0 -7.071-' 7.071-' 0 0 1.312" 4.821-'* 5.499-' -5.499-' -4.445-' 4.445-' 4.446-' 3.137" 2.968-2 -2.968-' 3.700-' -8.146'' -1.807-' 6.245-' -7.112-' 2.675-' l.101-3 -l.101'3 -1.185-2 -3.045-' -2.472-' -1.178-4 4.431-' 1.205''
0
0 0 0 0 -5.635-'
7.071-' 0 -7.071-' 0 0 0 1.5284 5.921-12 5.042-' -5.042-' -4.958-' 4.958-' 1.716-2 -1.716-2 4.264-' 3.344" 3.896-' -8.160-' 3.444-3 2.309-' -3.444-' 6.610-' -6.825-' -2.094-' 4.597-l 1.038-I -4.788-2 -3.193-' -1.963-' -7.503-'
0 0 1X2 0 0 1.6 0 -5.635-'
0 0 0 0 -7.071-' 7.071-I -5,002-' 5.002-I 1.552" 6.061-" 4.998-' -4.998-' 3.924-' -8.161-' 4.237-' 3.373-6 1.486-2 -1.486-2 -2.133-' 6.656-' -6.784-' 2.262-' -4.009-' 4.009-' -5.187-2 -3.196-' -1.919-' -7.216-' 4.613-' 1.022-'
2.5 8.8-2
= X - 1.478 X lo6 corresponding to kl = 7.667 X IO-'' cm3 molecule-' s-' and k-, = 3.519 in Table I$-'& = i5 for case 11 in Table IX. 'X, = k4 for case I1 in Table IX.
laxation processes can compete more favorably with reaction as Ycl/YBr decreases. Hence, the eigenvalues not dominated by reaction have increasingly smaller components of reaction as Ycl/Y,, decreases. This accounts for the smaller percentage difference for conditions 111 of Table X compared to conditions I and I1 of the same table. We can also rationalize the observation of quasi-steady states B in Tables I1 and 111. We note that the relative percentage differences in column 8 of Table X for p = 4 are much larger than those for p = 3. This means that A4 must change drastically during the course of the reaction and thus it would be impossible to observe a quasi-steady state associated with time-invariant A4. On the other hand, A, changes less drastically and more slowly than A4. For these cases, one can obtain a quasi-steady state associated with nearly constant A3. In particular, note that the eigenvalue spectrum for A, in Table V shows that X4 and As are quite close together (IA, - h,l/h, X 100% = 0.14%) but are well < 4. Hence, one might have separated from A,,, p > 5 and A,, expected that the phenomenological rate constants may not be constant to 2 or more significant figures so that quasi-steady state A of runs 9 and 11 in Table I11 is not observed. However, the eigenvalue spectrum of A,, in Table V shows that lAs - &l/& X 100% = 52% so that As is well separated from A4. Furthermore Table X shows that As for conditions I1 is largely associated with reaction and gives rise to rate constants equal to the equilibrium ones. This explains why quasi-steady state A is observed for conditions I1 in Table 111. We note that the percentage difference IA, - A41/A5 X 100% decreases with increasing YcI/YB, implying that A4 will contribute to the net rate earlier for conditions I. Consequently, the period over which quasi-steady state of type A exists will decrease as Ycl/YB, increases. The results in Table 11, 111, and IX confirm this prediction. In summary, the eigenvalues and eigenvectors obtained by linearizing the master equation about initial conditions help to clarify the numerical findings in Tables I1 and 111 as well as the earlier eigenvalue-eigenvector results in Tables VI-IX. F. Comparison of Efficiencies of the Various Methods. The results in Table I1 and 111 show that we can obtain quasi-steady states of type A and some quasi-steady states of type B with the
X
cm3molecule-' s-'.
=
x4 for case I
iterative technique even though it is a nonstiff method. To evaluate the relative efficiencies of the various methods, we will compare computational times for various runs. All times quoted are central processor times in minutes on the chemistry department's Digital Electronic Corporation VAXl 1/780 computer with floating point accelerator. The computational times for runs 1 and 2 are both 360 min, making a total of 720 min to run the initial conditions of run 1 to equilibrium. Run 3, which uses EPISODEonly, takes a total computational time of 996 min of which about 400 min is spent to reach quasi-steady state A and about 120 min is needed to reach quasi-steady state B. In contrast, the iterative method takes only 90 min to reach the first quasi-steady state but 165 min to reach the second quasi-steady state of run 2. Thus, the optimum method is to reach quasi-steady state A with the iterative method, and then switch to EPISODE to calculate the time evolution the rest of the way to equilibrium. This procedure was used for runs 9-12 where the computational times taken are 681, 650,908, and 901 min, respectively. In contrast to the numerical integration and iterative methods it takes only 2 min to obtain eigenvalues and eigenvectors for the master equation linearized about chemical equilibrium. If we assume rotational equilibrium we can calculate lumped eigenvalues and eigenvectors for the expansion of the master equation about equilibrium and about initial conditions in 0.2 and 0.3 min, respectively. The nonequilibrium internal-state distributions of HBr and HCl and the phenomenological rate coefficients can be calculated from the eigenvalues and eigenvectors, and x,, of the lumped analysis in about 0.5 min for each set of conditions. Therefore the final methods considered here involving the assumption of rotational equilibrium and an eigenvalue-eigenvector analysis for a vibrational model of the system linearized about equilibrium or about the initial conditions are about 3 orders of magnitude more efficient than numerically integrating the master equation to equilibrium. G. Consequencesfor Experiments. Rate constants are usually measured under conditions where back and side reactions are negligible. In the present treatment termolecular recombination reactions are neglected, but back reaction is not. We estimate for the present reaction that back reaction and termolecular re-
x,
The Journal of Physical Chemistry, Vol. 90, No. 12, 1986 2631
Effect of Vibrational-Rotational Disequilibrium
Since k-lf # k-leexcept at equilibrium, this implies that k-l # k-lc and as a result of eq 2.2, kl # kle. If internal-state relaxation were much faster than reaction for a wide range of Ycl/YBr,then we would have local equilibrium over this whole range, and we would observe the equilibrium rate constant and hence the same rate constant, independent of Ycl/YBr, over this whole range. In our case, in contrast, even though the phenomenological rate constants are well approximated by their equilibrium values under typical experimental conditions, local equilibrium is not maintained in the product HCl and also in general, in the reactant HBr. The nascent internal-state distribution of HCl is highly non-Boltzmann and the rate at which it relaxes to a Boltzmann form depends on the relative rates of reaction and vibrational relaxation which, in turn, depend on the ratio Ycl/YBr. From the discussion in section IV, we find that this results in a decrease of the extent of the steady-kinetics regime as Ycl/YB,increases. This has two consequences: (i) For larger YCl/YBr,an experimentalist would have to probe the reaction in its very early stages in order to observe a quasi-steady state of type A. (ii) In the event that an experimentalist does not have enough time resolution to yield the quasi-steady state of type A, time-dependent rate coefficients may be observed. When this happens in practice it is customary to attribute this time dependence to the occurrence of other competing reactions due to a significant build up of products. An example of this is found in Clyne and Stedman’s study of the reaction H HCl.76 The results found here suggest an alternative and/or additional consideration since time dependence of the rate constants may just correspond to the transition period from a quasi-steady state of type A to one of type B. Rate constants are usually measured under conditions where one of the reactants, B say, is in large excess. Under conditions where back reaction is neglected, pseudo-first-order rate constants
are determined for various initial concentrations of B. A plot of the pseudo-first-order rate constants vs. [B] generally yields a straight line whose slope gives a second-order phenomenological rate constant. Since the phenomenological rate coefficients characterizing quasi-steady states of type A do not depend on Ycl/YBr, and therefore on [HBr],, we would still obtain straight-line plots of pseudo-first-order rate constants by plotting kl[HBrIo vs. [HBrIo. The rate constants obtained here for quasi-steady states of type B do depend on composition, although they are constant in a given run. Experimentally rate constants are usually found to be independent of composition within the estimated experimental errors, although often only a limited range of composition is studied for a given reaction. When a composition dependence is observed [for example, Ambidge et al.77found that the rate of the reaction H + HCl increases by a factor of about 2 when Ycl/YBrincreases from 0(10-2) to O(l)], it can usually be rationalized in terms of side reactions, which are harder to minimize when wide ranges of composition are considered. The present results show a composition dependence for quasi-steady states of type B even though side reactions are excluded. Rate constants are sometimes measured under near-equilibrium conditions. The results in the previous section show that the phenomenological rate constant measured near equilibrium will differ significantly from that obtained in the absence of back reaction. The former is associated with a quasi-steady state of type C whereas the latter is associated with a quasi-steady state of type A. There are several characteristics that distinguish the two kinds of quasi-steady states observed in the present system. First, the nonequilibrium factor K for a quasi-steady state of type C is much less than unity. Second, quasi-steady states of type A depend on the initial distributions of molecules and on the extent of reaction whereas those of type C are independent of these two factors. Third, the phenomenological rate constants in quasi-steady states of type A do not depend on V-V energy transfer and/or Ycl/YBr,but those in quasi-steady states of type C do, since they are associated largely with vibrational relaxation. Last of all, the time interval for which quasi-steady states of type A exist depends on Ycl/YBr,whereas the time interval for which quasi-steady states of type C exist is almost independent of Ycl/YB,. If we can generalize the present example, it appears that we can use equilibrium theories of chemical kinetics to calculate rate constants even for relatively fast reactions since nonequilibrium effects are found to be negligible when back reaction is negligible, at least under typical experimental conditions. Nevertheless, observable phenomenological rate constants may show significant deviations from local equilibrium rate constants even when back reaction is negligible under certain conditions. These conditions can be inferred from eq 5.12. We will not observe equilibrium cannot 1 ) be neglected compared rate constants if ( k - l f / k - I e ) ( R 1 / R to unity even though R I / R l
