Full Iterative Solution of the Two-Dimensional Master Equation for

dimensional master equation, expressed in cintinuum rotation as. In eq 1, g(E,J) .... momentum elements at a fixed energy is a good compromise. The su...
0 downloads 0 Views 379KB Size
Download PDF
7090

J. Phys. Chem. 1996, 100, 7090-7096

Full Iterative Solution of the Two-Dimensional Master Equation for Thermal Unimolecular Reactions Stephen J. Jeffrey,† Kevin E. Gates,‡ and Sean C. Smith*,† Department of Mathematics and Department of Chemistry, UniVersity of Queensland, Brisbane, Qld 4072, Australia ReceiVed: NoVember 22, 1995; In Final Form: February 5, 1996X

The application of the Nesbet algorithm to the full iterative solution of the two-dimensional master equation for thermal unimolecular reactions is presented. This leads to an efficient algorithm for computing thermal rate coefficients in the falloff regime for reactions such as simple fission dissociations, radical-radical recombinations, and ion-molecule reactions wherein microcanonical dissociation rate coefficients are sensitive to both the vibrational energy and the angular momentum of the excited species.

I. Introduction It is well recognized that accurate prediction and modeling of rate coefficients for unimolecular dissociation, recombination, and collision-complex-forming bimolecular reactions require careful attention to the role of angular momentum in the kinetics of these elementary processes.1 The effect of angular momentum on the intrinsic (microcanonical) rate coefficients is most significant for potential surfaces where there is no barrier to the recombination process. Put simply, this is because the large change in geometry between the transition state and the reactant molecule (or collision complex) implies a corresponding large change of rotational energy. Hence, the dissociation rates depend on both rotational and vibrational energy, in contrast to dissociation over a chemical barrier where the rate depends primarily on the vibrational energy.2 These basic energetic effects, imposed by the conservation of angular momentum, have important consequences in relation to reaction rates and branching ratios, which can only be reliably computed by the development of sufficiently rigorous numerical methods. The separability of molecular interactions in the gas phase into unimolecular events and bimolecular events enables the overall thermal dissociation process to be modeled by the twodimensional master equation, expressed in cintinuum rotation as

dg(E,J) ) [M] ∫ ∫ dJ′ dE′ [R(E,J;E′J′)g(E′,J′) dt R(E′,J′;E,J)g(E,J)] - k(E,J)g(E,J) (1) In eq 1, g(E,J) is the statistical population of molecules with total energy E and total angular momentum J;R(E,J;E′,J′) is the bimolecular rate coefficient for collisional transition of molecules with initial energy and angular momentum (E′,J′) to final energy and angular momentum (E,J); k(E,J) is the total microcanonical rate coefficient for dissociation of the reactant molecule into all product channels. Computing the input quantities for eq 1 alone, Viz., the collisional energy and angular momentum transfer rate coefficients and the dissociation rate coefficients, presents a formidable problem. In recent years, significant progress has been made in improving the level of modeling and prediction of the collisional relaxation3-6 and reaction steps.7-11 †

Department of Chemistry. Department of Mathematics. X Abstract published in AdVance ACS Abstracts, April 1, 1996. ‡

0022-3654/96/20100-7090$12.00/0

In the present paper we describe a numerical algorithm for solution of the two-dimensional master equation itself, given the requisite input data. Our method extends the Nesbet algorithm,12 which was originally used by Gilbert and coworkers to solve the corresponding one-dimensional master equation,13 to the two-dimensional case. To the best of our knowledge, this represents the first general numerical solution of this problem that does not rely on any simplifying assumptions regarding the nature of R(E,J;E′,J′) or k(E,J). The program that we have developed runs comfortably on a workstation platform, typical compute times falling in the range of a few minutes for a rate coefficient at room temperature to a few tens of minutes for a rate coefficient at 1500 K. The outline of the paper is as follows. In section II we summarize previous work on the solution of the master equation, focusing on strategies that have been adopted for the twodimensional problem, and give a brief description of the Nesbet algorithm, which we have adapted for this purpose (further details of this method can be found in ref 13). In section III we discuss certain issues that arise in the solution of the twodimensional master equation that are different from its onedimensional equivalent, in particular, the formulation of detailed balance. In section IV, we present the results of illustrative falloff calculations that we have carried out for the methyl radical recombination reaction at 300 K and for the reverse process, dissociation of ethane, at 1500 K. These two examples serve to illustrate the main numerical performance features of the algorithm. II. Solution of the Master Equation Firstly, we deal with some relevant notation. For the purposes of modeling, it is convenient to factor the rate coefficient for collisional energy and angular momentum transfer into an overall rate coefficient for collisions between the molecule and bath gas (irrespective of the amount of energy transferred), z, and a collisional energy and angular momentum transfer probability distribution function (TPDF) P(E,J;E′,J′),

