J. Phys. Chem. 1994, 98, 6 5 3 0 4 5 3 6
6530
Comparison of Master Equation and Trajectory Simulation of the Relaxation of an Ensemble of Highly Vibrationally Excited Molecules Gyorgy Lendvay Central Research Institute for Chemistry, Hungarian Academy of Sciences, H-1525 Budapest, P.O. Box 17, Hungary
George C. Schatz’ Department of Chemistry, Northwestern University, Evanston, Illinois 60208-31 13 Received: January 25, 1994; In Final Form: May 9, 1994’
We present the results of master equation simulations of the collisional relaxation of highly excited CS2 in a CO bath, using a rate coefficient matrix that is derived from trajectory calculations simulating individual collisions. The energy distribution from the master equation solutions is compared with distributions derived from trajectory simulations, and the results are used to assess sensitivity of the energy relaxation process to details of the rate coefficient matrix. We find that excellent agreement between master equation and trajectory energy distributions is obtained when the rate coefficients are represented using a biexponential form as a function of the magnitude of the energy change This biexponential form consists of strong and weak collision components, with the average energy transfer ( AE)being mostly determined by the strong component. We find that it is very important to let the parameters in these two components vary with the energy in the excited molecule in order to obtain realistic energy distributions from the master equation. These parameters have a roughly linear dependence on energy, but if this is ignored, then the energy distribution can switch from unimodal to bimodal as the relaxation proceeds. Only unimodal distributions are found in the accurate master equation simulations. Master equation simulations which use rate coefficients based on the strong component of the energy transfer distribution (Le., a single rather than biexponential fupction with energy-dependent parameters) are in reasonable agreement with the biexponential results. Other single-exponential approximations are in poor agreement with accurate results.
la.
Introduction Collisional energy transfer from vibrationally highly excited moleculeshas received considerableinterest recently. The process is important in unimolecular reactions14 and in the dissipation of energy gained by a system from chemical reactions or by excitation from absorption of UV or IR radiation.s.6 The energy transfer has been studied in “chemical” and, more recently, in “physical”experiments.’ In the former, the competition between chemical reaction and energy loss from a vibrationally highly excited molecule is used to infer information about the average energy transferred per collision.1~*-11The rate coefficient for a thermal unimolecular reaction increaseswith increasingpressure (which corresponds to increasing collision rate), and the rate and shape of the “fall-off curve” are determined by the details of the “state-to-state” energy transfer rate. In related chemical and photochemical activation experiments a highly vibrationally excited molecule is formed in a chemical reaction or in an internal conversion process following UV light absorption, and it either reacts (resulting in product formation) or loses its excitation in collisions (leading to a stabilized molecule). The ratio of product formation to stabilization dependson pressure (i.e,, on the collision rate) or, more precisely, on the rate of energy transfer. Information concerning the “state-to-state” energy transfer rates is inferred from these chemical experiments through comparison of the measured rate constants with solutions to the master equation characterizing the system. In this equation, the energy transfer process is described by “state-to-state” energy transfer rate coefficients, k(E’,E), associated with going from level E to level E’. These rate coefficients are expressed as the product of an energy-independent collision rate, Z, and a probability factor, P(E’,E),that tells what the probability of an ~~
a Abstract
~~~
published in Advance ACS Abstracts, June 15, 1994.
E to E’transition is if a molecule with energy E undergoes one collision. .Typicallythe energy transfer probability distribution, P(E’,E), is approximated by simple analytical functions that depend on one parameter (or perhaps two) which is closely related to the average energy transferred per collision, (AE). This parameter is almost always assumed to be independent of initial energy. To determine this parameter, the master equation is solved to determine the pressure dependence of the thermal unimolecular rate coefficient or the product formation to stabilization ratio as a function of the parameter. The best fit to the experimental pressure dependence determines the parameter, and thus P(E’,E). However, the results are usually not sensitive to the functional form of the energy transfer probability distribution as long as the proper average energy transfer is used9 (exceptions in the case of supercollisions have been noted; see later). In the “physical” experiments another regime of the energy transfer process is ~tudied.~J2-~5 These experiments follow the relaxation of an initial distribution of excited molecules which is far from equilibrium. Thermalization of the ensemble results from collisions with molecules of the bath. The average energy and, sometimes, the width of the energy distribution are determined as a functionof time (which is convertedto a ‘collision number” scaleusing an assumed collision rate). In the experiments which provide the average energyvs collision number, the “average energy transferred per collision” is the change of the average energy in Yonecollision” (i.e., the ensemble average of (AE) characterizing the energy transfer probability distribution). Very recently, direct experiments have appeared which reveal more details of the energy di~tribution.’~J~ These are also analyzed in terms of the energy transfer probability densities (“stepsize densities”) using the master equation. As can be seen from the above discussion, most of the information derived from the experiments relies upon some
0022-3654/94/2098-6530~~4.50/0 Q 1994 American Chemical Society
Relaxation of Highly Excited Molecules assumed energy transfer probability distribution. In the analysis of thermal unimolecular and chemical activation experiments, four functional forms for P(E’,E) as a function of p’-4 have been used most commonly? exponential, Gaussian, stepladder, or Poisson. At the time these functional forms were originally suggested, there was no information on the possible shapes the distribution might have, so these functions were chosen on the basis of simple physical arguments. In particular, it has always been assumed that P(E’,E) falls off rapidly with E’- 4. While this assumption is probably correct in a coarse sense, there is still significant uncertainty in the selection of functional forms for P(E’,E). In addition, there is evidence from recent experiments that there can be a significant probability for transferring large amounts of energy in a single collision (so-called “supercollisions”l7). The chemical experiments in the group of OreflC20 show that as much as 30 kcal/mol can be transferred in one collision in systems where both the hot molecule and the collision partner are polyatomic,and the relaxation studies in G6ttingenI4 and in Ann Arbor16 with benzene and toluene and various collision partners can be successfully interpreted if it assumed that the energy transfer probability distribution has a long tail. Troe2’ and Bernshtein and Oref22 have studied the influence of these long tails on energy relaxation and unimolecular reaction rates, showing that a small fraction of collisions transferring large amounts of energy per collision can have a significant effect on overall rates. Classical trajectory methods have been used for several years tostudy the functionalform of P(E’,E). Perhaps themost detailed of the early work was that by Brown and Miller,23who studied the relaxation of H02 by a He bath. They found that P(E’,E) had a long tail and that the functionalform of P was well described by a double exponential (Le., the sum of two exponentials for both the “up” and “down” branches), with one of the exponent parameters being 10-20 times the other. They did not use this function to study master equation solutions. There have been several subsequent trajectory studies? stimulated in part by the appearance of the detailed energy relaxation studies, and the double-exponential form of P(E’,E) has occasionally been noted.24.25 In a few recent trajectory studies,2”30 the energy relaxation process was stimulated not by solving the master equation, but instead by using classical mechanics to determine the relaxation that occurs as the result of successive collisions that the excited molecule suffers with bath molecules. The results of these studies have generally been in good correspondence with experiment (within a factor of 2 for ( AE)usually), suggestingthat theenergy transfer distribution function (which was not explicitly determined in the trajectory calculations) is realistically described. In the present paper we consider one of these previously studied systems, namely CS2 relaxation in a CO bath, in more detail, with a determination of the distribution function P(E’,E) over the complete range of energies, 0-100 kcal/mol, and subsequent solution of the master equation todescribe relaxationof an initially excited molecule. Since we already have trajectory results for this relaxation,24it is possible for us todetermine what functional representation of P(E’,E) is sufficient to describe this relaxation using the master equation. In addition, we have tested a number of approximate forms for P(E’,E) to see what errors result from each approximation.
Methods We applied the single energy collision method developed earlier and described in our previous papers.26-30 According to this method, the energy transfer from a single molecule is studied by classical trajectory calculations. Collisions of a hot molecule with a specified energy are followed from a large initial separation until the partners separate so that the interaction is negligible again. The internal state of the hot molecule is chosen according
The Journal of Physical Chemistry, Vol. 98, No. 26, 1994 6531 100.0
80.0
3
.% E
f
80.0
40.0
L
a
20.0
0.0 0.0
20.0
40.0
80.0 80.0 Energy, kcaVmol
100.0
120.0
Figure 1. Energy transfer probabilitymatrix P(E’,E) versus final energy E’for the relaxation of vibrationallyhighly excited CSz in a bath of CO molecules. The individual histograms show the probabilities calculated from 10 OOO trajectories with each line on a logarithmicscale. Thecurves
are derived from the “fit-of-the-fit”described in the text. to a microcanonical sampling procedure. The collision partner is selected from a thermal ensemble. P(E’,E) is determined by running an ensemble of collisions of hot molecules having the same initial energy E and then binning the final energy E’for each trajectory in the ensemble. In the present case, we ran 10 000 collisions of hot CS2 with CO at energies E = 5 , 10, 15, ...,95 kcal/mol. The temperature of the bath was 2000 K. The physical quantity the classical trajectory calculations unequivocally determine is not the energy transfer probability but the energy transfer rate k(E’, E). The trajectory calculations, however, were performed with a well-defined total collision rate coefficient which is the hard-sphere collision rate corresponding tothemaximumimpact parameter (7.OA) usedinthecalculations. As a result, the division of k(E’,E)’s by the total collision rate will lead to unequivocal P(E’,E)’s. In making comparisonswith experiment, it may not be appropriate to use the hard-sphere expression for the overall collision rate.30 However, this simply leads to a renormalization of the inelastic part of P(E’,E) that is not significant to the present analysis as is amply discussed in ref 30. Other details of the calculations are given in our previous papers.2&30 Result9 and Discussion The Energy Transfer ProbabilityMatrix from TrajectoryData, The results of the trajectory calculations can be visualized by creating histograms corresponding to the energy transfer probability distributionat each initial energy in CS2. Thedistributions, P(E’,E), calculated with a bin width of 1 kcal/mol in the final energy, E‘, are shown in Figure 1 as a function of the final energy E’. The ordinate scale on the probability axis is logarithmic, and the histograms associated with different initial energies, E, in CS2 are shifted upward in proportion to E. The location of the elastic peak on the P(E‘,E) curves shows the initial energy. At high E, the distribution functions display a long tail both in the energy-decreasing and in the energy-increasing wings. The “down” wing seems to be longer at high excitation, showing that energy loss dominates the energy transfer. With decreasing E, the “down” wing shrinks and the “up” wing increases. The two wings appear to balance each other at about the E = 15 kcal/mol distribution (the average thermal energy in CS2 at 2000 K is 16 kcal/mol). At low excitation, the “up” collisions dominate. A basic requirement that the energy transfer probability function must fulfill is detailed balance. This relates the forward and reverse transition probabilities to ensure that in thermal
6532 The Journal of Physical Chemistry, Vol. 98, No. 26, 1994
Lendvay and Schatz
equilibriumtransitions between the different energy levels do not change the Boltzmann distribution. This implies
P(E’,E) g(E) exp(-E/RT) = P(E,E’) g(E’) exp(-E’/RT) (1) To see if the trajectory calculations obey detailed balance, we calculated the classical density of states
g(E) = E ” ’ / ( R T ) ~exp(-E/RT)
T
6.0
T
5.0 5
m
Y
(2)
as a function of energy and, using this, the ratios
x=P(E ’,E)/P(E,E f)
(3)
Y = g(E?/g(E) exp(-(E’- E ) / R T )
(4)
and
The ratio of x toy was generallyclose to unity and always between 0.2 and 5 except in cases when E and E’were so far away from each other that in one of the “boxes” there was only one or two (or none) collisions. Taking into account that the histograms are obtained from a limited number of collisions so, as a result, they contain fluctuations, this agreement seems to be encouraging. In order to use the distribution functions to model relaxation processes, however, it is desirable to have functions in which the statistical fluctuations are smoothed out. Such functions are also needed to obtain energy transfer probabilities for initial energiesthat have not been studied in the trajectory calculations. The results in Figure 1 indicate that none of the four functional formsused previously in modeling chemicalsystems (exponential, Gaussian, stepladder,or Poisson) can be used without modification to describe the distributionsprovided by the trajectory calculations. The simplest functional form that we found could describe the results in Figure 1 is the biexponential function. This result is consistent with the results of Brown and Miller,23as well as other studies.*4J6 Accordingly, we used the function P(z) = aoS(z)
+ a, exp(-z/aJ + a3exp(-z/a4)
(5)
where z = 1 4,to fit the calculated points, separately fitting the “up” and “down” wings of the distribution. The first term on the right-hand side describes elastic (z = 0) collisions and is included to make JP(z) dz = 1. To determine the parameters ai-4, Marquardt’s algorithmwas used,31 and Poisson error was assigned to the histogram points. Only the wings of the trajectory distribution (IAEl> 0.5 kcal/mol) are used in determining ~ 1 - 4 ; the elastic peak (IAEl C 0.5 kcal/mol) is used along with the normalization condition to determine ao. Using this procedure,we obtained fits to the individual P(E’,E) curves corresponding to each calculated initial energy. The resulting parameters are plotted as a function of the initial CSz energy in Figure 2 (the points with errors bars) for the “down” wing of the distribution. Because of statistical error, these sets of parameters do not form a “self-consistent” set that describes all E and and E’ and also obeys detailed balance. In order to obtain P(E’,E) for all E and E’we fit the energy dependence of the points in Figure 2 to straight lines as a function of E. This produced the lines shown in the figure. A similar picture can be obtained for the “up” wing; however, detailed balance will not necessarily be satisfied. One way to get around this problem would be to include detailed balance directly in the fitting procedure, but this turns out not to be feasible because the E grid in the trajectory calculations is too coarse (AE= 5 kcal/mol) to provide a statistically meaningful representation of P(E’,E) for a given E’. As an alternative, we only did fitting for the “down” wing and then calculated elements of the “up” wing by requiring detailed balance. Adjustment of the diagonal elements P(E,E) was also required to preserve probability conservation. The
0.01 0.0
=e=
-=a= ”
20.0
”
”
”
40.0 60.0 80.0 Initial energy, kcal/mol
1 1.0
T
i
0.00
0.0
20.0
x
m ~
m , -
I and Less Than (
