J. Phys. Chem. 1993,97, 2251-2256
2251
Quantum Effects in Elementary Surface Reactions: CO Diffusion on Ni( 111) August Calboun and Douglas Doren' Department of Chemistry and Biochemistry, University of Delaware, Newark, Delaware 1971 6 Received: August 4, 1992; In Final Form: October 29, 1992
Quantum and classical calculations of equilibrium rate constants for surface diffusion have been carried out using a model of CO adsorbed on Ni( 111). Quantum transition-state theory (TST) rates were calculated using a Monte Carlo path integral method. Arrhenius plots of both the quantum and classical rate constants are linear from 100 to 300 K. Over this range, the quantum activation energy is slightly lower than the classical value due to a difference between zero-point energies in the reactant and transition states. At temperatures below 100 K, quantum rates are significantly higher than expected from the high-temperature Arrhenius plot. This is the result of quantum tunneling effects. The density matrix elements used to derive the quantum TST rate are also used to illuminate the differences between classical and quantum mechanics for this model.
Introduction The diffusion of molecules adsorbed on a solid surface is an essential elementary step in many surface reactions. A growing body of experimental work has explored molecular surface diffusion rates. Recent classical molecular dynamics calculations have illuminated the microscopic dynamicsof diffusion in a model of CO/Ni( 11l).' Such simulationscomplement experimentsby providing a molecular scaleview of moleculesurface interactions. The classical simulations show strong correlations near the transition state between the frustrated rotational motion of CO (Le., bending of the molecular axis relative to the surface normal) and the translation of the CO center of mass across the surface.' It is easy to understand why this might happen. The CO molecule binds to the surface through the carbon atom, with the molecular axis normal to the surface. During the diffusion process, the molecular center of mass moves laterally between favorable binding sites and the barrier to diffusion results from less favorable binding at intermediate positions. As molecules cross the transition state, they simultaneously execute thermal bending motion. Those that are midway between sites are equally likely to have their axes directed toward the site they are leaving or toward the adjacent site. However, molecules with axes bent toward the new site are much more likely to continue moving toward that site, since they have weakened their bond to the old site and may have already started forming a bond to the new site. Molecules bent back toward the old site are more likely to return to that site. Thus, the molecules that successfully cross the barrier to lateral translation all tend to be bent in the same direction at the transition state. In this work, we explore the effects of quantum mechanics on this model of molecular surface diffusion. Simple arguments suggest that quantum effects may be significant. The bending mode for CO on a transition metal typically has a frequency of 200-300 cm-I. At many temperatures of interest, the spacing between bending levels is not small compared to keT (kBT = 300 cm-1 when T = 430 K). Unless several quantum levels are populated,classical mechanics may give a misleading description of the thermodynamicand dynamical propertiesof the model. At very low temperatures, tunneling can also make important contributions to diffusion rates. As the temperature decreases, the rate of activated barrier crossing decreasesexponentially,but the contribution of tunneling remains nearly constant. An estimate of the temperature at which tunneling dominates the rate is given by the "crossover temperature", TO= h w ~ / 2 u k (in e the weakdissipation limit), whereqis thecharacteristic frequency of the barrier.* For the present model, we have estimated' that the barrier frequency along the reaction coordinate on the static 0022-3654/93/2097-2251$04.00/0
surface corresponds to 176 cm-I, so that Tois about 40 K. This estimate is crude, but it indicates that tunneling is not likely to be important at temperatures of 175 K and above, where the classical simulations were performed.' However, this estimate is also a reminder that tunneling may be importantat temperatures far from 0 K, even for a molecule as heavy as CO. To understand more accurately the influence of quantum mechanics on the diffusion rate in this model, we have calculated quantum transition-state theory (TST) rate constants for a restricted version of the model of CO/Ni( 111) used in the earlier classical simulations. These rates were determined by Monte Carlo evaluation of Feynman path integrals. We find that quantum mechanics has a substantial influence on diffusion rates at temperatures below 100 K. This is a result of both the large spacing between energy levels and tunneling effects. At first thought, one might conclude that a TST calculation would provide no insight into dynamical questions such as coupling between degrees of freedom. However, in the course of this calculation, we determined diagonalelements of the density matrix or, equivalently,the probability to observe the system in a given configuration. We will show that quantum effects are more evident in the density than in the rates. At temperatures of 200 K and above, the classical and quantum densities are very similar. If classical mechanics can reproduce thermodynamic properties like the density, it suggests (though it is not a conclusive proof) that classical dynamics will also be qualitatively correct. In contrast, at temperatures below 100 K, classical mechanics does not even qualitatively reproduce the density, so classical dynamics is not likely to be correct in this case.
Model and Method The model of CO/Ni(l 11) used here is based on a potential that has been described in detail Classical simulations have shown that lateral translational motion is much more strongly correlated with bending motion than with any other degree of freedom. We want to determine the effect of this coupling on quantum rates. To simplify the calculations, we have restricted all coordinates to fixed values except a lateral center-of-masstranslation coordinate, labeled x, and the bending coordinate,0, that specifiesthe angle of the molecular axis relative to the surface normal. In this restricted model, the surface atoms are fixed in their equilibrium positions and the position of the molecular center of mass above the surface is fixed at the value that minimizes the energy on the most stable (atop) site. One lateral component of the center-of-mass position is fixed so that x = 0 corresponds to an atop site. Varying x corresponds to translation of the molecular center of mass over a bridge site to Q 1993 American Chemical Society
2252 The Journal of Physical Chemistry, Vol. 97, No. IO, 1993
Calhoun and Doren rate calculations. The reaction coordinatewill be taken as x , the lateral position of the center of mass. More complicated functions of x and 8 could be used, but this simple choice will illustrate the points we wish to make. A dividing surface that separates the reactant region from the product region is defined by a fixed value of the reaction coordinate, x = x*. The saddle point of this potential is at the bridge site with 8 = 0, so x* is chosen to correspond to the bridge site. TST estimates the rate of reaction as the equilibrium density of molecules moving toward products at the transition state times their mean velocity,314
where g(x) dx is proportional to the probability of finding the molecule between x and x + dx. The integral in the denominator of eq 1 is over values of x in the reactantregion. The average velocity in a classical ensemble at thermal equilibrium is Figure 1. Potential energy of the model as a function of lateral position of the center of mass ( x ) and angle (in deg) of the molecular axis with respect to the surface normal (e). More than a complete unit cell is shown in the x coordinate.
-27
where p is the reduced mass associated with the reaction coordinate. Classically, the position distribution function for a system of N degrees of freedom may be defined as
-28 n
0 Q,
/
E -29 > a
.
2.
/
)r
F Q)
-30
C
W
-31 -32 0
-0.2 -0.4 -0.6 -0.8
x
-1
-1.2 -1.4
(4
Figure 2. Cross sections of the potential at fixed values of 0. The x position of the potential maximum depends on the angle. Negativeangles indicate that the molecular axis is directed toward the initial (reactant) site; positive angles indicate the axis is directed toward the final (product) state. The average bending angle in classical simulations is 8 O .
an adjacent atop site. The potential is periodic, with a period of 2.48 A. The CO bond length is constant, and the azimuthal angle of the CO axis is fixed so that bending motion lies within the plane normal to the surface that contains the x axis. Figure 1 shows how the model potential depends on x and 8. On the lowest energy path from reactants to products, the molecular axis is normal to the surface at the bridge site ( x = 1.24 A). However, thermal bending motion makes it unlikely that the molecule will cross the barrier in this minimum energy configuration. As described in the Introduction, bent molecules may experience forces driving them toward (away from) the next site even before (after) their center of mass has crossed the midpoint between sites. This is illustrated in Figure 2, which shows cross sections of the model potential for fixed values of the bending angle. Near the transition state, the direction of the force on the center of mass clearly depends on the bending angle. Even for this restricted problem, a full quantum dynamical calculation for this problem would be difficult. However, we shall show that quantum TST provides substantial information on the role of quantum effects. We begin with a review of classical TST to stress the similaritiesbetween the classical and quantum
where q represents all the coordinates beside the reaction coordinate, pq and px are respectively the momenta associated with the coordinates q and x, pi is the reduced mass associated with the ith degree of freedom, and h is Planck's constant. Normalization of g(x) is arbitrary, but the choice in eq 3 is analogous to that used in the quantum version of this theory presented below. With this choice, the denominator of eq 1 is simply the reactant partition function. Equation 3 shows that the density is essentially determined by the Boltzmann-weighted average potential. Thus, the same classical distribution g(x) would arise in a one-dimensional system with an effective potential given by the potential of mean force, Vmdx), defined by
where G is an arbitrary normalization constant and VQmf determines the zero energy for the potential of mean force. This effective potential is simply the free-energy change associated with a change of position along the reaction coordinate. To establish the connection with quantum TST, we introduce the classical density,
which gives the probability to find the system with coordinates q, x. Then g(x) is simply the integral of pCl over q and
The quantum analogue to the TST rate is easily expressed in terms of the quantum density, p Q The probability to find the system with coordinates q, x is the Boltzmann weighted sum of
Quantum Effects in Surface Reactions
The Journal of Physical Chemistry, Vol. 97, No. 10, 1993 2253
the probabilities for each eigenstate:
function of the Fourier coefficients,
(7) where the E[s and q [ s are the eigenvalues and eigenfunctions of the Schrbdinger equation. This equation offers a very indirect route to the equilibrium density, because a great deal of statespecific information (i.e., all of the eigenstates of the system) must be known in order to obtain a thermally averaged property. A much more direct approach is based on the observation that the right-hand side of eq 7 is an element of the density matrix in the coordinate representation:
where@= l / k ~ T .Thequantumdensitymatrixcan beefficiently computed using the path integral formulation developed by 3 ~ path integral formalism is Feynman in the 1 9 4 0 ~ . ~ The equivalent to traditional quantum mechanics but describes matrix elements as a sum over all possible paths between the states involved. For example, for the 2 degree of freedom model described above, the diagonal density matrix elements are given by a two-dimensional path integral
so it is a simple matter to constrain the centroid to a particular value. Note that only the odd index terms of eq 1 la contribute to the path centroid. If all of the previous discussion in terms of thedensity is restated in termsof thecentroid density, thequantum TST of Gillan and Voth et al. results. In the following, we shall calculate only centroid densities and PQ(x,q)should be understood to denote the density calculated with path centroids constrained to the point (x,q). The path integral approach is efficient because the integrals can be evaluated with Monte Carlo methods. Computational effort can be expended primarily on paths with significantthermal weight. To apply Monte Carlo methods, it is convenient to evaluate the ratio of the density to the free particle density. In terms of the Fourier coefficients, this ratio is
Here x ( t ) and d ( t ) are time-dependent functions that define a path and S is the imaginary time action integral,
where where I = pr2 is the moment of inertia of the CO molecule. Note that the rotational kinetic energy term in eq 10 does not have the correct form for the path integral formulation of a rigid rotator.' The differences between this action and the true rigid rotor action are small when bending is restricted to a narrow angular range. Still, to be rigorous, our model should be described as two coupled oscillators, rather than an oscillator coupled to a constrained rotor. For diagonal elements of the density matrix, only cyclic paths need to be included in the sum over paths. That is, all of the paths that contribute to PQ(Xo960) begin at a point (xo,Od at r = 0 and return to the same point at a later time, t = @h. It is natural to expand such cyclic paths in Fourier series? x ( t ) = xo
krt + p a k sin( -)
@h
k=l
(1 l a )
The kinetic energy along any path can easily be written in terms of the Fourier coefficients, and the integration over paths can be expressed as an integral over values of the Fourier coefficients. Gillan9 and Voth et a1.I0 have demonstrated that a quantum analogue to eq 6 should not use the diagonal density matrix elements as we have defined them. Instead of evaluating PQ at fixed values of the initial point, X O , the centroid of the path, f
Fso
= 1
'h
x ( t ) dt
should be constrained. Near the transition state, this *centroid density", PQ(.%,q), can be qualitatively different from PQ(x0,q) when tunnelingcontributionsdominatetherate.9J0 In the Fourier path integral approach, the centroid of each path is an algebraic
In evaluating these integrals, the Fourier series was truncated after 25 terms, making use of the fact that high-frequency contributions to paths tend to cancel one another. The length scales associated with the Fourier coefficients, defined in eq 15, become shorter as the index increasesso that these terms represent short-range fluctuations. In fact, retaining 25 terms in the series is probably more than enough. Calculations performed with only 15 terms in the Fourier expansion gave very similar densities, even at the lowest temperaturestudied. Thedifferences in density from these two calculations have no measurable effect on the derived potentials of mean force. Two thousand sets of Fourier coefficients were chosen using the Gaussian weight functions defined by the kinetic energy terms for importance sampling. Five runs of 2000 Monte Carlo stepswith different random number seeds showed that sampling of paths was acceptable; distribution functions and potentials of mean force derived from the larger sample are indistinguishable from those derived from a single run at 200 K. As we shall show below, at 50 K,all 5 runs of 2000 steps are needed to accurately determine the potential of mean force near the transition state. We want to emphasize that we have used Monte Carlo integration to perform the path integrals, but not to integrate over x or 8. Instead, we have calculated the centroid density as a function of x and 8 and then used traditional numerical integration to evaluate the integrals in the quantum analogue of eq 6. As we shall illustrate below, the information available in the coordinate dependence of the centroid density is useful for understanding the role of quantum effects.
Results and Discussion The quantum centroid density and classical density at 50 K are shown in Figures 3 and 4. The quantum distribution is much broader than the classical one in the angular coordinate, but both
Calhoun and Doren
2254 The Journal of Physical Chemistry, Vol. 97, No. 10, 1993
-1 5-
1
I
-1 0-
0.8
-5v)
a,
EtT)
a, z a
0.6
0-
0.4
510-
0.2
1520-
I
I
I
I
I
-0.2 -0.4 -0.6 -0.8 -1 x
1
I
0 -1.4 -1.2
-1
-1.2 -1.4
-0.8 -0.B -0.4 -0.2
x
(4
Figure5. Quantum (filled squares) and classical (tilledcircles)distribution functions, g(x) and gcl(x), at 200 K. The symbols are larger than the la error bars. The curves are guides to the eye.
(4
Figure 3. Quantum centroid density, pg(x,B), at 50 K.
-
-1 5
-l0I
> T=50K
Classical Density (T=50 K)
-
lo] 15 I
I
I
I
I
I
-1.4 -1.2
-0.2 -0.4 -0.6 -0.8 -1 x
-1.2 -1.4
(4
Figure 4. Classical density, pc1(x,8), at 50 K. The density is calculated exactly at each point. The curves are jagged because the function is sampled only at the same finite set of points used to evaluate the quantum density.
distributions extend over a similar range in the x direction. This behavior is easily understood by approximating the system near the potential minimum as two independent harmonic oscillators. The bending mode has a much higher frequency (330cm-1 in this model) than the lateral translation mode (30 cm-I). At 50 K, several translational states have significant population, but only the ground state of the bend mode is populated. At this temperature and below, the angular dependence of the quantum density is that of the ground state of the bend mode. The classical density becomes concentrated about the potential minimum as the temperature decreases. When the classical thermal energy is much smaller than the zero-point energy, the classical approximationconcentrates the density in a much smaller region than the quantum ground state. At 50 K,the thermal energy is greater than the zero-point energy of the frustrated translation but much less than that of the bend. Thus, the classical density is a bad approximation to the angular dependenceof the quantum density. At higher temperatures the differences between the classical and quantum densities are less obvious. At 200 K,for example, the densities have the same qualitative shape and the quantum density is only slightly more diffuse than the classical density
-1
-0.8 -0.6 -0.4 -0.2
x
0
(4
Figure 6. Quantum potentialsof mean force at 50 and 200 K. Note that the activation energy is nearly unchanged. Except as shown near the transition state at 50 K, the symbols are larger than the la error bars. The curves are guides to the eye.
(see below). More subtledifferences are revealed in the quantum and classical distribution functions, g(x), derived by integrating the density over 6 (eq 4). The distribution functions at 200 K are shown in Figure 5. At all temperatures, the maximum of the classical distribution function is displaced from the energy minimum, but the maximum of the quantum distribution is at the energy minimum. The reason for this is simply that the classical probability is concentrated near the classical turning points, while low-lying quantum states have peak probabilities near the energy minimum.” The temperature dependenceof the centroid potential of mean force is illustrated in Figure 6. As usual, the potential of mean force decreases with temperature, but in this case, the lower temperature potential has higher curvature near the transition state. This is a quantum effect since we do not observe similar changes in the corresponding classical potentials. It is probably due to tunneling, though we cannot prove that from our calculations. The classical and quantum TST rates were calculated from eq 6 using the respective densities. In both cases, the classical value of 8 was used. This gives the correct quantum rate in the hightemperature limit9and for a parabolic barrier above the crossover temperature.I0 The Arrhenius plot in Figure 7 shows that the
Quantum Effects in Surface Reactions
The Journal of Physical Chemistry, Vol. 97,No. 10, 1993 2255
25
-25 -26
20
3 - -27
3E
15
-28
x=x
Y
Y
w
>.
c -
F
-29
E
10
w -30
-31
5
-32 -20 -15 -10 -5 0
2
4
6
8
10
12
p (mole/ kcal) Figure 7. Arrhenius plots of the quantum (filled squares) and classical (filled circles) rate constants. The lines are fit to all of the classical points (the slope is 2.5 kcal/mol) and to the quantum points above 100 K (the slope is 2.2 kcal/mol). Except as shown for the quantum rates at 50 and 75 K, the symbols are larger than the l a error bars.
classical and quantum rates have slightly different activation energies at high temperatures, but the rates at low temperatures are very different. Quantum effects are dramatic at 50 and 75 K. The rates are much higher than expected by extrapolating from the hightemperature quantum rates. We attribute this to an enhanced probability to find the molecule at the transition state caused by tunneling into the classically forbidden region. The crossover temperature for this model (below which tunneling dominates thermally activated barrier crossing) was estimated above to be 40 K. The path integral calculation shows that tunneling is important even at 50 and 75 K. These temperatures appear to be above the deep tunneling regime since the rates are still temperature dependent. These results show that our crude estimate of the crossover temperature is a useful guide to the onset of tunneling effects. Previous work9-10suggests that the centroid density differs significantly from the density when tunneling is important. We have observed such differences near the transition state at 50 and 75 K by direct calculation of both the density and centroid density. At higher temperatures, the differences are negligible. The high-temperatureclassical and quantum activation energies (2.5 and 2.2 kcal/mol, respectively) differ by 0.3 kcal/mol. They differ because the zero-point energy in the reactant region is smaller than that at the transition state of the quantum system. This makes the energy barrier lower (relative to the reactant state) than in the classical system, causing a higher diffusion rate. To support this interpretation, we have examined sections of the potential energy surface at the energy minimum (x = 0) and the transition state (x = x * ) , as shown in Figure 8. The potential well is significantly wider (as a function of the bending angle, 0 ) at the transition state than at the energy minimum. By using a harmonic oscillator approximation to the potential to estimate zero-point energies, we find that the zero-point energy in the reactant region is 0.7 kcal/mol, and at the transition state, it is 0.4 kcal/mol. Thus, the difference of 0.3 kcal/mol between the classical and quantum activation barriers is accounted for by the change in zero-point energy between the reactant state and transition state. The density matrix is a full description of the quantum system at thermal equilibrium. It contains much more information than the TST rate. Having used the density as an intermediate quantity
0
5
10
15
20
8 (degrees) Figure 8. Cross sections of the model potential as a function of 8 at the minimum energy (atop) site and the transition state (bridge site). The curvature and frequency at the transition state are lower than at the minimum energy site.
20 18 16 14
E
12
4
10
0
w
Q
8
6 4
2 0 -20 -15 -10 -5
0
5
10 15 20
8 (degrees)
Figure 9. Quantum centroid density at 200 K as a function of bending angle at the energy minimum, x = 0. The curves are guides to the eye.
to obtain a rate, we now discuss some of the additional information that it contains about the effects of quantum mechanics in this model. Figures 9 and 10 compare cross sections of the classical and quantum densities as a function of B at x = 0 and x = x* at 200 K. As stated above, the classical and quantum results are similar. At x = x*, the twodensities have the same shape, differing by a scale factor. At x = 0, the classical and quantum densities have different angular dependence (recall that in a full threedimensional model, a Jacobian factor of sin B will be included in integrals over the density, so large values of 0 contribute more to the rate). The difference between the angular dependence of the classical and quantum density at the energy minimum reflects the fact that the bending frequency at the minimum is larger than keTat 200 K so that classical mechanics does not reproduce all the details of the density. At the transition state, the bending frequency is lower, as shown in Figure 8, and classical mechanics more accurately describes the density. The similarity of the classical and quantum densities at 200 K has important consequences for understanding the coupling between bending and diffusive motion. As described in the Introduction, these 2 degrees of freedom are strongly correlated near the transition state in classical simulations. We have shown that classical mechanics provides a good approximation to the density for in these degrees of freedom at 200 K and above. The
2256 The Journal of Physical Chemistry, Vol. 97, No. 10, 1993
1
1
0.05
The coupling between diffusion and bending that motivated this work was observed in classical simulations at 175 K and above. This coupling depends primarily on dynamics near the transition state. Thus, the present work provides preliminary evidence that a full quantum dynamics calculation at 200 K may be expected to show correlations between bending and lateral translation similar to those already seen in classical dynamics simulations.
0.04 0.035 *" Z
d
0.03 0.025
v
Q
Calhoun and Doren
0.02
0.01 5 0.01
0.005
0 -20 -15 -10 -5
0
5
10
15 20
0 (degrees) Figure 10. Quantum centroid density a t 200 K as a function of bending angle on the dividing surface, x = x*. The curves are guides to the eye.
bending frequency is small enough at the transition state for thermal averaging to remove the distinctions between classical and quantum equilibrium properties at these temperatures. At the stable binding site, where the bending frequency is higher, thermal averaging is not as effective, but classical mechanics is qualitatively correct. At 50 K, classical mechanics provides a bad approximation for equilibrium properties since the quantum and classical densities are qualitatively different functions of x and 8 (see Figures 3 and 4). Of course, equilibrium properties cannot be used to prove anything about dynamics. Nevertheless, if the temperature is high enough that classical predictions of the density are accurate, it supports the hypothesis that classical mechanics is a good approximation for dynamics as well. Thus, a classicaldescription may be adequate to understand dynamics near the transition state at 200 K, but at 50 K,classical mechanics cannot be correct.
Acknowledgment. We cannot fully acknowledge,and may not fully be aware of, Dudley Herschbach's influenceon our thinking. This project seemed initially to have little in common with Dudley's own interests, yet this paper echoes many essential themes of his work, particularly concerning the relationship of classical and quantum mechanics. Acknowledgmentis also made to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for the partial support of this research. The work on which this material is based was also supported by the National Science Foundation under Grant CHE-9015368. We thank the University of Delaware and its Scienceand Engineering Scholars program for financial support during the initial part of this work. References and Notes (1) Dobbs, K.D.; Doren, D. J. J . Chem. Phys. l992,97,3722;J. Chem. Phys., submitted. (2) HBnggi, P.; Talkner, P.; Borkovec, M. Rev. Mod. Phys. 1990,62, 251. (3) Doren, D. J.; Tully, J. C. Lungmuir 1988,4, 256;J . Chem. Phys. 1991, 94,8428. (4) Pechukas, P. In Dynamics of Molecular Collisions; Miller, W. H., Ed.; Plenum: New York, 1976;Part B. ( 5 ) Feynman, R.P.; Hibbs, A. R. Quantum Mechanicsandfath Integrals; McGraw-Hill: New York, 1965. (6) Feynman, R.P. Statistical Mechanics; Benjamin: New York, 1972. (7) Kleinert, H. Path Integrals in Quantum Mechanics, Statistics and Polymer Physics; World Scientific: New York, 1990. ( 8 ) Doll, J . D.; Freeman, D. L.;Beck, T. L. Ado. Chem. Phys. 1990,78, 61. (9) Gillan, M. J. Phys. Rev. Lett. 1987,58, 563;J . Phys. C 1987,20, 3621. (10) Voth, G.A.; Chandler, D.; Miller, W. H.J . Chem. Phys. 1989. 91, 7749. (1 1) Pauling L.;Wilson, E. B., Jr. Introduction to Quantum Mechanics; McGraw-Hill: New York, 1935.