R(E′,J′;E,J) ) z P(E′,J′;E,J)

(2)

This factorization is to some extent artificial and requires a working definition of a “collision” for the purposes of calculating z.2 Further comment on this, in relation to the twodimensional master equation, is made below. Given the factorization of eq 2, the master equation becomes © 1996 American Chemical Society

Two-Dimensional Master Equation

J. Phys. Chem., Vol. 100, No. 17, 1996 7091

dg(E,J) ) ω ∫ ∫ dJ′ dE′ [P(E,J;E′J′)g(E′,J′) dt P(E′,J′;E,J)g(E,J)] - k(E,J)g(E,J) (3) where ω ) z[M] is the collision frequency with the bath gas. In discrete notation, eq 3 is written

dg ) Ug dt

(4)

where U is the (negative definite) matrix whose elements comprise the microcanonical rate processes of eq 3. After the relaxation of transient terms in the population during an induction period that is usually short on a chemically relevant time scale, the problem reduces to calculating the long-time population distribution of eq 4, which is the eigenvector of U with a corresponding eigenvalue closest to zero,

λ0g ) Ug

(5)

and the thermal unimolecular rate coefficient kuni is given by

∫ ∫ dE dJ k(E,J)g(E,J) ∫ ∫ dE dJ g(E,J)

kuni ) -λ0 )

(6)

The solution of the master equation, given requisite input data, is a difficult problem. A number of strategies of varying sophistication have been adopted for the approximate solution of this equation. The well-known strong-collision model14 represents perhaps the simplest of these, making the assumption that a single collision (statistically averaged) is sufficient to relax a molecule into an equilibrium distribution of energies and angular momenta. This assumption leads to a very simple solution for the statistical population g(E,J). The strong collision model provides a useful qualitative guide but is generally a poor approximation for quantitative purposes. Approximate analytic solutions to the master equation, which are valid in the lowpressure-limiting regime where the molecule dissociates immediately upon excitation above the reaction threshold, have been developed by Keck,15 Troe,16 and Borkovec and Berne.17 Troe has incorporated his low-pressure-limiting solution into a convenient approximate scheme for predicting falloff rate coefficient by interpolation between high- and low-pressure limits.18 The analytic solutions to the low-pressure-limiting master equation generally rely on either a diffusion approximation for the collisional transitions19-21 or an assumed exponential energy gap law for the transition probabilities.16 Numerical solutions to the one-dimensional master equation, formulated in terms of the vibrational energy alone, have been developed using modified power,22,18 Nesbet,13,22 and Davidson23 methods. Aside from these discretization methods, a stochastic approach has also been developed.24 Approximate numerical schemes for solving the two-dimensional master equation have in the past generally relied on making simplifying assumptions regarding the nature of the collisional energy and angular momentum TPDF’s that enable the two-dimensional problem to be reduced to one dimension. Included in this group are the fixed-V and fixed-J models of Penner and Forst,25 the weak-V/ strong-J relaxation model of Smith and Gilbert26 (weakcollisional vibrational energy relaxation/strong-collisional angular momentum relaxation), and a perturbative variant of this27 that approximately treated weak-V/weak-J relaxation. Recently, Wardlaw and co-workers28 have formulated a more generally valid approximation to the solution of the twodimensional master equation that is based upon transforming

the integral master equation, eq 3, into a corresponding differential equation, which is then solved by the use of a diffusion approximation. Green et al.29 have considered this model more generally in the context of a number of other diffusion approximations to the master equation. The method can be shown to be essentially exact for the case of a onedimensional master equation with an exponential down TPDF,

P(E;E′) ) N(E′) exp[-(E′ - E)/R], E′ g E

(7)

where N(E′) is the normalization factor and R is the average downward energy transfer. The accuracy of the method is related to the fact that a diffusion process gives the same longtime population distribution as the discrete jump process modeled by P(E;E′) if certain imbedding conditions are satisfied.29 In general, however, there appears to be no simple way to determine a diffusion model that corresponds to a given TPDF. Our algorithm for the solution of the full two-dimensional master equation, eq 3, is an extension of the Nesbet method utilized by Gilbert and co-workers13 for the one-dimensional case. The salient details are briefly summarized below. The matrix U of eq 4 is asymmetric but is readily symmetrized by the following similarity transformation,30

