J . Phys. Chem. 1984, 88, 5076-5083
5076
TABLE IV: Apparent First-Order Photodissociation Rate Constants for Acetone in the Trophosphere Near Sea Level Compared with CH20 and CHICHO data rate constant, s-l solar zenith 106k10'k106k-
angle, deg 0
20 40 60 78 86
(CH3COCH3) 1.17 1.08 0.785 0.365 0.057 0.016
(CHzO)' 6.0
(CH,CHO)'
5.8
3.5 3.4
4.7
2.2
2.7
1.o
0.6
0.13
0.1
"These represent total decomposition rates which are the sum of several processes; data for CHzO are from ref 29-31; data for CH,CHO are from ref 35 and 36. photodecomposition in air shows a pressure-sensitive generation of methyl and acetyl radicals which, for a given pressure, is relatively independent of the band of excitation wavelengths (280-313 nm) and temperature (-1 to 28 "C). Estimation of the Theoretical First-Order Photolysis Rate Constants for Acetone in the Troposphere. All of the acetone photodecomposition quantum yields observed in this work can be accurately summarized through the empirical equation 47 derived
=
aCo1= 0.0766 + 0.09415e-[a"l/(3.222x~o'*) (47)
from a least-squares fit of the data. Relation 47 was selected for = a + b/(c its simplified formulation. The representation of d[air]), expected in theory from the proposed mechanism, also fits the data well but requires a greater number of constants. The concentration of air, [air], is in units of molecules cm+. Equation 47 can be combined with solar irradiance and acetone absorption coefficient estimates to calculate the apparent first-order photo-
+
chemical decay constants for various solar zenith angles and altitudes within the t r ~ p o s p h e r e . ~However, ~ we may use the representative theoretical solar irradiances for the lower troposphere, published by Demerjian et al.,34 and our measured absorption coefficient data for acetone, to estimate the first-order acetone photodissociation constant for various solar zenith angles at a location near sea level. These results are summarized in Table IV, where they can be compared with the photolysis constants for C H 2 0 and C H 3 C H 0 . It is apparent that the photodissociation of acetone in the lower troposphere is much smaller than that for C H 2 0 and CH3CH0. The photodissociation lifetimes of the three ) 40" solar zenith angle are 14.8 days for compounds ( l / k d i s Wfor CH3COCH3, 5.9 h for C H 2 0 , and 5.3 days for CH3CH0. Since the rate constant for H O attack on CH3COCH3is also much smaller than those for C H 2 0 and CH3CH0, it is readily seen why the tropospheric [CH3COCH3]can build up to significant levels. Acknowledgment. This work was supported by an Interagency Agreement (DW49903 19-01-1) between the Environmental Protection Agency and the National Science Foundation. We thank Dr. William R. Stockwell of NCAR for his help in computer simulation of the acetone photooxidation results and Professor J. Alistair Kerr of the University of Birmingham for his evaluation of rate data used in the simulations. Registry No. CH3COCH3,67-64-1; CH,OH, 67-56-1; CH,O, 5000-0; COZ, 124-38-9;CO, 630-08-0. (33) Dr. Robert Chatfield of the National Center for Atmospheric Re-
search is working with us in carrying out these calculations which will be oublished smaratelv. (34) K. K.Demirjian, K. L. Schere, and J. T. Peterson, Adv. Enuiron. Sci. Technol., 10, 369 (1980). (35) A. Horowitz, C. J. Kershner, and J. G. Calvert, J . Phys. Chem., 86, 3094 (1982). (36) A. Horowitz and J. G. Calvert, J . Phys. Chem. 86, 3105 (1982).
A Classical plus Tunneling Model for Unimolecular Reaction Dynamics: The HNC HCN Isomerization
-
Boyd A. Waite Department of Chemistry, United States Naval Academy, Annapolis, Maryland 21 401 (Received: March 12, 1984)
A classical plus tunneling model is developed for describing mode-specific behavior in unimolecular reaction dynamics. Initial conditionsare chosen to correspond to various bond-specific preparations, the system then evolving according to the full classical Hamiltonian. At barrier turning points, the trajectories are stopped, and local tunneling probabilities are computed and averaged over to yield a survival probability as a function of time. The HNC-HCN unimotecular isomerization reaction is studied with a potential energy surface of Pearson et al. Results indicate a fair amount of mode specificity at low energies, but much less mode specificity at energies above the classical threshold for the reaction. A phenomenological analysis of mode specificity in unimolecular reactions is also presented, yielding results in qualitative agreement with the trajectory calculations.
I. Introduction Various experimental and theoretical aspects of unimolecular reaction dynamics have recently been reviewed.'" Previous work by this author'7 and others8 has focused on the question of mode (1) K. A. Holbrook, Chem. SOC.Reu., 12, 163 (1983). (2) W. L. Hase in "Potential Energy Surfaces and Dynamics Calculations", ACS Symp. Ser., D. G.Truhlar, Ed., Plenum, New York, 1981, p 1. (3) C. E. Dykstra, Annu. Rev. Phys. Chem., 32, 25 (1981). (4) B. A. Waite and W. H. Miller, J . Chem. Phys., 73, 3713 (1980). (5) B. A. Waite and W. H. Miller, J . Chem. Phys., 74, 3910 (1981). (6) B. A. Waite and W. H. Miller, J . Chem. Phys., 76, 2412 (1982). (7) B. A. Waite, S. K. Gray, and W. H. Miller, J . Chem. Phys., 78, 259 (1 983).
specificity in unimolecular reaction dynamics. Central to this idea is the possibility of controlling the rate of unimolecular reaction by Preparing the reactant in a specific energy configuration, e& via overtone e x c i t a t i ~ nO ~f ~a' particular ~ bond or mode by using a finely tunable laser. For the process A*(no)
k(E,no)
products
(8) W. P. Reinhardt, Annu. Reu. Phys. Chem., 33, 223 (1982). (9) See, for example, "Advances in Laser Chemistry", A. H. Zewail, Ed., Springer, New York, 1978. (10) See, for example, "Laser-Induced Processes in Molecules", K. L. Kompa and S. D. Smith, Ed., Springer, New York, 1979.
This article not subject to U S . Copyright. Published 1984 by the American Chemical Society
HNC-HCN Isomerization
The Journal of Physical Chemistry, Vol. 88, No. 21, I984
5077
RRKM theory
Etot
Figure 1. A schematic rate constant profile for the metastable states in a unimolecular reaction. Points represent decay rates of actual meta-
stable states corresponding to various preparations of the parent molecule. The solid line corresponds to the RRKM prediction for the rate constant as a function of energy. Wide dispersion of points on this diagram correlates with mode-specific behavior, while coalescence of the metastable state decay rates to a monotonically increasing function of energy correlates with RRKM (Le., statistical) behavior. where A*(no) represents a polyatomic molecule prepared in a specific energy arrangement h,the system is mode specific if the rate constant k ( E , b ) is dependent on both the total energy E and the initial energy distribution, no, and it is statistical (Le., RRKM-like) if the rate constant k(E,h) depends only on the total energy E . Elegant experiments by Reddy and Berry have measured the mode-specific character for both methyl' and ally112isocyanide isomerization reactions. Preparation of the parent molecules in highly excited vibrational states via overtone excitation of the various types of C H stretches and other modes and the subsequent analysis leading to the unimolecular microcanonical rate constant k ( E , b ) resulted in plots similar to the one schematically depicted in Figure 1. Dispersion of the specific decay rate constants is indicative of mode-specific behavior, while coalescence of the individual decay rate constants toward the depicted RRKM curve is indicative of statistical behavior. These same individual decay rate constants are obtained theoretically by calculating inverse lifetimes of the metastable states which correspond to the bond-specific excitation of the parent molecule. Such calculations have been done for several model systems4g5as well as for a few simplified molecular system^.',^ The principal technique used is the complex scaling method'j which enables metastable (or resonance) state energies and inverse lifetimes to be calculated directly. These essentially exact quantum mechanical calculations reveal how the various types of intramolecular coupling affect the mode-specific character of the metastable states. For example, strong coupling between the bound degrees of freedom and the reactive mode leads to more statistical-like rate constant profiles, while weaker coupling enhances mode specificity. In addition to the potential energy surface effects, studies have also been made of the effects of symmetry on mode-specific b e h a ~ i o r . ~ J ~ Exact calculations have been limited to systems containing only two degrees of freedom and relatively few metastable states. Though there have been attempts to treat scattering systems containing many degrees of freedom by focusing on only the pertinent degrees of freedom,I5 it is unclear at present how such an approach might be applied to unimolecular decomposition or isomerization. A semiclassical branching model has been developed6 and applied directly to the problem of mode specificity in unimolecular reaction dynamics which can perhaps be extended
K. V. Reddy and M. J. Berry, Chem. Phys. Lett., 52, 111 (1977). (12) K. V. Reddy and M. J. Berry, Chem. Phys. Lett., 66, 223 (1979). (13) See, for example, Int. J . Quantum Chem., 14, No. 4 (lq78), the entire volume of which is devoted to complex scaling. Also, see ref 8. (14) W. H. Miller, J . Am. Chem. SOC.,105, 216 (1983). (15) W. H. Miller and S . D. Schwartz, J . Chem. Phys., 77, 2378 (1982). ( 1 1)
Reaction Coordinates
Figure 2. A schematic depiction of the probability branching model described in section 11. Only the reaction coordinate profile is shown, although the system possesses transverse degrees of freedom. Cumulative
survival probabilitiesare given inside the well, and individual tunneling probabilities for successive hits are given outside the well. to more realistic polyatomic systems. Previous work by Waite and Miller4 described a classical model for treating the question of mode specificity. It is based on the use of classical trajectories to describe the intramolecular energy transfer and incorporates semiclassical tunneling to describe the actual passage from reactants to products. The energy states of interest for such a model are indicated in Figure 1, i.e., those with energies near the top of the potential energy barrier where tunneling may be significant. Several very important unimolecular systems have been postulated to proceed at energies in this barrier region.I6 The advantages of this approach are twofold. First, classical mechanics is easily applied to relatively complex polyatomic systems. Second, by allowing for unimolecular reaction (Le,, passage of the system from reactants to products), the model extends the work currently being doneI7 on intramolecular energy transfer in nonlinear systems and provides a means for incorporating these effects in actual unimolecular systems. Previous worP has demonstrated that the classical plus tunneling model describes very well the unimolecular reaction dynamics for simple model systems. Section I1 describes in detail this classical plus tunneling model. The unimolecular reaction chosen for study is the isomerization of hydrogen isocyanide: HNC*
-
HCN
(1.2)
This system was chosen because of its relative simplicity in terms of number of atoms, as well as the availability of a meaningful potential energy surface due to Pearson, et a1.I8 A description of the potential energy surface is given in section 111. In addition, this system was chosen because of the possibility of important tunneling effects on the overall reacfion rate at low temperatures.16b Section IV presents the dynamical results, which indicate that mode specificity is significant especially at energies below the potential energy barrier maximum. A discussion of the role of the C N stretch, which is essentially adiabatic throughout the reaction and therefore causes the most striking manifestation of mode specificity, is also given. Section V presents a phenomenological model based on a master equation formalism for unimolecular reactions, with direct application to the question of mode specificity. This model provides an approach for determining the significance of intramolecular (16) (a) S. K. Gray, W. H. Miller, Y. Yamaguchi, and H. F. Schaefer, J . Am. Chem. SOC.,103, 1900 (1981); (b) S . K. Gray, W. H. Miller, Y . Yamaguchi, and H. F. Schaefer, J . Chem. Phys., 73, 2733 (1980). (17) See, for example, (a) S . A. Rice, ref 9a, p 2; (b) R. A. Marcus, D. W. Noid, and M. L. Koszykowski, ref 9a, p 298. (18) P. K. Pearson, H. F. Schaefer, and U. Wahlgren, J. Chem. Phys., 52, 350 (1975).
5078
The Journal of Physical Chemistry, Vol. 88, No. 21, 1984
energy transfer on the overall decay rate constants. Finally, section VI summarizes the results and presents a discussion of how they relate to the studies of stochastic behavior in classical mechanics.I7 11. The Classical plus Tunneling Model
The classical plus tunneling model as described below follows closely the discussion given previ~usly.~It has as its basis a probability branching analysis similar to that proposed by Wigner,19 as depicted in Figure 2. The basic physical idea of the model is very simple: A classical trajectory is begun inside the potential well (the initial conditions are chosen to correspond to the particular excitation under investigation). If the total energy is below the classical threshold, the trajectory will oscillate forever inside the well. If the total energy is above the classical threshold, the trajectory may eventually escape from the reactant well, depending on the initial distribution of energy within the molecule and the intramolecular dynamics.20 In either case, each time the trajectory experiences a classical turning point along the barrier direction, it is allowed to tunnel through it with a probability computed from the local properties of the trajectory at that time. The probability that by time t the unimolecular reaction has not occurred (Le., the survival probability) is
where P k is the tunneling probability for the kth barrier turning point, and K ( t ) is the number of turning points which have occurred by time t. If and when the trajectory crosses over the barrier (for the case of total energy greater than the classical threshold), the ”tunneling” probability for that “turning point” is set to one, resulting in a survival probability at all subsequent times equal to zero. The survival probability in eq 2.1 is for a single trajectory, and this must be averaged over the particular distribution of trajectories appropriate to the initial energy distribution being investigated. Such an averaged survival probability (Ps(t)) may decay exponentially, in which case the unimolecular rate constant k(E,no) is obtained as the negative slope of a plot of In (Ps(t)) vs. t. As shown in section V, however, this averaged survival probability may exhibit a more complicated decay behavior which in some manner is indicative of intramolecular energy exchange among the various vibrational modes of the molecule. It is computationally more convenient to calculate a cumulative for each trajectory, where Ptun(t) tunneling probability, Ptun(t), is related to the trajectory’s survival probability by Ptun(t) = 1 - Ps(t)
(2.2)
This cumulative tunneling probability is easily shown to be given by KV)
Pt,n(t) = C P k ( 1 -Pk-I)(l - Pk-Z)..*(l - P2)(1 k= 1
- PI)
(2.3)
This quantity is then averaged over appropriate initial conditions to give (Ptun(t)),and the averaged survival probability is then obtained from (ps(t)) = 1 - (Ptun(t))
.,
(2.4)
(a) Calculation of the Tunnelinp Probabilities. PG.The most \
,
I
. I
rigorous semiclassical approach for calculating tunneling probabilities is to integrate the classical equations of motion along complex time contours for which the particle tunnels through the barrier.21 This technique is very difficult to apply in practice:2 however, largely due to instabilities which arise in the propagation (19) J. 0. Hirschfelder and E. P. Wigner, J . Chem. Phys., 7, 616 (1939). (20) R.J. Wolf and W. L. Hase, J . Chem. Phys., 73, 3779 (1980). Also, see W. L. Hase and R. J. Wolf in “Potential Energy Surfaces and Dynamics Calculations”, ACS Symp. Ser.,D. G. Truhlar, Ed., Plenum, New York, 1981, n 37. r
(21) (a) T. F. George and W. H. Miller, J . Chem. Phys., 56,5722 (1972); (b) 57, 2458 (1972); (c) J. D. Doll, T. F. George, and W. H. Miller, J . Chem.
Phys., 58, 1343 (1973). (22) C. Cerjan, private communication.
Waite of the complex valued trajectories. Approximate versions of this rigorous approach based on vibrationally adiabatic motion during the actual tunneling are much simpler to implement, and, in the calculations which follow, the simplest of these approximations is employed, i.e., the vibrationally adiabatic zero curvature (VAZC) approximatio11.2~ For the particular application to mode specificity, it has been shown that the results are insensitive to the level of approximation used for calculating the tunneling probabilities, P k a 4 Within the VAZC approximation, the classical Hamiltonian used to describe the tunneling is given by
+
+ +
H@,,s,n) = p,2/2m Vo(s) (n 1/2)-w(s) (2.5) where s is the reaction coordinate (the steepest descent path from the saddle point), n is the vector of vibrational action variables for the transverse degrees of freedom, Vo(s) is the potential energy along the reaction coordinate, and w(s) is the vector of transverse vibrational frequencies which are functions of the reaction coordinate, thus maintaining the effects due to frequency modulation. It is essential to note that the classical trajectories for intramolecular motion prior to reaching the barrier turning point are computed with the full classical Hamiltonian, the approximate Hamiltonian of eq 2.5 being used only for the purpose of determining an approximate tunneling probability for the branching model. At time tk, Le., the time of the kth turning point for motion along the s direction, the current values of the vibrational action variables n are obtained. (Since the full Hamiltonian is nonseparable, the actions n must be approximated as the zeroth-order action variables pertinent to the Hamiltonian.) The tunneling probability is then given within the VAZC approximation by the semiclassical tunneling e x p r e s ~ i o n : ~ ~ P k = exp(-2rk)/[l exp(-2rk)] (2.6) where
+
rk =
1;
ds (2m[V0(s)+ (n
+1/2)4)
- E])’/2
(2.7)
s< and s, are the classical turning points inside and outside the potential barrier, respectively. (b) Initial Conditions. The averaging over initial conditions of the various trajectories must correspond to the particular excitation prokess under investigation. To specify the average over initial conditions corresponding to a given total energy E and initial distribution of energy among the transverse modes (defined by the initial zeroth-order action variables, no), it is convenient to specify the transverse coordinates and momenta in terms of their respective action-angle variables, which for the case of Morse oscillatorsz5 are given by xJ = x,”
+ l/a
In [1/X2(1 - (1 - X2)’l2)cos qJ] (2.8)
PI = 4mDA/kkf(axJ/aqJ)
(2.9)
where kM = 2(2mD)’/2/a, X = 1 - (2n, + l)/kM, nJ and qJ are the action-angle variables, and a and D are the Morse parameters. For state-specific excitation, n, is set initially to the given set of integers and the angle variables q are averaged over by standard Monte Carlo methods,z6 Le., each qJ is chosen randomly from a uniform distribution between 0 and 21r. The transverse coordinates and momenta are then calculated from eq 2.8 and 2.9. The reaction coordinate motion is always initiated at s = 0 (Le., the bottom of the reactant well), with the reaction coordinate momentum being determined by energy conservation: H(ps,s=O,nO)= E (2.10) The average over the phase of the reaction coordinate motion is (23) R. A. Marcus, J . Chem. Phys., 46, 959 (1967). (24) N. Froman and P. 0. Froman, “JWKB Approximation”, NorthHolland, Amsterdam, 1965. (25) C. C. Rankin and W. H. Miller, J . Chem. Phys., 55, 3150 (1971). (26) See, for example, R. N. Porter and L. M. Raff in “Modern Theoretical Chemistry”, Vol. 2, Part B, W. H. Miller, Ed., Plenum, New York, 1976, p 1.
The Journal of Physical Chemistry, Vol. 88, No. 21, I984
HNC-HCN Isomerization
5079
N
Figure 4. The coordinate system employed for the HNC-HCN unimo-
lecular isomerization.
As stated previously, the correctly averaged survival probability for this bond specific excitation is then given by (Psurv(t) )E,Q = - (ptun(t) )E,% (2.18)
time t
Figure 3. The function Ptu,,’(t),the cumulative tunneling probability for a single trajectory i, as a function of time. See eq 2.12. This is to be averaged over the time of the first half-oscillationin the reaction coordinate motion as described in section IIb. taken by averaging over the time T ~ for / motion ~ from the bottom of the well to the first turning point along the reaction coordinate direction (actually one-fourth of a full s oscillation). This is more clearly demonstrated by considering the function Ptu:(t), the cumulative tunneling probability for the ith trajectory whose initial conditions are given by q, = ti
(between 0 and 27)
no,= specified set of integers s=o
p s from eq 2.10
(2.1 1)
Let tk and Pi ( k = 1,2, ...) denote the times of the s turning points and the tunneling probabilities at these times, respectively. Figure 3 depicts a typical Ptu,’(t), an explicit expression for which is given bv m
Pt,,i(t) =
4Lh(t
-ti)
(2.12)
the exponential decay of which yields the state-specific rate constant k(E,n,,). Though this model was originally developed4 for treatment of classically forbidden processes, it has been easily extended to include energies higher than the classical threshold. Whenever a particular trajectory passes over the barrier, its survival probability for all subsequent times is simply set to zero. It has been shown4 that such a procedure leads to the standard methodology for determining purely classical unimolecular rate constants. A possible refinement of this extended model is to allow for semiclassical “reflection” for trajectories which are classically allowed to react. Such a refinement would probably not be significant, however. 111. Potential Energy Surface for HNC-HCN Isomerization The coordinate system used for the study of the unimolecular isomerization of H N C is shown in Figure 4. In this system, rl is the distance between the carbon and nitrogen nuclei, r2 is the distance from hydrogen to the center of mass of the C N fragment, and 0 is the angle between the vectors rl and r2. In these coordinates, the Hamiltonian takes the simple form H = P12/2P1 + PZ2/2PZ + be2/2>[1/PlrlZ + 1/Pzr221 + V(r,,rz34 (3.1) where p1 and p 2 are the conjugate momenta of the coordinates rl and r2, K~ and l 2are reduced masses given by
k=l
PI
where the quantity 42 is defined by 4ki
= Pk(1 - Pk-li)(l
- Pk-i)..,(l
- P2)(1 - PI) (2.13)
and h is the unit step function defined by for x > 0 h(x) = 1
=O
forx tki
t
tk
(2.16b)
so that the averaged tunneling probability is finally given by
= mCmN/[mC + mN1
(3.2a)
(3.2b) = m H h C + mN)/[mH + mc + mN1 and pe is the momentum (angular) corresponding to the reaction coordinate 0. For this Hamiltonian, it has been implicitly assumed that the total angular momentum J is identically zero. The results presented in section IV are based on J = 0. The potential energy function V(rl,rz,O) in eq 3.1 to be used in this study is an analytic function which was fit to ab initio points of the HNC-HCN surface calculated by Pearson et a1.18 The form of the function is given by the following: ~2
where the various Morse parameters are considered to be functions of the reaction coordinate 0 and are given by the simple cosine expansion 6
&N, DH-CN,aCN,(YH-CNlrle, rz? =
k=O
ck
COS
(ko) (3.4)
where the coefficients are chosen so as to give the best fit to the data points from the ab initio calculation. Values for these coefficients for the various Morse parameters are given in Table I. The potential energy profile, VO(O),along the reaction coordinate is depicted in Figure 5 . Notice that the barrier lies toward the H C N region of this isomerization reaction, and also the low frequency associated with the bending motion of the H N C side of the barrier. The form for the potential energy surface in eq 3.3 is very conducive to obtaining the zeroth-order action variables at the
Waite
5080 The Journal of Physical Chemistry, Vol. 88, No. 21, 1984 TABLE I: Coefficients for Potential Surface Morse Parameters, Eq 3.3 and 3.4 Darameter DHXN,
c1
Cl7
DCN,kcal/mol
212.0 91.67 2.39 1.776 1.174 1.372 27.2
kcal/mol
A-' A-' rlcI rZe, A ~ C N I
Vo(0),kcal/mol
C?
C7
0
0 0 0 0.0423 -0.0135 0.08 11 -0.579
9.67 0 0.2382 -0.0139 0.261 -18.78
cs
c4
0 0 0 -0.0093 0.0033 -0.009 2 -7.65
C6
0 0.333 0 0.0228 0.00031 0.0067 0.7649
0 0 0 -0.00098 0.0005 0.0026 0.935
0 -1.67 0 -0.0808 -0.0012 -0.030 -1.89
-.lo -.14 -.16 0
.1 .
.2
.4 2.5
1.3
.6
U .7
t (picoseconds) a.
8
(degrees)
Figure 5. The potential energy profile along the reaction coordinate 0 for the HNC-HCN unimolecular isomerization. 0 = 0 corresponds to the H C N arrangement, while 0 = 180° corresponds to the parent H N C arrangemerit. See ref 18. Also, Table I gives the parameters used for V,(O) as defined in eq 3.4. c. Y CI
turning point encounters. This is effected by simply assuming that, at the turning point, the rl, rz, and 0 motions are separable, using the Morse oscillator inversion formula for obtaining the action variables from the coordinate-momentum representationz5 ni = -yZ
+ (kM/2)11
a
n
-c
-.a
- (-p~'/~pLi- ~ M O ~ S ~ ( ~ I ) ) ' / ~(3.5) /~''*]
-1.0
where RM, p,, and D are defined as in eq 2.8 and 2.9. The tunneling integral of eq 2.7 is approximated as rk
-1.1
=
l/h
le> do ~ H - C N ( ~ ) ( ~ ~ H [ V+O ( ~ ) + O
