J . Phys. Chem. 1985,89, 2139-2144
2139
Quantum Mechanical Rate Constants via Path Integrals: Diffusion of Hydrogen Atoms on a W( 100) Surface Ralph Jaquet* and William H. Miller Department of Chemistry, and Materials and Molecular Research Division, Lawrence Berkeley Laboratory, University of California, Berkeley, California 94720 (Received: August 8, 1984)
Using Monte Carlo path integrals methods described recently by Miller, Schwartz, and Tromp [ J . Chem. Phys., 79,4889 (1983)], we have calculated quantum mechanical rate constants for diffusion of H atoms on a model W(100) surface. The most interesting aspect of the present work is that it includes the effect of surface atom motion, i.e., phonons, on the H atom tunneling dynamics. Although the description is fully quantum mechanical, the phonons enter the present path integral treatment in a way very similar to how they appear in the classical generalized Langevin models; the present formulation thus provides the correct quantum mechanical analogue to classical frictional effects. The principal qualitative result shown by the calculations is that surface atom motion increases the rate of H-atom tunneling from one stable site on the surface to another.
I. Introduction There is an ever-growing effort, both experimentally and theoretically, to illucidate the details of gassurface interaction dynamics.’ The theoretical problem of chemisorption on a solid surface has been approached by a variety of techniques to describe elastic, inelastic, and reactive processes, but the ability to model the reactive and diffusive dynamics on metal surfaces is severely hindered by the lack of knowledge of the correct interaction potential. In addition, approximations are necessary for the investigation of the dynamics because of the large number of degrees of freedom. Exact close-coupling calculations are possible on a rigid surface,’,* but in order to allow for energy exchange between the adsorbates and the solid surface the inclusion of surface atom motion is compulsory. In the case of a nonrigid surface the dynamics has usually been described by classical mechanics, often using response functions and generalized Langevin methods. Adelman and Doll,3 and Tully and co-workers: have formulated a generalized Langevin approach to gas-surface dynamics that derives from earlier work of Z ~ a n z i gMori,6 ,~ and Kubo’ on linear response theory. Several approximate forms for the friction and random forces have been proposed to model the response to the lattice.* We are interested in describing chemical reactions on surfaces, for example the recombination and desorption of hydrogen H surface’ ‘H
-
surface t ~p
(1.1)
and, in particular, wish to study the effects of surface atom motion, Le., phonons, on the H-atom dynamics. In this paper, though, we begin with the even simpler process of diffusion of a single hydrogen atom on a model of the W(100) surface. Since the dynamics is that of hydrogen atoms, we utilize quantum mechanical methods to determine the rate, which is especially necessary in this case because the diffusion is primarily by tunneling from one local potential minimum to another. To include many degrees of freedom quantum mechanically is not a trivial matter, and we have employed the path integral methods as recently formulated by Miller, Schwartz, and Tromp.9 One of the interesting features in the application of this approach to the present problem is that the phonon degrees of freedom, though treated fully quantum mechanically, enter in a way that is very similar to the classical generalized Langevin m0de1;~+~ one thus has the correct quantum analogue to the classical frictional effects that appear in the Langevin picture. The particular model system we have treated is a LEPS potential constructed by McCreery and Wolkenlo to describe H atoms on a W( 100) surface. Although new experimental infor-
* Permanent address: Lehrstuhl fur Theoretische Chemie, RuhrUniversitat Bochum, D-4630 Bochum, Federal Republic of Germany. 0022-3654/85/2089-2139$01.50/0
mation” indicates that this potential is actually not an accurate description of H/W-for one reason, because the surface structure changes with H-atom coverage-it nevertheless provides a useful model in order to study the effect of surface atom motion on H-atom diffusion. The principal qualitative result that our calculations show is that coupling to the phonons increases the rate of H-atom tunneling from one site to another. This can be understood by noting the generic form of the tunneling rate
k
(w/2*)P
(1.2)
where w is a frequency and P a tunneling probability. For the rigid surface, w is the frequency of H-atom vibration in its local minimum, and P is the probability of tunneling through the barrier from one local minimum to another. In the extreme limit of large-amplitude surface atom motion, the tunneling probability is so much larger when the surface atoms are close together that tunneling will occur essentially only from this configuration; The frequency w in this case is then the frequency of the (heavy) (1) (a) J. C. Tully, Annu. Rev. Phys. Chem., 31,319 (1980); (b) G. Ehrlich and K.Stolt, Annu. Rev. Phys. Chem., 31,603 (1980); (c) S. T. Ceyer and G. A. Somorjai, Annu. Rev. Phys. Chem., 28,477 (1977); (d) R. F. Madix and F. Bensiger, Annu. Reu. Phys. Chem., 29,285 (1978); (e) “Vibrational Spectroscopyof Adsorbates”, R. F. Willis, Ed., Springer-Verlag, Berlin, 1980; (f) “Vibrations at Surfaces”, R.Caudano, J. M. Gilles, and A. A. Lucas, Ed., Plenum Press, New York, 1982; (g) ‘Surface Physics of Materials”, J. M. Blakely, Ed., Academic Press, New York, 1975. (2) G. Wolken, Chem. Phys. Lett., 21, 373 (1973). (3) (a) S. A. Adelman and J. D. Doll, J . Chem. Phys., 61,4242 (1974); (b) S. A. Adelman and J. D. Doll, J . Chem. Phys., 63, 4908 (1975); (c) S. A. Adelman and J. D. Doll, J. Chem. Phys., 64, 2375 (1976); (d) S. A. Adelman and J. D.Doll, Acc. Chem. Res., 10, 378 (1977). (4) (a) M. Shugard, J. C. Tully, and A. Nitzan, J . Chem. Phys., 66, 2534 (1977); (b) J. C. Tully, G. N. Gilmer, and M. Shugard, J . Chem. Phys., 71, 1630 (1979). (c) J. C. Tully, J. Chem. Phys., 73, 1975 (1980). (5) R. Zwanzig, Annu. Rev. Phys. Chem., 16, 67 (1965). (6) H. Mori, Prog. Theor. Phys., 33,423 (1965). (7) R. Kubo, Rep. Prog. Phys., 29, 255 (1966). (8) (a) S. A. Adelman and B. J. Garrison, J . Chem. Phys., 65, 3751 (1976); (b) J. D. Doll and D. R.Dion, J . Chem. Phys., 65,3762 (1976); (c) B. J. Garrison, and S. A. Adelman, Surf.Sci., 66, 253 (1977); (d) A. C. Diebold and G. Wolken, Surf.Sci., 82, 245 (1979). (9) W. H. Miller, S. D. Schwartz, and J. W. Tromp, J. Chem. Phys., 79,
4889 (1983). (10) (a) J. H. McCreery and G. Wolken, J . Chem. Phys., 63,2340 (1975); (b) J. H. McCreery and G. Wolken, J. Chem. Phys., 63,4072 (1975); (c) J. H. McCreery and G. Wolken, J . Chem. Phys., 64, 2845 (1976); (d) J. H. McCreery and G. Wolken, J . Chem. Phys., 65, 1310 (1976); (e) J. H. McCreery and G. Wolken, J. Chem. Phys., 66, 2316 (1977); ( f ) J. H. McCreery and G. Wolken, J . Chem. Phys., 67, 2551 (1977); (g) G. Wolken and J. H. McCreery, Chem. Phys. Lett., 54, 35 (1978). (11) (a) L. D. Roelofs and P. F. Estrup, Surf.Sci., 125, 51 (1983); (b) Y. F. Chabal and A. J. Sievers, Phys. Rev. Lett., 44, 944 (1980); (c) R. F. Willis, Surf. Sci., 89, 457 (1979); (d) P. Thirey in ref If, p 231.
0 1985 American Chemical Society
2140 The Journal of Physical Chemistry, Vol. 89, No. 11, 1985
surface atom motion and P is the tunneling probability from this more favorable distorted surface atom geometry. Compared to the rigid surface, therefore, the effect of phonons is to decrease w and to increase P. Because P varies exponentially in the tunneling regime of interest in the present application, the increase in it overwhelms the decrease in w. This picture shows that the entropic and the enthalpic characters of the diffusion rate are both affected by coupling to phonons. Section I1 first summarizes the quantum mechanical methods used to carry out the present calculations. There are some minor (but significant) modifications of the earlier methodology9in order to make it more applicable to systems with many degrees of freedom. The application to diffusion of H atoms on W ( 100) is described in section 111, where the primary focus is on the qualitative effect that coupling to phonons has on the tunneling rates. Section IV concludes. 11. Path Integral Calculation of the Flux Correlation Function Here we review and slightly modify the method presented by Miller, Schwarz, and Trompg for determing Boltzmann-averaged quantum mechanical rate constants for a chemical reaction. The rate constant is given by the time integral of a flux-flux autocorrelation function
k = Z , - l I m0 d r Cdt)
[Fe-(B/z-il/h)HFe-(B/2+Jr/h)r
Ps2 H(pS,s) = 2m
(2.2)
“tr” denotes a quantum mechanical trace, p = (kBT)-I,H i s the full Hamiltonian of the system, and F is the symmetrized flux operator
+ Vo + AV(s)
where AV(s) is typically a potential barrier. When eq 2.3 is used for the flux operator, the imaginary time correlation function, eq 2.4, is given more explicitly by
where sN = s, = so = 0 after the differentiation; if the reactants and products are equivalent (Le., V(-s) = V ( s ) ) ,then this can be written more simply as
(2.6b) with so = s, = 0 after differentiation. The two Boltzmann operator matrix elements in eq 2.6 are now written out as a standard Feynman path integral,Is with the special feature that the same imaginary time increment AT I h@/Nis used in both path integrals, and for this to be possible the (imaginary) time T is chosen to be one of the discrete values T ,
(2.1)
where Z, is the partition for reactants, and the correlation function is defined by
cdt) = tr
Jaquet and Miller
Tn I
np( If. - 1) N
2
n = 0 , ..., N i.e., one of the intermediate time values; this leads to the following result for the path integral representation of the product of the two matrix elements: (,yNle-(B/2-T/h)HISn) (~,le-(B/2+r/h)H1sO) = I d s , ...I d s , , I d s , + ] . . .I d s N - ] ( sNle-ArH/h(sN-l) ...
...
(2.3)
( s,+, le-ArH/hls,) ( ~ , l e - ~ ~) ~( s1 / le-A7H/hlso) ~ l ~ , ~ (2.7)
where s is the coordinate normal to the “dividing surface” (s = 0 is the equation of the surface) through which the flux is computed. In order to be able to deal with systems with many degrees of freedom, the Boltzmann operator/time evolution operator e-@’/2*a/h)H is represented via a Feynman path integral expression and the path integral evaluated by a Monte Carlo random walk method. It is, in general, not feasible to do this for real values of the time 1, however, because the integrand of the path integral would be oscillatory. (Although see recent work by Do11I2 and Wolynes et aLI3which suggest otherwise.) We thus first calculate Cf for real values of T = it, Le., pure imaginary time
Equation 2.7 has a very useful structure: the two path integrals appear as one overall path integral from initial position so to final position sN with total (imaginary) time increment h@,the only modification being that the intermediate position s, at time T,, is not integrated over. One thinks of the matrix element (s,l exp[-((@/2) + ( ~ ~ / h ) ) f Z las l spropagation ~) from position so at time -hp/2 to position s, at time T , (time increment = T , + (h@/2)), and the other matrix element, (sNl exp[-(P/2) - ( T , / h))qIs,), as propagation from s, at time T , to sN at time +(h@/2) (time increment = (h/3/2) - 7,). The net time inverval for propagation from so to sN is h& We use the random walk algorithm described before9 to evaluate the two path integrals in eq 2.7, and the resulting expression for the complex time correlation function is (from eq 2.6b)
F = ‘/2[W@s/m) + @ , / m ) W l
Cf(T) = tr [Fe-(B/2-‘/h)HFe-(B/2+r/h)r
(2.4)
and then use these calculated values to generate a numerical analytic continuation to obtain C, for real r.l4J5 Since Cf is an even function of t , and thus of 7 , one uses the values computed for real T to construct a Pad6 approximant in the variable r2;T~ is then replaced by -t2 to obtain C d t ) . We note that similar expressions for rates in terms of flux correlation functions have been derived by Yamamoto16 and used by W~lynes.’~J’ A . One Degree of Freedom. We summarize the new, modified calculational procedure we use first for the case of only one degree of freedom. The Hamiltonian in this case is (12) J. D.Doll, preprint. (13) E. C. Bchrmann, G. A. Jongeward, and P. G. Wolynes, J . Chem. Phys., 79, 6277 (1983). (14) L. Schlessinger, Phys. Rev., 167, 1411 (1968). (15) G. A. Baker, Jr., “Essentials of Pad6 Approximants”, Academic Press, New York, 1975. (16) T. Yamamoto, J . Chem. Phys., 33, 281 (1960). (17) P. G.Wolynes, Phys. Rev. Lett., 47, 968 (1981).
(2.8) where the positions (si)are given in terms of the integration variables (w,S),i = 1, ..., n -1, n + 1, ..., N - 1 by 1/2
s. =
n-i+l
mN n - i + 1
z(w;)
(2.9a)
(18) R. P. Feynman and A. R. Hibbs, “Quantum Mechanics and Path Integrals”, McGraw-Hill, New York, 1965, p 63.
Diffusion of H on W ( 100)
The Journal of Physical Chemistry, Vol. 89, No. 11, 1985 2141
for i = 1, ..., n - 1, and
( N - i)si-l + SO si = N-i+ 1
+(mN
2ah2P
)
N -i ' I 2 z(wf) (2.9b) N-i+ 1
+
for i = n 1, ..., N - 1; so and s, are specified values, and the function z(w) is related to the inverse of the error function, an explicit representation for which was given before? The integrative over the variable Iw:) in eq 2.8 is done by straightforward Monte Carlo. The important feature of the above prescription, eq 2.8-2.9, is that the two path integrals, for the two matrix elements, are combined effectively into one path integral. This is even more important as we now consider the case of more than one degree of freedom. B. Two Degrees of Freedom. The case of two degrees of freedom, Le., one degree of freedom (P,Q) in addition to the reaction coordinate s, is characterized by the Hamiltonian
ps2 P2 H(p,,s,P,Q) = - + - + Vo 2m 2m
1 + -mWQ2Q2 + V(s,Q) 2
(2.10)
for i = 1, ...,N - 1. By combining the two path integrations into one overall Monte Carlo path integral, the Q degree of freedom thus contributes precisely the way it would is one were simply evaluating its partition function. Because of the coupling potential AV(s,Q), however, its contribution is not simply the harmonic partition function Ze. C. Two Degrees of Freedom Linearly Coupled to Many Harmonic Modes. The Hamiltonian is now generalized to
the same as before, eq 2.10, with the addition of harmonic degrees of freedom (qk]which are linearly coupled to the s and Q degrees of freedom. The generalization of eq 2.1 1 for the imaginary time function is
and the generalization of eq 2.6b and 2.7, for example, is
(") I ...
Cd7,) =
2
a2
-
2m asoas, dsl I d s , , I d s , + l . ..
IdsN-ll d Q o I d Q l...
With regard to the Q degree of freedom, combining the two path integrals into one overall path integral thus has the structure of a simpler trace over that degree of freedom. Incorporating the harmonic part of the potential in the Q degree of freedom into the importance sampling, as before? leads to the following expression (which is to be evaluated by Monte Carlo):
The special simplifying feature that arises here is that the can be carried out anaintegration over variables qo, ql, ..., lytically because these harmonic degrees of freedom are coupled only linearly to the s and Q degrees of freedom. With the expression in Feynmann and Hibbs,'* it is a straightforward calculation to show that the integration over the (qi]variables gives (2.17) where hU/2
A = C(Zmhwk)-' J h p 1 2 d T ~ ~ ud7'fk(7) ,2 fk(T') k
(coth (Uk/2) cash [Wk(T' - 7)]
+ sinh [Wk(7' - T)]) (2.18)
and Zk is the partition function for mode k z k
= [2 sinh (Uk/2)]-' uk = Auk,
with fk(7) = fk(s(T),Q(T))
(2.12) where ZQis the harmonic partition function for the Q degree of freedom
ZQ= [2 sinh(u/2)]-' u = hwQP
(2.13)
The coordinates (si)are given in terms of the Monte Carlo variables (w,S) as above, by eq 2.9, and the coordinates {ei),i = 0, ..., N 1, are given in terms of the Monte Carlo variables (wiQ]by =
Qi
= sinh
(
N-i+l
sinh ( u / N ) sinh sinh
(
where S ( T ) and Q(7) are the discretized paths (sO,s1,...,sN-1) and i = 0 , ..., N - 1, as Qo,Ql,...QN-1) at time T~ = h/3((i/N) given in terms of the Monte Carlo variables (wf)and {wiQ]via eq 2.9 and 2.14. It is not hard to show that A, of eq 2.18, can also be written in the more symmetrical form A = E(Zmhwk)-'
and finally, if the integrals over T and T' are evaluated by the trapezoid rule over the 71grid, then eq 2.20 becomes
+ u) (NU)
N-
N-i+ 1
u)
X
k
z(woQ) (2.14a)
mwQ tanh (u/2)
(
(2.19)
]
where z(w,Q) (2.14b)
2142 The Journal of Physical Chemistry, Vol. 89, No. 11, I985
The final expression for the correlation function, eq 2.16, is then
A)’ 2m
(y)4 m
2
4
-
112 asoasn
x
Jaquet and Miller
TABLE I: Sites of H on W(100) site
BE, eV
r,.“ an
ICN 2CN 5CN
3.00
3.12 3.70
2.74 1.72
Z.,b an 3.12
2.19
3.04
0.056
w.‘
cm-’
1232, 372 (2) 1232, 738, 3721’ 1232, 7381’ (2)
“Distance of H from the nearest neighbor W atom. bHeight of H above the 1CN plane of W atoms. ‘Harmonic frequencies of H on the rigid W(100) surface.
1 exp[
Jldw’SldwQ 0
exp
+A
1
(2.22)
with A given by eq 2.21. The only difference of eq 2.22 from the earlier expression for the case of only the two degrees of freedom s and Q, eq 2.12, is the extra term A in the exponent of the integrand. Thus it is a little more difficult to carry out the Monte Carlo path integral here, eq 2.22, than the two-dimensional case. To the extend that the coupling between the (qk)and (s,Q) degrees of freedom is h e a r in qk, this provides a powerful method for including truly many degrees of freedom. Finally, we note that the above expressions are essentially a quantum version of the classical generalized Langevin approach. The extra action A in the exponent of the integrand of eq 2.22 is the quantum counterpart to the “memory” terms and “random force” in the classical generalized Langevin e q ~ a t i o n . ~It , may ~ be though of as a quantum mechanical “friction” that the s and Q degrees of freedom experience from their coupling to the oscillator bath modes. 111. Application to H on W(100) A . The Model. The LEPS potential function devised by McCreery and Wolkenloa was used to descrbe the interaction of H on W(100). This model is based on a valence bond picture and was originally developed to describe the interaction of two hydrogen atoms with a surface; it should be straightforward, therefore, to extend the present studies to deal with the recombination/desorption process, eq 1.1. The solid surface is assumed to be rigid, providing a static background potential in which the gas molecules move. The Morse parameters for H-W( 100) adsorption were estimated from extended Hiickel calculations of Anders et al.,I9 and in order to account for variation in H-atom binding energy in the plane of the surface (the x-y plane), the Morse parameters (well depth, curvature, and equilibrium position) are assumed to be functions of x and y . The specific form of the potential function is given in the Appendix. Figure 2 of ref 10a shows the geometry of the tungsten surface and the various symmetric sites for an H atom on it; these are designated l C N , 2CN, and 5CN according to the coordination number of the hydrogen atom at the site. Table I gives the geometry and energetics of an H atomm at these three sites and also the three harmonic frequencies from a vibrational analysis of the LEPS potential function. The on-top site (1CN) is the stable location, the bridge site (2CN) is the transition state between two such stable sites, and site 5CN is a ”hypertransition state”, Le., it has two imaginary frequencies. Figure 1 shows a contour plot of the H rigid surface potential in the variables x,z (with y = 0)8 and also shows the minimum energy path through the transition state between two stable minima. As noted in the Introduction, this potential for H-atom adsx.ption is in contradiction to recent experimental work” which indicates that the stable minimum is actually the bridge (2CN) site. The interpretation of the experiments is clouded, though, because the tungsten surface reconstructs several times with (19) L. W. Anders, R. S . Hansen, and L. S. Bartell, J . Chem. Phys., 59, 5277 (1957); R. Gomer, R. Wortman, and R. Lundy, J . Chem. Phys., 26, 1147 (1957); H. Froitzheim, H. Ibach, and S . Lehwald, Phys. Rev. Lett., 36, 1549 (1976); P. Thiry in ref l f , p 231.
Z (a.u.1 Figure 1. Equipotential contours for a H atom on the (rigid) W(100) surface, from the potential function of ref loa. The dashed curve indicates the reaction path from the stable site 1CN through the transition state (site 2CN).
varying H-atom coverage. This is one reason that no reliable data yet exist for hydrogen atom diffusion on W(100). In the present paper, therefore, the McCreery-Wolken potential is viewed as a model system to illustrate qualitative effects. Better systems for future work, for which the experimental situation is clearer, or Ni( are hydrogen on Cu( Following McCreery and Wolken,Iof the LEPS potential function is generalized to include surface atom motion as follows:
V(7d = VLEPS(?) + VR(q) + V C O R R ( 7 d
(3.1)
where r‘ = (x,y,z) denotes the three coordinates of the H atom, and q = {qk)are the normal-mode coordinates for displacement of the surface atoms from their equilibrium positions. VR is the harmonic potential of the surface atoms vR(q)
=
C’/2mwk2qk2 k
(3.2)
and Vco, describes the change in the gassurface potential when the surface atoms are distorted from their equilibrium positions. ( VcoM(F,q) = 0 for (qk]= 0.) Following the spirit of section IIc, as linear in the surface atom displacements we approximate V, (qkl VCORR(r‘d -uk(?)qk (3.3a)
-
k
where
The potential function Vof eq 3.1 then has the form as that in section IIc (except for the fact that there are two Q-like degrees of freedom rather than the one denoted there). (20) A. Gelb and M. F. Cardillo, Surf.Sci., 59, 128 (1976); A. Gelb and M. F. Cardillo, Surf.Sci., 64, 197 (1977); A. Gelb and M. F. Cardillo, Surfi Sci., 75, 199 (1978). (21) T. N. Upton and W. A. Goddard, Phys. Reu. Lett., 42,472 (1979); W. Ho. N. J. DiNardo, and E. W. Plummer. J . Vac. Sci. Technol., 17. 134 (1980); P. Thiry in ref If, p 231; M. J. Push, R. M. Nieminen, M. Manninen, B. Chakraborty, S . Holloway, and J. K. Nsrskov, Phys. Rev. Lett., 51, 1081 (1983); B. I . Lundquist, B. Hellsing, S. Holmstrom, P. Nordlander, M. Persson, and J. K. Nstskov, Int. J . Quant. Chem., 23, 1083 (1983).
Diffusion of H on W ( 100)
The Journal of Physical Chemistry, Vol. 89, No. 11, 1985 2143
TABLE II: Tunneling Coefficients for Diffusion of H on the Rigid W(100)Surface T, K kTST,Os-l omT? cm2/s 2 300 1.96 X lo8 6.99 x 10-7 1.16 200 140 100
3.78 X lo4 4.60 X 6.27 X
1.06 X lo6 1.29 x 103 0.176
1.27 1.35 4.3
‘The transition-state theory rate, eq 3.9. *Equation 3.7 with kTST. The interaction potential VCoRR is specified initially in tezms of the atomic positions {RiJof the surface atoms, VCORR(F,{Ri)), so that the normal-mode problem for surface atom qotion must be solved. This provides th_e transformation matrix L i , k relating surface atom coordinates {R$and the normal-mode coordinates, {qkl
mw‘I2(& - &,q) = &Ck
k
(3.4)
so that the function f k ( 3 , eq 3.3b, is then given by I
1 I
2
3
wK/566 crn-l
As in ref lOf, we have taken V , , to be a sum of pair interactions between the H atom and the surface atoms
where Wis a Morse potential with parameters CY = 0.5361 and D = 0.0503 (au). In the present, preliminary study we have included only one surface phonon mode, the relative stretch along the x axis of the two tungsten atoms above which the H atom is tunneling. This is the surface motion that clearly interacts most strongly with tunneling of the H atom,and is therefore of most interest to investigate initially. The normal-mode coordinate qk for this mode is
where {Xi)is the x coordinate of surface atom i. The frequency of this photon mode was chosen as w k = 566 cm-I, which corresponds to the maximum in the density of phonon states in the classical model of Shugard et al.4a B. Results. With the usual assumption that the diffusion is via uncorrelated “hops” from one stable minimum on the surface (Le., site 1CN) to another-which is reasonable in the tunneling regime where the residence time at a given stable site is many vibrational periods-the diffusion coefficient D is given by22 a2
D=-k 2d
(3.7)
where a = 5 . 9 7 ~is ~the lattice constant (the distance between stable sites), d = 2 is the dimensionality, and k (s-l) is the rate of “hopping” from one stable site to another. It is the rate k that is calculated by the method described in section 11, and it is useful to express it as k = rkTsT
(3.8)
where kTST is the rate given by transition-state theory (3.9) and I’is the correction factor for tunneling, nonseparability of the various degrees of freedom, etc. Z,, Z,, and Z , are the (harmonic) vibrational partition functions of the reactant, i.e., (22) G. A. Sornorjai, ‘Principles of Surface Chemistry”, Prentice Hall, Englewood Cliffs, NJ., 1972; H. S. Johnston, ‘Gas Phase Reaction Rate Theory“, Ronald, New York, 1966, p 44; W. H.Miller, “Potential Energy Surfaces and Dynamics Calculations”, D. G. Truhlar, Ed., Plenum, New York, 1981, p 265.
Figure 2. r, the correction factor to the transition-state theory rate, for H atom diffusion on a W(100) surface, at T = 300 K. The horizontal
dashed line is for the case of a rigid surface, and the solid curved lines are the results for a moveable surface. wk is the frequency of the surface phonon, and the various curves correspond to the strength of coupling to the phonon; X1 is the strength function as calculated from the potential function of section IIIa, and X2, X3, and X5 correspond to increasing it by a factor of 2, 3, and 5 , respectively. the stable site l C N , and Z,,* and Z,* the (harmonic) partition functions for the two real frequency vibrations at the transition state (site 2CN). The partition functions for the phonons are the same in numerator and denominator and thus cancel each other. Table I1 gives the rate over the temperature range 100-300 K for the rigid surface case. Because the imaginary frequency at the transition state is small, 3721’ cm-’, the tunneling factor r is not much greater than unity except at the lowest temperature. For the above calculations via the methods in section 11, a value of N (the number of divisions of the hb ”time” interval) N 15-20 with -5000-10000 random walks was sufficient to give convergence. The effect of surface atom motion was studied by use of the one-phonon model described in section IIIa with the method described in section IIc. The resulting rate is also expressed in the form of eq 3.8, so the interesting result is to see how r depends on the phonon coupling. Because of the modelistic nature of the present study, the calculations were carried out not only for the phonon frequency tdk = 566 cm-I but also for larger and smaller values to illustrate how r varies with it. Similarly, calculations were carried out with the coupling functionfk(3 as calculated from eq 3.3-3.6 and then also for the case that this was increased by a factor of 2, 3, and, in some cases, 5 to show how r varies with the strength of coupling to phonons. Figure 2 shows how r varies with phonon frequency and with coupling strength for T = 300 K , and Figure 3 gives similar results for T = 140 K. In all cases, coupling to the phonon increases r, and thus the diffusion rate, and the increase is greater the larger the coupling strength. r varies inversely with phonon frequency. For the “standard” values used for the phonon frequency and coupling strength, the effect of the phonon interaction is to increase r by 10% at 300 K and by -20% at 1 4 0 K. As seen, though, the increase in r becomes greater rapidly if a lower phonon frequency is used.
-
IV. Concluding Remarks The results of this admittedly preliminary study are quite interesting and encouraging. First, it is encouraging that one is, in fact, able to carry out rigorous quantum mechanical calculations for the dynamics of an adsorbed species on a nonrigid surface. The interesting qualitative result is that coupling to the phonon
Jaquet and Miller
2144 The Journal of Physical Chemistry, Vol. 89, No. 1I, 1985
Chemical Sciences Division of the US.Department of Energy under Contract EA-AC03-76SF00098. The calculations were carried out on a Gould 8780 minicomputer funded by the National Science Foundation Grants CHE-79-2018 1 and CHE-83-20487. R.J. thanks the Alexander von Humboldt Foundation for a Feodor Lynen grant and the whole of the Miller group for many helpful discussions.
Appendix: LEPS Potential for Atom/Diatom and Solid Surface The diatom-solid surface potentialIoa has the form VLEPS=
U, + U2
+ U3 - [AI2+ (A2 + A# - Al(A2 + A3)]’/2
where Vi describes a Morse potential
u. = I
Di 4(1
+ Ai)((3 + Ai) exp[-2ai(ri (2
- ria)] -
+ 6AJ
exp[-ai(ri - ria)])
and Ai an anti-Morse potential Ai = I
I
I
2 ~ ~ 1 5 cm-‘ 6 6
Di 4(1
+ Ai) ((1 + 3Ai) exp[-2ai(ri - ria)] (6
3
Figure 3. Same as Figure 2, but for T = 140 K.
+ 2Ai) exp[-cri(ri - ria)])
where by introducing a Sato parameter A additional flexibility is obtained. Di, ai, and rio are the dissociation energy, Morse parameter, and equilibrium distance for the ith two-body interaction. Symmetry in the x-y plane is explicitly built in by defining
degrees of freedom tends to increase the tunneling (and thus diffusion) rate. D(X,Y) = DOL1 + &!(X,Y)l Although this major qualitative result, Le., that the coupling increases the tunneling rate, seems at variance with other ro(x,y) = z,[l + W X , Y ) l which concludes that frictional effects quench tunneling, there where is actually no inconsistency. This is because in the present microcrospic model the “coupling” induces a shift in the effective Q(x,y) = k[ cos cos potential the tunneling coordinate sees as well as dissapative effects on the tunneling degrees of freedom. In these other ~ o r k s ~ ’ . ~ ~ a phenomenological model was used, where only the dissipative A[cos - 1][ cos -11 effect of the coupling was considered and not a shift in the effective potential. There are many avenues that are suggested for future studies. First, it would be useful to pursue the present model by allowing P ( x , y ) = k[ cos + cos for many phonons. To do this in a rigorous fashin would require carrying out a normal-mode analysis for a cluster of surface atoms in order to generate the _distribution of frequencies (wk)and the B [ cos - 1 1 cos -11 transformation matrix Li,k needed to determine the coupling functions lfk(7));this is certainly feasible. The path integral The following parameters are calculated by fitting to experimental calculation via section IIc should not be much more difficult to and theoretical results to reproduce the data in Table I. carry out for the case of many phonons than it is for one. It would also be of interest to parameterize the McCreery-Wolken potential ~ H - w ( x , ~=) [0.02894 a u / D ( ~ , y ) ] ” ~ function to describe diffusion of H atoms on other, better unDo = 2.74 eV, 6 = 0.0474, z , = 2.19 au derstood surfaces. Of even more interest, though, would be to apply these methods A = 0.08, c = 0.21121, A = 1.463, B = 0.653 to the dynamics of two H atoms on a surface, in order to be able to describe the simplest chemical reaction, eq 1.1, and such work For the H-W interaction we use a Morse potential is indeed contemplated. vH-W = D(x~y)[e~P(-~aH-W(z - ZO)) - exp(-aH-W(z - ZO))] Acknowledgment. This work has been supported by the Diwith parameters rector, Office of Energy Research, Office of Basic Energy Sciences, Do = 4.7466 eV, Ro = 1.40083 au, CYH-H = 1.04435 au
($)
+
(r)]
(T)
(%)
($)
(T)]
(y)
(23) A. 0.Caldeira and A. J. Legett, Phys. Reu. Lett., 46, 211 (1981).
Registry No. H, 12385-13-6; W, 7440-33-7.
[ (7 )