B ) SUS-1; Sij ) δijfi-1/2

(8)

In eq 8 and below, we use a multi-index i (or j, etc.) to represent indices for both the energy and the angular momentum of a given element in the two-dimensional grid. In eq 8, fi ) Fi exp(-Ei/kBT)/Q, where Q is the partition function for the molecular species and Fi the rovibrational density of states for a given energy and angular momentum.10a Defining the transformed eigenvector

c ) Sg

(9)

λ0c ) Bc

(10)

we have

The matrix B is real symmetric negative definite. The Nesbet algorithm perturbatively improves an initial guess for the eigenvector c(0) according to the prescription

c(n+1)i ) c(n)i + ri/(λ0 - Bii)

(11)

where λ0 is understood to be the present approximation to the smallest eigenvalue (i.e., based on c(n)) and ri is the ith element of the residual vector

r ) (B - λ0)c(n)

(12)

Note that there is some flexibility as to how frequently one updates the elements ci during one pass over the vector. Updating every element ci before moving on to the next is termed sequential relaxation, and updating all of the ci’s together at the end of a single pass is termed simultaneous relaxation. Sequential relaxation gives better convergence properties, but is computationally more intensive. For our model twodimensional calculations, we have found that updating the elements at the end of each (inner) loop over the angular momentum elements at a fixed energy is a good compromise. The success of the Nesbet algorithm is dependent on having an initial guess for the eigenvector that is reasonably good. This is provided by the population distribution predicted by the strong-collisional model.

7092 J. Phys. Chem., Vol. 100, No. 17, 1996

Jeffrey et al.

As noted by previous workers,13,23 there are certain superficial resemblances between the smallest eigenvalue problem of the master equation and the computation of the ground state electronic energy in quantum chemistry. Indeed, the Nesbet algorithm was originally derived for use in the latter context.12 More recently, the Davidson algorithm,31 a subspace expansion technique that uses the simultaneous Nesbet correction vectors to iteratively build up a basis set, has been shown to provide superior convergence properties for the quantum chemical problem.32 However, we have found that at lower temperatures, where the sought-after eigenvalue is extremely small, there are critical numerical instabilities associated with the formation of the n × n representation of B in the Davidson basis, and furthermore, the accurate diagonalization of this matrix, which prevent the Davidson method from being stable. The same difficulties arise in applying the Lanczos method to determine the smallest eigenvalue. The Nesbet algorithm, by virtue of being a single-vector iterative method, avoids these difficulties. There is considerable current interest in time-dependent aspects of the master equation, for example, with a view to calculating incubation times in shock tubes.33 Within this context we note that because it is a perturbative technique, the Nesbet algorithm is unsuitable for determining higher eigenvalues of the matrix B (or U) unless a suitably good initial approximation for the corresponding eigenvectors is available. The Lanczos method is likely to be more suitable for this purpose, and we are currently investigating this possibility. III. Detailed Balance and Related Issues In this section we discuss the implementation of detailed balance. Firstly, however, it is necessary to consider in more detail the choice of variables for formulating the master equation. In addition to the total rovibrational energy E and total angular momentum J, the projection of the angular momentum onto the spatially fixed z axis, M, is also conserved. Since the microcanonical dissociation rate coefficients k(E,J) do not depend on M, it is not necessary to explicitly include it in the master equation. However, it is informative to see how the more detailed master equation collapses to eq 3 and to note the implications of this with regard to formulating detailed balance. The master equation in terms of E,J and M is written

dg(E,J,M) ) dt

equilibrium is maintained). To avoid confusion, it should therefore be emphasized that as a consequence of eqs 13-15 the population g(E,J) includes the (2J + 1) degeneracy factor associated with the M manifold. Likewise, in the formulation of the detailed balance relations below, the density of states F(E,J) will include this (2J + 1) degeneracy factor. In our formulation of the two-dimensional master equation, instead of the total energy E, we utilize the variable

 ) E - CJ2

where C is the smallest rotational constant of the unimolecular reactant. This is essentially a convenience, since it allows one to take advantage more easily of the fact that, for dissociation channels proceeding over a barrier (i.e., a “tight” transition state), the microcanonical rate coefficient can be quite accurately formulated in terms of  with a threshold determined by the J-dependent effective potential.2 Indeed, the same simplification occurs with loose transition states at sufficiently high J values.10a  is commonly referred to as the “active” energy. We note at this point that the choice of variables outlined above assumes that, on a time scale defined by the average lifetime of the excited unimolecular species with respect to dissociation, energy is randomized among all rotational and vibrational degrees of freedom, subject only to conservation of E, J, and M. Hase34 has recently discussed the possibility that the projection K of the total angular momentum onto the bodyfixed z axis (the A principal axis) of a prolate or nearly prolate top might also act as a good quantum number. This will depend on the extent of Coriolis coupling present in the highly excited molecule, not only in the region of its equilibrium geometry but also as it distorts on the way to the transition state. The implementation of detailed balance for P(,′) in the solution of the one-dimensional master equation is straightforward and has been discussed previously [e.g., ref 2]. The standard practice is to specify probabilities for downward transitions and to obtain the probabilities for upward transitions by using the detailed balance relationship

P(;′)F(′) e-′/kBT ) P(′;)F() e-/kBT

P(′;) )

P(E′,J′,M′;E,J,M)g(E,J,M)] - k(E,J)g(E,J,M) (13)

g(E,J) ) ∫ dM g(E,J,M)

P (′;) ,  g ′ N()

∫ ∫ dM dM′ P(E,J,M;E′,J′,M′)g(E′,J′,M′) P(E,J;E′,J′) ) ∫ dM′ g(E′,J′,M′) ∫ dM′ P(E,J;E′,J′,M′)g(E′,J′,M′) ) (15) ∫ dM′ g(E′,J′,M′) In other words, as usual, one sums over the final states and averages over the initial states. In the absence of external fields, one can without loss of generality assume that g(E,J,M) is evenly distributed across the M manifold (i.e., a microcanonical

(18)

then

∫0 d′ P (′;) N() )  [1 - ∫ d′ P(′;)] max

(14)

(17)

Normalization is carried out iteratively from the top energy down. If the downward probability is specified by the unnormalized function P (′,), i.e.,

ω ∫ ∫ ∫dE dJ dM [P(E,J,M;E′,J′,M′)g(E′,J′,M′) -

where again we have used continuum notation. Integrating eq 13 over M and M′ leads back to eq 3 with the following definitions,

(16)

(19)

We now turn to the two-dimensional case. In previous work,16,25-27 for want of further information, P (′,J′;,J) was taken to be a separable function. For example, a separable exponential down model for collisional relaxation of both the active energy  and the minimum rotational energy R ) CJ2 has been utilized by Troe16 and by Smith et al.27 If the rovibrational density of states F(,J) is also approximated as a separable function,35,2 then P(′,J′;,J) becomes separable and the normalization procedure is carried out independently for  and J in a manner entirely analagous to that outlined above. The physical model, which generates a separable density of states F(,J), associates the total angular momentum J with the

Two-Dimensional Master Equation

J. Phys. Chem., Vol. 100, No. 17, 1996 7093

Figure 1. Schematic indication of the procedure for implementing detailed balance constants in two dimensions. Given an initial energy and angular momentum (,J), an unnormalized probability is integrated over all final states (′,J′) in the unshaded region. Normalized probabilities for transitions to final states in the shaded region are computed using the detailed balance relationship, eq 20, and previously normalized probabilities for the reverse transition. Normalization constants are recursively generated using eq 21, starting from the element in the top left-hand corner of the figure and moving sequentially across the J manifold for successively decreasing .

two-dimensional “external” rotations about the B and C principal axes and treats the rotation about the A principal axis of the unimolecular reactant (approximated as a prolate top) as an “internal” or “active” degree of freedom that makes no contribution to the total angular momentum. We have recently developed efficient methods for computing F(,J) that account correctly for the contribution of the rotation about the A principal axis to the total angular momentum.10a The resulting density of states F(,J) is seen to be nonseparable at low J values. The implication of this is that even if the unnormalized downward TPDF P (′,J′;,J) is assumed to be separable, the normalized TPDF P(′,J′;,J) will be nonseparable through the detailed balance relation -E′/kBT

P(,J;′J′)F(′,J′) e

-E/kBT

) P(′,J′;,J)F(,J) e

(20)

Thus, it is necessary to formulate detailed balance more generally for the nonseparable, two-dimensional case. We have implemented a “top down, left to right” scheme for imposing detailed balance, illustrated in Figure 1. The unnormalized TPDF P (′,J′;,J) is defined for all transitions with ′ <  and if ′ ) , for J′ g J. Probabilities for transitions where the final energy ′ >  and, if ′ ) , for J′ < J (the shaded region in Figure 1) are obtained through the detailed balance relation, eq 20. The normalization constants N(,J) are determined sequentially from the top down (outer loop over ) and left to right (inner loop over J) as indicated in Figure 1,

N(,J) )

∫0 d′ ∫0J dJ′ P (′,J′;,J){1 - δ(′,)[1 - h(J′,J)]}  J {1 - ∫ d′ ∫0 dJ′ P(′,J′;,J)[1 - δ(′,)h(J′,J)]} max

max

max

(21) where δ represents the Dirac delta function and h the Heaviside step function. Note that this more general prescription requires one to define the functional form of the TPDF for downward and upward transitions in J if ′ < . The unnormalized TPDF P (′,J′;,J) can be computed using well-established classical simulation techniques (see, e.g., refs

Figure 2. Convergence of the methyl radical recombination rate coefficient as a function of the number of Nesbet iterations. Convergence is expressed as the log of the percentage error, (kn - kconv) × 100/kconv, where kn is the computed rate coefficient after n iterations and kconv is the converged result. Calculations are carried out at T ) 300 K and a pressure of 1 × 10-4 Torr.

3-6). The reliability of such methods has been improved recently with progress on the problematic issue of zero-point energy conservation.36,37 For the present illustrative calculations, we have chosen an exponential gap model for the TPDF as follows:

P (′,J′;,J)) B(2J′ + 1) e-(-′)/R e-(R-R′)/γ, R > R′ -1

) B(2J′ + 1) e-(-′)/R e-(R′-R)[γ

+ (kBT)-1]

, R′ > R (22)

Note that the Jacobian B(2J′ + 1) is included to ensure that the normalization integral is the same in (,R) space and in (,J) space. Equation 22 incorporates an approximate reversibility relationship between upward and downward transitions in J via the rotational energy exponential gap model of Troe.16 As discussed above, this is not necessary, since eqs 20 and 21 will ensure overall detailed balance regardless of the functional form utilized in eq 22. However, since the density of states F(,J) is separable for medium to large J values, the separable functional form of eq 22 will imply a corresponding separability in P(′,J′;,J) for medium to large J values, and hence, independent detailed balance relationships for collisional relaxation of  and J would be implied. Such considerations will change, of course, if more detailed energy transfer studies reveal significant correlations between relaxations of  and J with an implied nonseparability of the function P (′,J′;,J). Finally in this section, we make an additional comment on the factorization of the bimolecular rate coefficient for energy and angular momentum transfer, R(′,J′;,J), into an effective collision rate coefficient z and the TPDF P(′,J′;,J), eq 2. It should be noted that the phenomenological cross section for rotation-translation (R-T) energy transfer in collisions is typically significantly greater than that for vibration-translation (V-T) energy transfer. Since a single phenomonological rate coefficient for collisions is factored out in eq 2, one should bear in mind the limitations of this convention2 carefully when interpreting the modelling of data by solution of the twodimensional master equation.

7094 J. Phys. Chem., Vol. 100, No. 17, 1996

Jeffrey et al.

IV. Calculations We have carried out test calculations for the recombination of methyl radicals in argon at 300 K and the reverse process, dissociation of ethane, at 1500 K. The microcanonical rate coefficients k(E,J) used in the calculations have been evaluated using the microcanonical variational Rice-Ramsperger-Kassel-Marcus (µVRRKM) algorithm developed by Smith.10 Input parameters describing a model interaction potential for the methyl/methyl system were taken from Wardlaw and Marcus.7 The high-pressure-limiting rate coefficient for recombination of methyl radicals calculated at 300 K with these parameters is 7.0 × 10-11 cm3 s-1.7,38,10 In the illustrative calculations of Figures 2-5, we have used an average downward energy transfer, R, of 200 cm-1 and an average downward rotational energy transfer, γ, of 200 cm-1. The bimolecular rate coefficient, z, for collisions of ethane with the bath gas, argon, was calculated in the standard manner using Lennard-Jones parameters.39 Grain sizes for convergence at 300 K were 100 cm-1 and 10 quanta for  and J, respectively. At 1500 K (Figure 6) the grain sizes were 200 cm-1 and 20 quanta. Figure 2 shows the convergence of the rate coefficient at 300 K and a pressure of 10-4 Torr. Convergence is achieved rapidly after only five iterations. This behavior is typical of the performance of the algorithm at the relatively low temperatures required for modeling atmospheric reactions (i.e., 300 K and below). This calculation takes approximately 2 min on an IBM RISC6000/370 workstation and requires around 5 megabytes of memory. By way of comparison of the relative amount of computation involved, it may be noted that the µVRRKM calculation requires around 15 min on the same machine. In Figure 3, the convergence of the eigenvector corresponding to λ0, from which one obtains the nonequilibrium population distribution, is illustrated by plotting the relatiVe residual as a function of the active energy  for an angular momentum J ) 20,

|

|

Figure 3. Plot of the relative residual (eq 23) as a function of energy for a fixed angular momentum J ) 20. The residual function is plotted for successively increasing numbers of Nesbet iterations, from 0 (i.e., the relative residual for the strong-collision population distribution) to convergence at 5. Calculations were carried out for T ) 300 K and a pressure of 1 × 10-4 Torr. For convenience, unity is added to the residual before taking the log so that the plotted function goes to zero for the very small residuals at low energies.

(∑Bijcj) - λ0ci

rreli )

j

ci

(23)

The threshold for dissociation at J ) 20 is ca. 304 (in units of 100 cm-1 grains). We have plotted the relative residual instead of the usual residual because the population of molecules at energies near and above the dissociation threshold is many orders of magnitude less than that at the average thermal energy. Hence, the usual residual is entirely insensitive to the convergence of the rate coefficient. The lines have been plotted successively for increasing numbers of Nesbet iterations, from zero (i.e., the starting vector obtained from the strong-collision model) to convergence at five iterations. Figure 4 shows a threedimensional profile of the converged nonequilibrium population distribution g(,J), plotted as the ratio of g(,J) to the Boltzmann population f(,J). The falloff of the population at successively lower energies with increasing angular momentum is a reflection of the threshold energy 0(J), which decreases approximately linearly with increasing rotational energy.16,26 In Figure 5 is plotted the falloff curve for the recombination rate coefficient as a function of pressure of the argon bath gas. In order to obtain the recombination rate coefficient as a function of pressure, we have used the well-established fact that the forward and reverse nonequilibrium rate coefficients vary in proportion to one another, with the proportionality constant being given to an excellent approximation by the equilibrium constant. This relationship has been derived by Troe from the pseudo-steady state approximation40 and by Smith et al.41 by formal solution

Figure 4. Three-dimensional profile of the converged nonequilibrium population distribution g(,J) at T ) 300 K and a pressure of 1 × 10-4 Torr, plotted as a function of the active energy  and angular momentum J. The function plotted is the ratio of the nonequilibrium population, g(E,J), to the Boltzmann population, f(,J).

of the multichannel recombination/chemical activation master equation. Included also in Figure 5 is the falloff curve obtained for the rate coefficient krecSC predicted by the simple strongcollision approximation. The low-pressure-limiting ratio βSC ) krec/krecSC obtained for the present energy transfer parameters (R ) γ ) 200 cm-1) is 0.24. Figure 6 shows the convergence of the rate coefficient for dissociation of ethane, calculated at 1500 K and a pressure of 1 Torr. Again, monotonic convergence of the lowest eigenvalue of the matrix is achieved. However, a significant increase in the number of iteration required occurs. The degradation in the performance of the Nesbet algorithm at higher temperatures can be attributed to two factors. (1) The starting vector, obtained by using the strong-collision approximation, is a much poorer

Two-Dimensional Master Equation

J. Phys. Chem., Vol. 100, No. 17, 1996 7095 Nesbet and Davidson algorithms has been well established (see, e.g., ref 42). Finally, we note that convergence with respect to grain size is more difficult to achieve at lower temperatures and that this requires a certain amount of testing to ensure accurate results for any given molecular system and set of reaction conditions. For the present illustrative calculations, grain sizes of 100 cm-1 for the energy and 10 quanta for the angular momentum sufficed to produce results with precision less than 2% at 300 K. The same precision was achieved at 1500 K with grain sizes of 200 cm-1 and 20 quanta, respectively. V. Conclusion

Figure 5. Falloff plot for the methyl radical recombination reaction at T ) 300 K. Parameters are as outlined in text. The solid line is the solution of the two-dimensional master equation with a weak-collisional TPDF. The broken line is the prediction of the strong-collision approximation.

Figure 6. Convergence of the unimolecular dissociation rate coefficient, kuni, as a function of the number of Nesbet iterations. The convergence is expressed as the log of the percentage error, (kn - kconv) × 100/kconv, where kn is the computed rate coefficient after n iterations and kconv is the converged result. Calculations are carried out at T ) 1500 K and a pressure of 1 Torr.

approximation to the true eigenvector than it is at lower temperatures. This can be seen by noting that the initial rate coefficient at zero iterations, corresponding to the strongcollision rate coefficient, is a factor of 258 times greater than the converged weak-collisional rate coefficient. (2) At higher temperatures, the matrix B is significantly less diagonally dominant because the TPDF becomes broader, that is, collisions of a molecule in a given (,J) element have a chance to access a wider range of final states (′,J′). Indeed, this factor has been underestimated in our illustrative calculations because we have maintained the average downward internal and rotational energy transfer parameters at 200 cm-1. Despite this, the change in the detailed balance relationship at the higher temperature still has the effect of broadening the TPDF. The importance of diagonal dominance in relation to the convergence behavior of

We have presented the first full solution of the twodimensional master equation for thermal unimolecular and recombination reactions in the gas phase. The method computes falloff rate coefficients directly without the need for interpolation schemes and requires no simplifying assumptions as to the nature of either the microcanonical dissociation rate coefficients or the collisional energy and angular momentum TPDF. Our algorithm follows that of Gilbert and co-workers,13 which was successfully utilized for the solution of the one-dimensional master equation. However, a number of features, notably the implementation of detailed balance constraints, require generalization and extension for the case of the two-dimensional master equation. The full solution is quite efficient and will easily facilitate the testing of simpler and faster approximations. Simpler approximations to the solution of the two-dimensional master equation16,28 will continue to play a crucial role in modeling and correlating experimental data, where typically many hundreds of calculations are required in order to establish reliable fitting parameters. The full solution to the two-dimensional master equation, although more numerically intensive than approximate models, provides the advantage of being able to produce the correct result for arbitrary TPDF’s. This generality will enable rigorous investigation of the consequences of features of collisional relaxation recently observed experimentally, such as the dependence of average collisional energy transfer parameters on total internal energy,43 supercollisions,44-46 and correlations between energy and angular momentum transfer.46 In addition, the effect of radiative emission47,48 could also be readily incorporated into the TPDF. The calculations presented in this paper represent primarily an advance in algorithms for the computational modeling of unimolecular reactions. The ability to readily perform such calculations will enable a clearer focus to be placed on the fundamental problems of computing the requisite input data for the master equation, Viz., rate coefficients for collisional energy and angular momentum transfer and the dissociation dynamics of excited unimolecular species. Acknowledgment. This work has been carried out with the assistance of the Australian Research Council Large Research Grants program (Grant No. A29532744). We are grateful to Dr. Neville Bofinger for helpful interaction during the early stages of this project. References and Notes (1) Aubanel, E. E.; Wardlaw, D. M.; Zhu, L.; Hase, W. L. Int. ReV. Phys. Chem. 1991, 10, 249. (2) Gilbert, R. G.; Smith, S. C. Theory of Unimolecular and Recombination Reactions; Blackwell Scientific Publications: Cambridge, MA, 1990. (3) Lim, K. F. J. Chem. Phys. 1994, 100, 7385.

7096 J. Phys. Chem., Vol. 100, No. 17, 1996 (4) Miklavc, A.; Markovic´, N.; Nymann, G.; Harb, V.; Nordholm, S. J. Chem. Phys. 1992, 97, 3348. (5) Clarke, D. L.; Oref, I.; Gilbert, R. G.; Lim, K. F. J. Chem. Phys. 1992, 96, 5983. (6) Lendvay, G.; Schatz, G. C. J. Chem. Phys. 1993, 98, 1034. (7) Wardlaw, D. M.; Marcus, R. A. AdV. Chem. Phys. 1988, 70, 231. (8) Klippenstein, S. J. J. Chem. Phys. 1992, 96, 367. (9) Klippenstein, S. J.; East, A. L. L.; Allen, W. D. J. Chem. Phys. 1994, 101, 9198. (10) (a) Smith, S. C. J. Phys. Chem. 1993, 97, 7034. (b) Smith, S. C. J. Phys. Chem. 1994, 98, 6496. (11) Miller, W. H.; Hernandez, R.; Handy, N. D.; Jayatilaka, D.; Willets, A. Chem. Phys. Lett. 1990, 172, 62. (12) Nesbet, R. J. Chem. Phys. 1965, 43, 311. (13) Gaynor, B. J.; Gilbert, R. G.; King, K. D. Chem. Phys. Lett. 1978, 55, 40. (14) Robinson, P. J.; Holbrook, K. A. Unimolecular Reactions; WileyInterscience: New York, 1972. (15) Keck, J. C. J. Chem. Phys. 1967, 46, 4211. (16) Troe, J. J. Chem. Phys. 1977, 66, 4745, 4758; Z. Phys. Chem. 1987, 154, 73. (17) Borkevic, M.; Berne, B. J. J. Chem. Phys. 1986, 84, 4327. (18) Troe, J. J. Phys. Chem. 1979, 83, 114. (19) Keck, J. C.; Carrier, G. J. Chem. Phys. 1965, 43, 2284. (20) Nikitin, E. E. Theory of Elementary Atomic and Molecular Processes in Gases; Clarendon: Oxford, 1974. (21) Troe, J. J. Chem. Phys. 1982, 77, 3485. (22) Gilbert, R. G.; Luther, K.; Troe, J. Ber. Bunsenges. Phys. Chem. 1983, 87, 169. (23) Schranz, H. W.; Nordholm, S. Chem. Phys. 1983, 74, 365. (24) Barker, J. R. Chem. Phys. 1983, 77, 301. (25) Pe´nner, A. P.; Forst, W. Chem. Phys. 1975, 11, 243; 1976; 13, 51. (26) Smith, S. C.; Gilbert, R. G. Int. J. Chem. Kinet. 1988, 20, 307, 979. (27) Smith, S. C.; McEwan, M. J.; Gilbert, R. G. J. Chem. Phys. 1989, 90, 1630.

Jeffrey et al. (28) Robertson, S. H.; Shushin, A. I.; Wardlaw, D. M. J. Chem. Phys. 1993, 98, 8673. (29) Green, N. J. B.; Robertson, S. H.; Pilling, M. J. J. Chem. Phys. 1994, 100, 5259. (30) Montroll, E. W.; Shuler, K. E. AdV. Chem. Phys. 1, 361. (31) Davidson, E. R. J. Comput. Phys. 1975, 17, 87. (32) Davidson, E. R. Comput. Phys. Commun. 1989, 53, 49. (33) Shi, J.; Barker, J. R. Int. J. Chem. Kinet. 1990, 22, 187. (34) Zhu, L.; Chen, W.; Hase, W. L.; Kaiser, E. W. J. Phys. Chem. 1993, 97, 311. (35) Forst, W. Theory of Unimolecular Reactions; Academic Press: New York, 1973. (36) Miller, W. H.; Hase, W. L.; Darling, C. L. J. Chem. Phys. 1989, 91, 2863. (37) Lim, K. F.; McCormack, D. A. J. Chem. Phys. 1995, 102, 1705. (38) Klippenstein, S. J.; Marcus, R. A. J. Phys. Chem. 1988, 92, 3105. (39) Hippler, H.; Wendelken, H. J.; Troe, J. J. Chem. Phys. 1983, 78, 6709. (40) Troe, J. Annu. ReV. Phys. Chem. 1978, 29, 223. (41) Smith, S. C.; McEwan, M. J.; Gilbert, R. G. J. Chem. Phys. 1989, 90, 4265. (42) Morgan, R. B.; Scott, D. S. SIAM J. Sci. Stat. Comput. 1986, 7, 817. (43) Brenner, J. D.; Erinjeri, J. P.; Barker, J. R. Chem. Phys. 1993, 175, 99. (44) Hassoon, S.; Oref, I.; Steel, C. J. Chem. Phys. 1988, 89, 1743. (45) Luther, K.; Reihs, K. B. Ber. Bunsenges. Phys. Chem. 1988, 92, 442. (46) (a) Mullin, A. S.; Park, J.; Chou, J. Z.; Flynn, G. W.; Weston, R. E., Jr. Chem. Phys. 1993, 175, 53. (b) Mullin, A. S.; Michaels, C. A.; Flynn, G. W. J. Chem. Phys. 1995, 102, 6032. (47) Barker, J. R. J. Phys. Chem. 1992, 96, 7361. (48) Herbst, E. Chem. Phys. 1982, 65, 185. Herbst, E.; Dunbar, R. C. Mon. Not. R. Astron. Soc. 1991, 253, 341.

JP953430Y