- Email: [email protected]
The energetics and dynamics of H2 dissociation on Al(110)
Z...
(:,:,, .....‘.~.~,,.;;,.,li,,,
:“:“:::::::~:.:‘:::::::::::::~:::::::~:~:::~:~:~~~~~~~:~~~~~~
“.‘....V.‘i . . . . . ...._..............~....,,.,,~,,,~~,,,,,,, ~3iQii~,~~~:~:i:::::::~:.:.:.:.~.: .i................:
:.:.:.:.:.:i.
.i. ): :.::::::
<..
science .surface ... .....:::i:;:
::::::::::::::::~.:.:$ i:::::~.:,.:.:.:........ ................,.i,.(......_., ,,..., ,,,(,, ..__. :‘,‘;:“.“.. .............._,,,, ...._.........~. “““‘.“A.. ....i.. ..,.,.... ...........:.x.: ....... .........,,., ..~,_.,,, _,:, ,, ~:.,. .:~.‘.““.~.‘.~~~~ ‘i”’ ‘“‘:.:.:x .i......... ..: .,.....,.,: .. ,,,,.,‘.,.,:,!,,,
ELSEVIER
Surface Science 304 (1994) 131-144
The energetics and dynamics of H, dissociation on Al( 110) Kent Gundersen,
Karsten W. Jacobsen, Jens K. Norskov *
Physics Department, Technical University of Denmark, DK-2800 Lyngby, Denmark
Bj@rk Hammer Fritz-Haber-Institut der Max-Planck-Gesellschaft, D-14195 Berlin (Dahlem), Germany
(Received 3 September 1993; accepted for publication 16 November 1993)
Abstract We present a modeling of the sticking dynamics of H, on Al(110). The modeling is based on an ab initio calculation of the H,/Al(llO) potential energy surface. The calculation is done both within the local density approximation (LDA) and using non-local corrections via the generalized gradient correction (GGA). We find that the GGA increases the barrier for dissociation substantially, and that the inclusion of gradient corrections greatly influences the dissociation dynamics. The dissociation dynamics is simulated by first doing a full quantum mechanical calculation of the dynamics in two of the six H, coordinates, the distance of the molecule from the surface and the intra-molecular bond length. The other four degrees of freedom are then included using the hole model. This implies that these coordinates are treated classically, and in the sudden approximation. Test calculations in three dimensions, where full quantum calculations are feasible, show the hole model to work well for the onset of
sticking. Finally, the full six-dimensional sticking calculations experiments, both for pure and seeded beams.
1. Introduction In spite of the importance of the sticking and dissociation process in surface science and more generally in an understanding of reactions at surfaces it is only with the advent of molecular beam studies of the dissociation process that enough experimental data has become available for a detailed testing of the theoretical description of these processes. Hydrogen dissociation has turned out to be the standard test case. Beam experi-
* Corresponding
author.
are compared
to the results of molecular
beam
ments have been performed for a number of systems [l-.51 and the H, dissociation process has been by far the most studied theoretically. Calculations of the interaction potential have been performed based on the jellium model [6,7], larger and larger clusters [g-lo], approximate energy methods [11,12], and most recently slab [13] and Greens function methods [141. Simulations of the dynamics have generally been based on model potentials and usually only a few degrees of freedom (typically the Hz-surface distance, Z, and the H-H bond length, b) have been included [15-171. This approach has developed a very important conceptual understanding of the depen-
0039-6028/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved XSDZ 0039-6028(93)E0876-V
132
K. Gundersen et al. /Surface
dence of the dissociation and the coupling between dissociation and vibrations of the incoming molecule on the topology of the potential in the Z-b plane. Only the approximate potentials have been able to provide knowledge about the full six-dimensional potential of the H, molecule outside a metal surface, and these calculations show that the potential is strongly corrugated in the other four degrees of freedom indicating that a full understanding of the dynamics and direct comparisons to experiment are only meaningful if all degrees of freedom are treated [12]. In the present paper we present a detailed calculation of the potential energy surface for H, outside an Al(110) surface using a Car-Parrinello-like total energy method with planewave pseudopotentials in a slab geometry [18]. Preliminary accounts of these calculations have been presented elsewhere [13,19]. Here we complement these by calculations for a system where the HZ-H, distance is larger, and which is therefore expected to be closer to the zero-coverage limit. We then study the dynamics on this potential including all six degrees of freedom. We do a full quantum treatment in two and three dimensions and test in this way the so-called “hole model” [20] proposed to include approximately the four degrees of freedom outside the Z-b plane. Having established that the hole model works well at least for the onset of sticking, we proceed to model the sticking including all six degrees of freedom of the H, molecule and compare this to experiment.
2. The potential energy surface In the present work the dynamics of the H, dissociation has been studied using an ab initio interaction potential determined from total energy calculations using planewave pseudopotentials. First, the reaction path geometry was determined in Car-Parrinello-like calculations involving both electronic and ionic degrees of freedom, second, potential energy surfaces were mapped and vibrational frequencies were determined by a number of calculations for fixed ionic coordinates. The most important approximation in the
Science 304 (1994) 131-144
total energy calculations is in the modeling of the electronic exchange and correlation energy. Widely used for this purpose is the local density approximation (LDA). However, it has been found by the present authors, that the inclusion of non-local effects in the exchange-correlation description drastically changes the energy barrier compared to the LDA result. This was found using the generalized gradient approximation (GGA) but is to be expected for all non-local descriptions improving on the LDA [19]. The GGA has been shown in a number of studies to consistently improve binding energies of atoms and molecules [21,22]. For solid systems the LDA is known generally to overbind and the dominating effect of the GGA is to weaken the binding leading to improved cohesive energies [23-251. However, in some cases, like Pd where the LDA overbinding is absent, the GGA overcorrects with too small cohesive energies and bulk moduli as a result [26]. In the present work we choose to discuss the dynamical interaction of H, with AK1101 on potential energy surfaces derived within both of the two exchange-correlation approximations. This allows us to comment on the sensitivity of the final estimates of the sticking coefficient to the exchange-correlation approximation. Besides the exchange-correlation approximation we consider the geometrical constraints periodic boundary conditions - to cause the only other non-vanishing error in the potential energy surfaces. The influence from this factor is therefore also monitored by calculating potential energy surfaces for two different boundary conditions. In Fig. 1 we show the geometries of the H,/Al(llO) system used for the mapping of the potential energy surfaces (PESs). Prior work has shown [13,27] that the minimum energy reaction path (defined as a steepest descent from the smallest possible dissociation barrier) with either of the exchange-correlation approximations is the one where the molecule lies parallel to the surface, centered over and aligned along the so-called long bridge of the (110) surface. In determining this path no assumptions of certain symmetries were made. The relaxation of the surface was neglected, as the surface atoms are much heavier
K. Gundersen et al./Surface Science304 (1994) 131-144
than the H atoms and hence do not have time to move appreciably during the dissociation event. The minimum energy reaction path only involves two ionic degrees of freedom, the distance, 2, of the molecular center of mass from the outermost (IlO)-layer of Al nuclei and the H-H bond-length, b. The reaction path is illustrated in Fig. la. The final state of the dissociation is that of the H atoms chemisorbed on-top of the outermost Al atoms. This position implies a high frequency of the vibration of H perpendicular to the surface. We calculate the frequency to be 1775 cm-’ (220 meV) in good agreement with the experimental value of 1810 cm-’ [28]. The potential energy surfaces are mapped by performing a number of total energy calculations for fixed ionic coordinates. Ten different b-values combined with thirteen different Z-values where chosen for the mapping. Contour plots based on these total energies are shown in Fig. 2. Figs. 2a and 2b respectively show the LDA and GGA
2.25 2.00 1.75 03
1.50
N
1.25 1.oo 0.75 2.25 2.00 1.75
03
1.50
N
1.25 1.oo 0.75
1.0
1.5
2.0
2.5
3.0
Fig. 2. Contour plot of the Z-b plane potential energy surfaces. (a) LDA potential for the supercell of Figs. la and lb, (b) GGA potential for the supercell of Figs. la and lb, (c) GGA potential for the supercell of Fig lc. The thick lines indicate the reaction paths. The contour lines run from 0.1 to 1.0 eV with intervals of 0.1 eV.
(a>
09
133
@>
Fig. 1. (a) Perspective view of the H, /AKIlO) system. Small spheres represent H atoms, large spheres Al atoms. The dark line is the minimum energy reaction path. (b) The (1 x2) surface supercell of an orthorhombic supercell. (c) The surface supercell of a monoclinic supercell.
based potential energy surface for the supercell of Figs. la and lb. Fig. 2c displays the GGA based potential for the supercell of Fig. lc. It is seen that all three potential energy surfaces of Fig. 2 exhibit very similar topologies, it is only the barrier heights that differ. For the LDA calculation we get a barrier height of 0.25 eV, for the GGA it is 0.54 eV and 0.70 eV for the two
K. Gundersen et al. /Surfnce
134
different supercells respectively. The physics behind the barrier increase from the LDA description to the GGA description has been addressed in Ref. [13]. We find that the increase in barrier is related to an increased energy cost due to orthogonalization of the surface metal electron states to the H, states which are more tightly bound than in the LDA when non-local exchange-correlation effects are included. The difference between the barriers as calculated within the GGA using the two different supercells of Figs. lb and lc is 0.16 eV, suggesting a substantial HZ-H, interaction, at least in the cell of Fig. lb. HFre the HZ-H, separation is 2.8 A, while it is 4.9 A in the cell of Fig. lc. While we cannot be certain that the inter-molecular interactions are negligible even in the last cell, we suggest this result is close to the zero-coverage limit. Further test of this would require larger unit cells, which is not possible computationally at the moment, or a method like the one suggested by Feibelman [14]. The vibrational energies at the transition state have been determined through a normal mode analysis. The total energy was fitted to a quadratic function around the saddle point. The fit was done over energies typical to the ground-state motion. In case the directions of the modes were not given by symmetry a proper diagonalization of the dynamical matrix was performed. In Table 1 the result is given. The frequencies are seen to
Table 1 The vibrational
frequencies
of the normal
modes
at the transition
Science 304 (1994) 131-144
be quite independent of whether the LDA or the GGA based potentials are used. We find only small differences between the two geometries of Figs. lb and lc with respect to the frequencies. All the frequencies of the normal modes are sizable. This means that the barrier for dissociation depends strongly on the other four degrees of freedom and therefore that they should be included in the investigation of the dynamical behavior of H, on Al(110). The coming sections of this paper will focus on how to include all six dimensions without having to map the potential energy surface entirely and perform a full 6D dynamic simulation. Details of the total energy calculations now follow - the reader might prefer to skip directly to the next section for the description of the model of sticking. Aluminum slabs of five (110) layer thickness separated by five layers of vacuum are used. The theoretical lattice constant, 3.96 A, is used. H, molecules approach both sides of the slabs and the final coverage after dissociation is one H per surface Al for the supercells of Fig. lb as well as for Fig. lc. As stated, the total energy calculations are performed within either the LDA (Ceperley-Alder [29]) or the GGA (The GGA-II version from 1991 by Perdew and Wang - see Ref. [21] and references [26,27] therein). Aluminium is represented by a pseudopotential of the Kerker type [30] in the Kleinman-Bylander form [31]. This ab
states of the minimum
energy
paths
calculated
within
the LDA
and the GGA
Reaction
path
Mode in the 6-Z
plane
Frustrated translation in the Y-8 plane Frustrated rotation in the Y-8 plane Rotation parallel to the surface, 4 Translation
along the close-packed
a) Calculated with b, Calculated with All numbers are coordinate for the
rows, X
the supercell of Fig. lb. the supercell of Fig. lc. in meV. Y is the coordinate molecule.
LDA
GGA a’
GGA b,
ho*
1lOi
1lOi
1lOi
gw, ho, gw, h%
120 80 180 70
110 60 170 80
110 _
h%
40
70
perpendicular
to the close-packed
rows.
_ _
19 is the out-of-surface-plane-rotation
K. Gundersen et al. /Surface
initio non-local pseudopotential was generated within the LDA. As the adsorption properties of the surface are expected not to depend on the aluminium core levels which are most affected by the non-local exchange-correlation effects, we use this potential for the GGA study too. This procedure is justified by work by Garcia et al. [26]. The hydrogen atoms are represented by the bare Coulomb potential and all non-local effects are therefore included here. The bare Coulomb potential naturally sets a rather high demand for the basis set required to represent the electronic states. We choose to include planewaves with a kinetic energy up to a cutoff of 500 eV and find that this causes an absolute error less than 0.02 eV/molecule on the barrier heights. As the GGA requires quantities like 1V,,I/n one might be concerned that a rather dense FFT grid would be required to obtain convergence of the exchange-co~elation energy. However, direct tests have shown that the FFT grid which includes all G - G’ vectors of basis set planewave vectors, G and G’, suffices to converge the exchange-correlation energy to less than 1 meV per molecule. We choose k-point grids according to Monkhorst and Pack [32], with (4, 8) and (6, 6) k-point grids for the surface Brillouin zones of the supercells of Figs. lb and lc respectively. The error in the total energies caused by this choice of k-point grids is negligible compared to the error due to the finite basis set. The solution of the Kohn-Sham equations is found by means of an all-bands minimization technique [27] representing a slight generalization of the band-by-band minimization technique of Teter et al. [33]. Occupation numbers for these solutions are found by the Gaussian smearing method of Fu and Ho [34] as modified by Needs et al. [35]. The smearing width is 200 meV. All calculations both within the LDA and the GGA are iterated to self-consistency. In order to improve the stabili~ of the calculations involving the GGA, all gradient dependent terms in the exchange-correlation evaluation are neglected for densities smaller than 0.00035 A-“. This is done purely for numerical reasons and represents no further approximation.
Science 304 (1994) 131-144
135
3. The model of sticking The dynamical problem of H, dissociating on metal surfaces concerns six ionic degrees of freedom (neglecting all the degrees of freedom of the surface ions). In the particular case of H, dissociation on Al(llO), where there is a sizable corrugation of the PES in all of the six dimensions [13], all of these have to be included in a modefing of the sticking coefficient [12]. Due to the small mass of the H, the dynamical problem preferentially should be studied quantum mechanically. Up to now this has not been possible due to ~mputational constraints [36]. Therefore some approximations to the six-dimensional quantum dynamical treatment must be made. We describe below an extension of the hole mode1 of Karikorpi et al. [20], that includes a full quantum dynamical description in two dimensions while only estimating the dynamical effects of the remaining four dimensions. The two degrees of freedom that are treated in a full quantum calculation are the HZ-surface distance, 2, and the H-H bond length, b. As noted above, these degrees of freedom contain the entire minimum energy reaction path. Also contained is the vibrational dimension of the gas phase molecule. This allows for a proper representation of vibrationally ground state and excited state molecules, which is needed in the discussion of the coupling of the sticking to the vibrational state. These “first” two degrees of freedom are highly coupled as can be judged from the contour plots of Fig. 2. However, a change of coordinates into one along the reaction path and one perpendicular to it, will allow for a somewhat simpler representation of the PES of these two dimensions. Such a representation is given in Fig. 3, where the thin line is the potential energy along the reaction path, R. The vibrational energy of the dissociating molecule can approximately be described as the energy levels in a harmonic fit to the potential energy in the coordinate perpendicular to the reaction path. Such vibrational energy levels along the reaction path have been determined for a molecule in the ground state, v = 0, and in the first excited state, Y = 1. In Fig. 3 these are added to the potential
136
K. ~~der~e~
et al. /Surface
0.8
r 2 0.6 6 g 0.4
w 0.2 0.0
0.8 F 5 0.6 $ Gi 0.4 5 0.2
O.OL’. : -2.0 -1.0
,/ ’
1
’
’
0.0 1.0 R [Al
’ 2.0
’
3.0
Fig. 3. The potential energy per H, molecule (thin line) along the minimum energy path, the effective ID PES of a molecule in the vibrational ground state (thick line) and of a molecule in the first vibrational state (dotted line). The vibrational energy is determined at the levels, $Aw(R) and $ho(R), in a harmonic fit to the potential energy perpendicular to the reaction path, R. The latter is defined as the sum of (AR1 in A, where R is the six-dimensional vector of the Cartesian Soordinates of the two hydrogen atoms. R is zero for Z-2.3 A. wa is the intra-molecular vibrational frequency of the gas phase H, molecule. Each subfigure refers to a different PES: (a) the LDA based PES of Fig. Za, (b) the GGA based PES of Fig. Zb, Cc>the GGA based PES of Fig. 2~.
energy along the reaction path and provide adiabatic, quantum mechanical estimates of the effective 1D potential energy surfaces (thick and dotted lines) [16]. The effective 1D PESs illustrate
Science 304 (19941 131-144
the “bootstrapping” effect [15]; as the vibrational frequency of the molecule is softened as a result of the molecule-surface interaction, a molecule will encounter an effectively smaller barrier for dissociation. The remaining four degrees of freedom, (the two molecular rotations and the two translations parallel to the surface), we now have to model approximately. The first approximation is to consider these degrees of freedom individually, which implies that they only couple weakly to the first two, and that they couple only weakly with one another. The second approximation is to neglect any velocities in these dimensions. This means that we model only molecules impinging normal to the surface having no rotational energy (J = m = 0). Next we have to choose how to model the dynamical behavior in these dimensions. Candidates for simple models are the extreme cases, i.e. either the adiabatic or sudden limit. And for each of these in either a quantum mechanical or a classical picture. As the molecule approaches the surface the potential in the remaining four degrees of freedom increases from zero frequency to some finite value. If, on the timescale of the vibrational period, the frequency increase is small compared with the frequency itself, then we are in the adiabatic limit. In a quantum mechanical picture, this adiabatic limit implies that the molecule has time to “build up” the zero-point energy, ihw;, in each of the four degrees of freedom. The energy has to come from the kinetic or the initial intra-molecular vibrational energy of the scattering molecule, as illustrated in Fig. 3 in the Z-b plane. A similar energy “build up” takes place in a classical picture. Here the vibrational energy, E, divided by the frequency, o, for each of the additional degrees of freedom is an adiabatic i~~~~ia~t 1371. Thus when the frequency increases the energy “bound” in the corresponding mode increases proportionally to the frequency - like in the quantum mechanical picture. The other extreme is the sudden limit. If only a fraction of a vibrational period can be completed, while the frequency change takes place, the molecule does not have time to react on the change in the potential. The molecule will there-
K Gundersen et al. /Surface
fore interact with unchanged coordinates in the third (and higher) degrees of freedom. Classically, this implies that the molecule encounters a surface with a distribution of barrier heights. The effect of the extra degrees of freedom must therefore be included by integrating out the coordinates while weighting the sticking probability for each set of coordinates with the likelihood for the molecule of having these coordinates. For the purpose of the present sticking model, we adapt a classical description within the sudden approximation. That this is adequate can be judged from the following argument: The typical kinetic energies of impinging molecules are about 0.4 eV (for those that stick) and the potential strength (vibrational frequency) in each of the four degrees of freedom rises towards the transition state to values in the range 40-180 meV as given in Table 1 of the previous section. This gives of the order a few tenths of a period of vibration for a classical trajectory in the time the molecule spends in the region of appreciable molecule-surface interaction before it reaches the transition state. Now for the details of the hole model. We start out by considering the sticking coefficient, S,,(E), in the two dimensions of the Hz-surface distance, Z, and the H-H bond length, b. The S,,(E) is chosen to be the sticking coefficient resulting from a 2D quantum dynamical simulation. This we do as a split-time-propagation of a wavepacket [38] in the 2D potential energy surfaces of molecule-surface interaction described in the previous section. The wavepacket is initiated at some large Zvalue. First, the quantum mechanical solution of one vibrational state, v, is found in the 1D intersection of the PES at this large Z-value. This is done as a propagation in imaginary time. If a vibrationally excited state is desired a proper out-projection of lower lying states is done during this initial propagation. Next, the solution is given a Gaussian shaped kinetic energy, E, of width 5-15% of E in the Z-direction. The following 2D propagation is then performed until the products have reached the asymptotic regions of space. In the exit channel there is an absorbing potential making sure that we are not disturbed by back
137
Science 304 (1994) 131-144
reflection due to the periodic boundary conditions that are imposed for simplicity. The sticking coefficient is given by the fraction of the wavepacket that moves down the exit channel. Tests have been carried out to ensure that the sticking coefficient is well converged with respect to the finite time step of the wavepacket propagation and with respect to the finite phase space grid-density on which the potential and the wavepacket are represented. If the additional four ionic degrees of freedom, 13,4, X and Y, where all steered such that the lowest possible barrier was encountered by all impinging molecules then the 2D sticking coefficient, S,,(E), would also be the 6D sticking coefficient. However as argued above, for the kinetic energies of interest we are in the sudden limit, where the instantaneous ~9,4, X and Y-coordinates of the gas (or molecular beam) phase remain essentially unchanged during the reaction event. Therefore to get from S,,(E) to the sticking coefficient, S,,(E), of all six ionic degrees of freedom we now subtract from the energy argument of S,,(E) whatever is needed to overcome the barrier increase, AE,(8, 4, X, Y), due to the values of 13, 4, X and Y. The S&E AE,(8, 4, X, Y)) is then averaged over all configuration space to yield S,,(E). The barrier increase, AE,(B, 4, X, Y), we estimate from a decomposition of the considered four ionic degrees of freedom into normal mode coordinates, ui(8, 4, X, Y>, at the transition state. Each of these normal modes, ui, is associated with a frequency, oi, that describes to second order the behavior of the potential. The wi were determined in the previous section. The resulting barrier increase reads
i=3 (1)
where mu is the proton mass (assuming u properly normalized). We notice that the expansion above will be best for energies not too far from the lowest activation energy. The present description is therefore best suited for describing the onset of sticking. Now only the average over
K. Gundersen et al. /Surface Science 304 (1994) 131-144
138
configuration space remains. This is done by calculating the volume of the “hole” in configuration space, on,(E), available for sticking at kinetic energy E and properly normalizing by the volume of the entire configuration space, 0,. The volume of the “hole” becomes a,(E)
= 2j$,,(E
the final 6D sticking curves and compare with experiment, it is very instructive to consider separately the output from the 2D wavepacket simulations. Also some support will be given to the hole model in Section 4.2, concentrating on a comparison between 3D wavepacket simulations and 3D hole models.
- A&,(0, $7 X, Y)) 4.1. Sticking curves for 20 PES
X
sin 0 d0 d+ dX dY,
(2) where the integration is over the 4D configuration space and sin 0 is the Jacobian with 8 being the polar angle. The factor of two in front of the integral is due to the inversion symmetry of the molecule which means that there are actually two minimum energy reaction paths in all configuration space whereas AE, in Eq. (1) is expanded as if only one minimum energy reaction path was present. The total volume of possible configurations is: %=j
sin 0 dtl d+ dX dY= 47rA,,
(3)
4D
where A, is the area of the (1 x 1) surface cell. The sticking probability in this 6D hole model is finally &n(E)
=%(E)/%*
(4)
Here the hole model has been described for all six dimensions. It can, however, easily be modified to yield a 3D sticking curve, where besides b and 2 only 1 degree of freedom is taken into account. If - as will be the case below - the X-coordinate is considered, the sticking coefficient is given by S,,(E)
= ;
jS,,(
E - $2muw;X”)
dX,
(5)
0
where R, is now the length of the (1 X 1) surface cell in the X-direction.
4. Results and discussion In this section we present and discuss the sticking curves that emerge from the sticking model of the last section, when applied to the PESs of Section 2. However, before we present
Much understanding of the H,/Al(llO) system can be gained already at the level of the first part of the model, the 2D quantum dynamical simulation of the reaction event. In Fig. 4 are shown the 2D sticking curves for the PES of Fig. 2. For each PES is shown the sticking curve for an impinging molecule, which is in its vibrationally ground state, v = 0, and in its vibrationally first excited state, v = 1. The sticking curves for the LDA based PES (lowest barrier) is shown in Fig. 4a. The sticking of vibrational ground state molecules is seen in general to increase with increasing energy. However a pronounced suppression of the sticking takes place in the kinetic energy range from 0.1 to 0.3 eV. This suppression can be understood in terms of the impinging molecules being trapped in the metastable v = 1 states in front of the barrier (see Fig. 3a). The trapping of the molecules in such metastable states strongly affects the sticking behavior and may cause the observed decreased sticking probability. The effect is discussed in detail by Halstead and Holloway [39]. That the effect is observed in the right energy range can be judged from the v = 1 curve of Fig. 3a. The bottom of this dip corresponds to a kinetic energy, 0.09 eV, of a ground state H, and the curvature of the dip along the reaction path coordinate corresponds to a vibrational frequency of 0.25 eV. A metastable level is therefore expected for kinetic energies around 0.22 eV. The results of the 2D dynamical simulations on the two GGA PESs are shown in Figs. 4b and 4c. Again the sticking probability increases with the kinetic energy for small kinetic energies. For larger kinetic energies the curves bend over and the sticking decreases with the kinetic energy. This is the same effect of interaction with
K. Gundersen et al. /Surface
1.0
i$
I
’
’
’
I
’
’
I
0.8
(d
% 0.6 ii 2 0.4 Ii .g 0.2 c/J 0.0
1.0
I
’
g 0.8 tu % 0.6 2 p 0.4 E .g 0.2 cn 0.0
139
H, in the Z-b plane is considered at different kinetic energies. The probability distributions are taken from the wavepacket simulations on the GGA based potential of Fig. 2b. At high energies the reflected beam shows excited H-H vibrations in the interaction region. Unlike the LDA based results, the GGA potentials with a high barrier for dissociation, give rise to differences between the vibrational ground state (v = 0) and the first excited state (v = 1) sticking at small energies; the v = 1 sticking curves has an 0.15 eV lower onset. We will return to this considerable coupling between the vibrational state of the wavepacket and the sticking coefficient when comparing with experiment below. 4.2. Verification of hole model - 30 study
1.0 $ 0.8 Cu “0 0.6 lc p 0.4 P .g 0.2 cn 0.0
Science 304 (1994) 131-144
0.0
0.2
0.4
0.6
0.8
Kinetic Energy [eV] Fig. 4. 2D sticking curves for the ground state (solid line) and the vibrationally first excited state (dotted line). Also shown is the 3D wavepacket simulated sticking curve for the ground state (dashed line). Each subfigure refers to a different PES: (a) the LDA based PES of Fig. 2a, (b) the GGA based PES of Fig. 2b, (c) the GGA based PES of Fig. 2c.
metastable states in front of the barrier as for the LDA PES, but now shifted to higher energies due to the more repulsive character of the two GGA PESs. The trapping in the metastable states implies that impinging molecules may become vibrationally excited due to the interaction with the surface. An indication that the interaction is indeed strong enough for this to happen is given in Fig. 5. Here the probability distribution for the
We have investigated the accuracy of the hole model by comparing 3D quantum calculations with the hole model in its 3D version, treating two degrees of freedom quantum dynamically and one degree of freedom in the classical sudden approximation through Eq. (5). The 3D potential has been determined within the LDA description. The 3D GGA potentials have been obtained by adding the difference between the 3D and 2D LDA potentials to the 2D GGA potentials. This means that the 3D potentials used are, to some extent, no more than model potentials. Since they are only used to ascertain the general applicability of the hole model this introduces no further approximations. The third dimension is chosen to be the translation along the close-packed rows, a very smooth mode which remains a normal coordinate all the way out from the surface. The vibrational frequency in this mode is 40 meV. In order to address the range of applicability of the hole model we also modify the mode to yield a vibrational frequency of 180 meV. Hereby the whole range of potential strengths of the four modes to be included is covered. The 3D quantum calculations are performed in a similar manner to the 2D calculations. The initial wavepacket is chosen to be the previously mentioned 2D wavepacket multiplied by a planewave in the third dimension.
140
(a) At barrier, E=0.2 eV
(b) Reflected wave, E=0.2 eV
(c) At barrier, E=0.4 eV
(d) Reflected wave, E=0.4 eV
(e) At barrier, E=O.6 eV
(f) Reflected wave, E=0.6 eV
Fig. 5. Probability density in the Z-b plane at different incoming energies. The PES is that of Fig. 2b. (a), (b) The scattering wavepacket for E = 0.2 eV. It is seen that the wavepacket is adiabaticaIly reflected from the potential with only very little tunneling. Cc), (d) The scattering wavepacket for E = 0.4 eV. There is a large barrier penetration and no vibrational excitation of the beam in the reaction zone. (e), (f) The scattering wavepacket for E = 0.6 eV. The wavepacket is clearly excited in the interaction region of the potential. Only a small barrier penetration is seen. Notice the excited tail of the outgoing wavepacket indicating that it has been trapped in an excited metastable state.
K. Gundersen et al. /Surface
In Fig. 6 are given the quantum dynamical 2D and 3D sticking curves. The sticking curves are seen to be very sensitive to the added third dimension. This confirms that all of the six molecular degrees of freedom must be included in a proper description of the sticking coefficient. Also included in Fig. 6 are the 3D hole model estimates of the sticking curve. It can be seen that the hole model captures the large decrease of the sticking probability due to the third degree of freedom, and that it is almost quantitative at the lowest energies (up to about 0.4 eV) where it is supposed to work best. The deviations are largest
141
Science 304 (1994) I31 -144
for the high frequency mode, the reason being that here the steering of the molecule towards the lowest energy barrier is strongest. As a consequence, the hole model systematically underestimates the sticking probability for this mode. The slight underestimate for the low frequency mode is, on the other hand, most likely due primarily to the second-order expansion of Eq. (1). 4.3. Results of the 60 hole model In Fig. 7 we show the 6D sticking curves obtained from our 2D quantum curves, using the
r
0.60
I
I
(b)
’
I
I
2D
1 0
0.50
i t .= h =
I
0.40
I-
% s z
I 0.30
i-
.-F
i
% z
0.20 0.005 0.10
Kin. Energy [eV] Fig. 6. The hole model strengths of the potential
applied to 3D sticking in the third dimension.
Kin. Energy [eV]
curves plotted together with the quantum (a) PES of Fig. 2b, (b) PES of Fig. 2c.
3D sticking
curve
for two different
K Gundersen et al. /Surface
Kinetic Energy [eV] Fig. 7. 6D sticking curves obtained from Eqs. (l)-(4) using the v = 0 ZD sticking curves of Fig. 4, (a) corresponds to the LDA PES while (b) and (c) correspond to the two GGA PESs (Figs. 2b and 2c respectively). The difference between (a) and (b) is caused by the inclusion of the non-local exchange-correlation effects, (b) and (cf differ due to the different Ha-H, separations with Cc) being the most converged result.
hole model, Eqs. (l)-(4). Comparing to the 2D sticking curves of Fig. 4 it is seen that the inclusion of all six degrees of freedom is essentia1 for a proper description of the sticking process. The large difference between the LDA and GGA results is striking. The GGA includes, in an approximate way, some non-local exchange-correlation effects not included in the LDA. One would therefore expect the GGA to be more accurate, but there is no reason to believe that it gives a
2
0.0008
r----
i-’
/Level
d
-7
of atomic H
-~E(T~YT~)I
in the beam ”; g
0.0006 0.0004
1 L
2 %
0.0002
L
g g
0.0000
i 0.0 Kinetic Energy [eV]
Fig. 8. The pure and seeded beam calculated from the 2D quantum curves of Fig. 4c using the hole model, Eq. (61, including the v-dependence of the sticking up to v = 2. The experimental result of Berger et al. [3] for the pure beam is also included. Also shown is the level of atomically adsorbed hydrogen found in the experiment. This has been subtracted from the raw experimental data to give the molecular sticking probability and thus gives an indication of the uncertainty.
Science 304 (19941 131-144
quantitative description. The difference between the LDA and GGA curves thus gives an indication of the accuracy by which “ab initio” sticking probabilities can be determined today. In order to compare the calculated sticking curves with the results of the beam experiments of Berger et al. [3] we must include the effect of vibrational excitations in the beam. In a beam experiment, both the translational energy and the degree of vibrationa excitation of the molecules in the beam, are given by the temperature of the nozzle. If the kinetic energy of the beam is changed by seeding techniques, the nozzle temperature and the beam energy can be considered independent variables, and the sticking probability as a function of normal kinetic energy E and nozzle temperature T,, is given by S( E, T,) = [ S,(E) + exp
X[I+exp[-$j
j-$hj”‘(“) +...I +...]-I,
(6)
where wO is the vibrational frequency for H, and S,(E) is the (6D) sticking probability - from Eq. (4) - for a molecule in vibrational state Y at the kinetic energy E. These have been summed with their Boltzman weights at the nozzle temperature to give the total sticking probabili~. For an unseeded beam, the relation between the translational energy and the nozzle temperature is E(T,)
= ;kBTn
(7)
and the sticking probability of the pure beam is therefore given by S(E(T,), T,). The model calculation, using what we believe to be our best results based on the Iargest HZ-H, separation possible and the GGA, is compared to the experimental result in Fig. 8. Considering the uncertainties in the tota energy calculation and in the description of the dynamics and the uncertainties in the experiments (see below) the agreement between theory and experiments is surprisingly good. In Fig. 8 we have also shown the sticking curve for a seeded beam with a fixed nozzIe tempera-
K. Gundersen et al. /Surface
ture of T,, = 1850 K calculated according to Eq. (6). Here we have included up to the second excited vibrational state. Higher excited states do not contribute, on the scale of interest, because their Boltzman factors are less than 10p4. It is clear that there is a difference between the pure and the seeded beam sticking curves but it is small on the scale of the experimental accuracy. The large contribution from atomically adsorbed hydrogen drowns the effect from the excited states and this is thus not seen in the experiment. It is therefore not surprising that Berger et al. suggest that there is no vibrational coupling for the H,/Al(llO) system [3]. In Ref. [13] we used a simplified version of the hole model to get a sticking probability from the LDA potential. This calculation differs from the present one in three respects. First, it was based on the LDA potential. Second, it had only a rough WKB estimate of the quantum dynamics in the Z-b plane, which among other things did not include the coupling between the H-H vibration and the sticking. And finally, it did not include the thermal averaging of Eq. (6). The present result is therefore substantially different from the earlier version.
5. Conclusion In the present paper we have presented a set of ab initio calculations of the H, potential outside an Al(110) surface. Further we have modeled the fully six-dimensional H, dynamics on this potential energy surface. We have shown that the treatment of exchange correlation effects greatly influences the dynamics. We have shown that a full quantum treatment of the dynamics in the Z-b plane combined with the hole model gives a good description of the dynamics, at least for the onset of sticking. Finally, we have compared the hole model calculation of the sticking probability directly to the beam experiments of Berger et al. [3] and shown that a reasonable agreement between theory and experiment can be obtained considering both the theoretical approximations and the experimental uncertainties.
Science 304 (1994) 131-144
143
6. Acknowledgements We are grateful to M.C. Payne for giving us the CASTEP (Cambridge Serial Total Energy Package) code. We are indebted to Stephen Holloway for pointing out an initial misinterpretation of the 2D sticking curves. We thank the Center for Applications of Parallel Computers (CAP) for providing computer resources. The research was done in part under the Center for Surface Reactions of the Danish Research Councils and in part under the Center for Atomic Scale Materials Physics (CAMP).
7. References 111 B.E. Hayden and C.L.A. Lamont, Phys. Rev. Lett. 63 (1989) 2823. 121 H.F. Berger, M. Leisch, A. Winkler and K.D. Rendulic, Chem. Phys. Lett. 175 (1990) 425. [3] H.F. Berger and K.D. Rendulic, Surf. Sci. 253 (1991) 325. [4] H.A. Michelsen and D.J. Auerbach, J. Chem. Phys. 94 (1991) 7502. 151 C.T. Rettner, D.J. Auerbach and H.A. Michelsen, Phys. Rev. Lett. 68 (1992) 1164. [6] H. Hjelmberg, J.K. Norskov and B.I. Lundqvist, Phys. Ser. 20 (1979) 192; P. Johansson, Surf. Sci. 104 (1981) 510. 171 J.K. Nsrskov, A. Houmoller, P. Johansson and B.I. Lundqvist, Phys. Rev. Lett. 46 (1981) 257. 181 J. Harris and S. Andersson, Phys. Rev. Lett. 55 (1985) 1583. [91 J.E. Miiller, Phys. Rev. Lett. 59 (19871 2943. [lOI J.E. Miiller, Surf. Sci. 272 (1992) 45. [ill J.K. Norskov, J. Chem. Phys. 90 (1989) 7461. 1121 C. Engdahl, B.I. Lundqvist, I-J. Nielsen and J.K. Norskov, Phys. Rev. B 45 (1992) 11362. [131 B. Hammer, K.W. Jacobsen and J.K. Norskov, Phys. Rev. Lett. 69 (1992) 1971. [14] P.J. Feibelman, Phys. Rev. Lett. 67 (1991) 461. 1151 J. Harris, Surf. Sci. 221 (1989) 335. I161 M.R. Hand and S. Holloway, J. Chem. Phys. 91 (1989) 7209. [171 S. Kuchenhoff, W. Brenig and Y. Chiba, Surf. Sci. 245 (1991) 389. [181 M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias and J.D. Joannopoulos, Rev. Mod. Phys. 64 (1992) 1045. I191 B. Hammer, K.W. Jacobsen and J.K. Norskov, Phys. Rev. Lett. 70 (1993) 3971. [201 M. Karikorpi, S. Holloway, N. Henriksen and J.K. Norskov, Surf. Sci. 179 (1987) L41. [211 J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson,
144
K. Gundersen et al. /Smface
M.R. Pederson, D.J. Singh and C. Fiolhais, Phys. Rev. B 46 (1992) 6671. [22] G. Ortiz and P. Ballone, Phys. Rev. B 43 (1991) 6376. [23] X.J. Kong, CT. Chart, K.M. Ho and Y.Y. Ye, Phys. Rev. B 42 (1990) 9357. [24] G. Ortiz, Phys. Rev. B 45 (1992) 11328. 1251 M. Kiirling and J. ~glund, Phys. Rev. B 45 (1992) 13293. [26] A. Garcia, C. Elsiisser, J. Zhu, S.G. Louie and M.L. Cohen, Phys. Rev. B 46 (1992) 9829. [271 B. Hammer, PhD Thesis, Physics Department, Technical University of Denmark (1993). [28] H. Kondoh, M. Hara, K. Domen and H. Nozoye, Surf. Sci. 287/288 (1993) 74. [29] D.M. Ceperly and B.J. Alder, Phys. Rev. Lett. 45 (1980) 566; J.P. Perdew and A. Zunger, Phys. Rev. B 23 (1981) 5048; J.P. Perdew and Y. Wang, Phys. Rev. B 45 (1992) 13244. [30] G.P. Kerker, J. Phys. C 13 (1980) L189.
Science 304 (1994) 131-144
[31] L. Kleinman and D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. 1321H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13 (1976) 5188. 1331M.P. Teter, M.C. Payne and D.C. Allan, Phys. Rev. B 40 (1989) 12255. [34] C.-L. Fu and K.-M. Ho, Phys. Rev. B 28 (1983) 5480. [35] R.J. Needs, R.M. Martin and O.H. Nielsen, Phys. Rev. B 33 (1986) 3778. [36] U. Nielsen, D. Halstead, S. Holloway and J.K. Nsrskov, J. Chem. Phys. 93 (1990) 2879. [37] L.D. Landau and E.M. Lifshitz, Quantum Mechanics (Pergamon, Oxford, UK 1977). [SS] J.A. Fleck Jr., J.R. Morris and M.D. Feit, Appl. Phys. 10 (1976) 129. [39] D. Halstead and S. Holloway, I. Chem. Phys. 93 (1990) 2859.
(:,:,, .....‘.~.~,,.;;,.,li,,,
:“:“:::::::~:.:‘:::::::::::::~:::::::~:~:::~:~:~~~~~~~:~~~~~~
“.‘....V.‘i . . . . . ...._..............~....,,.,,~,,,~~,,,,,,, ~3iQii~,~~~:~:i:::::::~:.:.:.:.~.: .i................:
:.:.:.:.:.:i.
.i. ): :.::::::
<..
science .surface ... .....:::i:;:
::::::::::::::::~.:.:$ i:::::~.:,.:.:.:........ ................,.i,.(......_., ,,..., ,,,(,, ..__. :‘,‘;:“.“.. .............._,,,, ...._.........~. “““‘.“A.. ....i.. ..,.,.... ...........:.x.: ....... .........,,., ..~,_.,,, _,:, ,, ~:.,. .:~.‘.““.~.‘.~~~~ ‘i”’ ‘“‘:.:.:x .i......... ..: .,.....,.,: .. ,,,,.,‘.,.,:,!,,,
ELSEVIER
Surface Science 304 (1994) 131-144
The energetics and dynamics of H, dissociation on Al( 110) Kent Gundersen,
Karsten W. Jacobsen, Jens K. Norskov *
Physics Department, Technical University of Denmark, DK-2800 Lyngby, Denmark
Bj@rk Hammer Fritz-Haber-Institut der Max-Planck-Gesellschaft, D-14195 Berlin (Dahlem), Germany
(Received 3 September 1993; accepted for publication 16 November 1993)
Abstract We present a modeling of the sticking dynamics of H, on Al(110). The modeling is based on an ab initio calculation of the H,/Al(llO) potential energy surface. The calculation is done both within the local density approximation (LDA) and using non-local corrections via the generalized gradient correction (GGA). We find that the GGA increases the barrier for dissociation substantially, and that the inclusion of gradient corrections greatly influences the dissociation dynamics. The dissociation dynamics is simulated by first doing a full quantum mechanical calculation of the dynamics in two of the six H, coordinates, the distance of the molecule from the surface and the intra-molecular bond length. The other four degrees of freedom are then included using the hole model. This implies that these coordinates are treated classically, and in the sudden approximation. Test calculations in three dimensions, where full quantum calculations are feasible, show the hole model to work well for the onset of
sticking. Finally, the full six-dimensional sticking calculations experiments, both for pure and seeded beams.
1. Introduction In spite of the importance of the sticking and dissociation process in surface science and more generally in an understanding of reactions at surfaces it is only with the advent of molecular beam studies of the dissociation process that enough experimental data has become available for a detailed testing of the theoretical description of these processes. Hydrogen dissociation has turned out to be the standard test case. Beam experi-
* Corresponding
author.
are compared
to the results of molecular
beam
ments have been performed for a number of systems [l-.51 and the H, dissociation process has been by far the most studied theoretically. Calculations of the interaction potential have been performed based on the jellium model [6,7], larger and larger clusters [g-lo], approximate energy methods [11,12], and most recently slab [13] and Greens function methods [141. Simulations of the dynamics have generally been based on model potentials and usually only a few degrees of freedom (typically the Hz-surface distance, Z, and the H-H bond length, b) have been included [15-171. This approach has developed a very important conceptual understanding of the depen-
0039-6028/94/$07.00 0 1994 Elsevier Science B.V. All rights reserved XSDZ 0039-6028(93)E0876-V
132
K. Gundersen et al. /Surface
dence of the dissociation and the coupling between dissociation and vibrations of the incoming molecule on the topology of the potential in the Z-b plane. Only the approximate potentials have been able to provide knowledge about the full six-dimensional potential of the H, molecule outside a metal surface, and these calculations show that the potential is strongly corrugated in the other four degrees of freedom indicating that a full understanding of the dynamics and direct comparisons to experiment are only meaningful if all degrees of freedom are treated [12]. In the present paper we present a detailed calculation of the potential energy surface for H, outside an Al(110) surface using a Car-Parrinello-like total energy method with planewave pseudopotentials in a slab geometry [18]. Preliminary accounts of these calculations have been presented elsewhere [13,19]. Here we complement these by calculations for a system where the HZ-H, distance is larger, and which is therefore expected to be closer to the zero-coverage limit. We then study the dynamics on this potential including all six degrees of freedom. We do a full quantum treatment in two and three dimensions and test in this way the so-called “hole model” [20] proposed to include approximately the four degrees of freedom outside the Z-b plane. Having established that the hole model works well at least for the onset of sticking, we proceed to model the sticking including all six degrees of freedom of the H, molecule and compare this to experiment.
2. The potential energy surface In the present work the dynamics of the H, dissociation has been studied using an ab initio interaction potential determined from total energy calculations using planewave pseudopotentials. First, the reaction path geometry was determined in Car-Parrinello-like calculations involving both electronic and ionic degrees of freedom, second, potential energy surfaces were mapped and vibrational frequencies were determined by a number of calculations for fixed ionic coordinates. The most important approximation in the
Science 304 (1994) 131-144
total energy calculations is in the modeling of the electronic exchange and correlation energy. Widely used for this purpose is the local density approximation (LDA). However, it has been found by the present authors, that the inclusion of non-local effects in the exchange-correlation description drastically changes the energy barrier compared to the LDA result. This was found using the generalized gradient approximation (GGA) but is to be expected for all non-local descriptions improving on the LDA [19]. The GGA has been shown in a number of studies to consistently improve binding energies of atoms and molecules [21,22]. For solid systems the LDA is known generally to overbind and the dominating effect of the GGA is to weaken the binding leading to improved cohesive energies [23-251. However, in some cases, like Pd where the LDA overbinding is absent, the GGA overcorrects with too small cohesive energies and bulk moduli as a result [26]. In the present work we choose to discuss the dynamical interaction of H, with AK1101 on potential energy surfaces derived within both of the two exchange-correlation approximations. This allows us to comment on the sensitivity of the final estimates of the sticking coefficient to the exchange-correlation approximation. Besides the exchange-correlation approximation we consider the geometrical constraints periodic boundary conditions - to cause the only other non-vanishing error in the potential energy surfaces. The influence from this factor is therefore also monitored by calculating potential energy surfaces for two different boundary conditions. In Fig. 1 we show the geometries of the H,/Al(llO) system used for the mapping of the potential energy surfaces (PESs). Prior work has shown [13,27] that the minimum energy reaction path (defined as a steepest descent from the smallest possible dissociation barrier) with either of the exchange-correlation approximations is the one where the molecule lies parallel to the surface, centered over and aligned along the so-called long bridge of the (110) surface. In determining this path no assumptions of certain symmetries were made. The relaxation of the surface was neglected, as the surface atoms are much heavier
K. Gundersen et al./Surface Science304 (1994) 131-144
than the H atoms and hence do not have time to move appreciably during the dissociation event. The minimum energy reaction path only involves two ionic degrees of freedom, the distance, 2, of the molecular center of mass from the outermost (IlO)-layer of Al nuclei and the H-H bond-length, b. The reaction path is illustrated in Fig. la. The final state of the dissociation is that of the H atoms chemisorbed on-top of the outermost Al atoms. This position implies a high frequency of the vibration of H perpendicular to the surface. We calculate the frequency to be 1775 cm-’ (220 meV) in good agreement with the experimental value of 1810 cm-’ [28]. The potential energy surfaces are mapped by performing a number of total energy calculations for fixed ionic coordinates. Ten different b-values combined with thirteen different Z-values where chosen for the mapping. Contour plots based on these total energies are shown in Fig. 2. Figs. 2a and 2b respectively show the LDA and GGA
2.25 2.00 1.75 03
1.50
N
1.25 1.oo 0.75 2.25 2.00 1.75
03
1.50
N
1.25 1.oo 0.75
1.0
1.5
2.0
2.5
3.0
Fig. 2. Contour plot of the Z-b plane potential energy surfaces. (a) LDA potential for the supercell of Figs. la and lb, (b) GGA potential for the supercell of Figs. la and lb, (c) GGA potential for the supercell of Fig lc. The thick lines indicate the reaction paths. The contour lines run from 0.1 to 1.0 eV with intervals of 0.1 eV.
(a>
09
133
@>
Fig. 1. (a) Perspective view of the H, /AKIlO) system. Small spheres represent H atoms, large spheres Al atoms. The dark line is the minimum energy reaction path. (b) The (1 x2) surface supercell of an orthorhombic supercell. (c) The surface supercell of a monoclinic supercell.
based potential energy surface for the supercell of Figs. la and lb. Fig. 2c displays the GGA based potential for the supercell of Fig. lc. It is seen that all three potential energy surfaces of Fig. 2 exhibit very similar topologies, it is only the barrier heights that differ. For the LDA calculation we get a barrier height of 0.25 eV, for the GGA it is 0.54 eV and 0.70 eV for the two
K. Gundersen et al. /Surfnce
134
different supercells respectively. The physics behind the barrier increase from the LDA description to the GGA description has been addressed in Ref. [13]. We find that the increase in barrier is related to an increased energy cost due to orthogonalization of the surface metal electron states to the H, states which are more tightly bound than in the LDA when non-local exchange-correlation effects are included. The difference between the barriers as calculated within the GGA using the two different supercells of Figs. lb and lc is 0.16 eV, suggesting a substantial HZ-H, interaction, at least in the cell of Fig. lb. HFre the HZ-H, separation is 2.8 A, while it is 4.9 A in the cell of Fig. lc. While we cannot be certain that the inter-molecular interactions are negligible even in the last cell, we suggest this result is close to the zero-coverage limit. Further test of this would require larger unit cells, which is not possible computationally at the moment, or a method like the one suggested by Feibelman [14]. The vibrational energies at the transition state have been determined through a normal mode analysis. The total energy was fitted to a quadratic function around the saddle point. The fit was done over energies typical to the ground-state motion. In case the directions of the modes were not given by symmetry a proper diagonalization of the dynamical matrix was performed. In Table 1 the result is given. The frequencies are seen to
Table 1 The vibrational
frequencies
of the normal
modes
at the transition
Science 304 (1994) 131-144
be quite independent of whether the LDA or the GGA based potentials are used. We find only small differences between the two geometries of Figs. lb and lc with respect to the frequencies. All the frequencies of the normal modes are sizable. This means that the barrier for dissociation depends strongly on the other four degrees of freedom and therefore that they should be included in the investigation of the dynamical behavior of H, on Al(110). The coming sections of this paper will focus on how to include all six dimensions without having to map the potential energy surface entirely and perform a full 6D dynamic simulation. Details of the total energy calculations now follow - the reader might prefer to skip directly to the next section for the description of the model of sticking. Aluminum slabs of five (110) layer thickness separated by five layers of vacuum are used. The theoretical lattice constant, 3.96 A, is used. H, molecules approach both sides of the slabs and the final coverage after dissociation is one H per surface Al for the supercells of Fig. lb as well as for Fig. lc. As stated, the total energy calculations are performed within either the LDA (Ceperley-Alder [29]) or the GGA (The GGA-II version from 1991 by Perdew and Wang - see Ref. [21] and references [26,27] therein). Aluminium is represented by a pseudopotential of the Kerker type [30] in the Kleinman-Bylander form [31]. This ab
states of the minimum
energy
paths
calculated
within
the LDA
and the GGA
Reaction
path
Mode in the 6-Z
plane
Frustrated translation in the Y-8 plane Frustrated rotation in the Y-8 plane Rotation parallel to the surface, 4 Translation
along the close-packed
a) Calculated with b, Calculated with All numbers are coordinate for the
rows, X
the supercell of Fig. lb. the supercell of Fig. lc. in meV. Y is the coordinate molecule.
LDA
GGA a’
GGA b,
ho*
1lOi
1lOi
1lOi
gw, ho, gw, h%
120 80 180 70
110 60 170 80
110 _
h%
40
70
perpendicular
to the close-packed
rows.
_ _
19 is the out-of-surface-plane-rotation
K. Gundersen et al. /Surface
initio non-local pseudopotential was generated within the LDA. As the adsorption properties of the surface are expected not to depend on the aluminium core levels which are most affected by the non-local exchange-correlation effects, we use this potential for the GGA study too. This procedure is justified by work by Garcia et al. [26]. The hydrogen atoms are represented by the bare Coulomb potential and all non-local effects are therefore included here. The bare Coulomb potential naturally sets a rather high demand for the basis set required to represent the electronic states. We choose to include planewaves with a kinetic energy up to a cutoff of 500 eV and find that this causes an absolute error less than 0.02 eV/molecule on the barrier heights. As the GGA requires quantities like 1V,,I/n one might be concerned that a rather dense FFT grid would be required to obtain convergence of the exchange-co~elation energy. However, direct tests have shown that the FFT grid which includes all G - G’ vectors of basis set planewave vectors, G and G’, suffices to converge the exchange-correlation energy to less than 1 meV per molecule. We choose k-point grids according to Monkhorst and Pack [32], with (4, 8) and (6, 6) k-point grids for the surface Brillouin zones of the supercells of Figs. lb and lc respectively. The error in the total energies caused by this choice of k-point grids is negligible compared to the error due to the finite basis set. The solution of the Kohn-Sham equations is found by means of an all-bands minimization technique [27] representing a slight generalization of the band-by-band minimization technique of Teter et al. [33]. Occupation numbers for these solutions are found by the Gaussian smearing method of Fu and Ho [34] as modified by Needs et al. [35]. The smearing width is 200 meV. All calculations both within the LDA and the GGA are iterated to self-consistency. In order to improve the stabili~ of the calculations involving the GGA, all gradient dependent terms in the exchange-correlation evaluation are neglected for densities smaller than 0.00035 A-“. This is done purely for numerical reasons and represents no further approximation.
Science 304 (1994) 131-144
135
3. The model of sticking The dynamical problem of H, dissociating on metal surfaces concerns six ionic degrees of freedom (neglecting all the degrees of freedom of the surface ions). In the particular case of H, dissociation on Al(llO), where there is a sizable corrugation of the PES in all of the six dimensions [13], all of these have to be included in a modefing of the sticking coefficient [12]. Due to the small mass of the H, the dynamical problem preferentially should be studied quantum mechanically. Up to now this has not been possible due to ~mputational constraints [36]. Therefore some approximations to the six-dimensional quantum dynamical treatment must be made. We describe below an extension of the hole mode1 of Karikorpi et al. [20], that includes a full quantum dynamical description in two dimensions while only estimating the dynamical effects of the remaining four dimensions. The two degrees of freedom that are treated in a full quantum calculation are the HZ-surface distance, 2, and the H-H bond length, b. As noted above, these degrees of freedom contain the entire minimum energy reaction path. Also contained is the vibrational dimension of the gas phase molecule. This allows for a proper representation of vibrationally ground state and excited state molecules, which is needed in the discussion of the coupling of the sticking to the vibrational state. These “first” two degrees of freedom are highly coupled as can be judged from the contour plots of Fig. 2. However, a change of coordinates into one along the reaction path and one perpendicular to it, will allow for a somewhat simpler representation of the PES of these two dimensions. Such a representation is given in Fig. 3, where the thin line is the potential energy along the reaction path, R. The vibrational energy of the dissociating molecule can approximately be described as the energy levels in a harmonic fit to the potential energy in the coordinate perpendicular to the reaction path. Such vibrational energy levels along the reaction path have been determined for a molecule in the ground state, v = 0, and in the first excited state, Y = 1. In Fig. 3 these are added to the potential
136
K. ~~der~e~
et al. /Surface
0.8
r 2 0.6 6 g 0.4
w 0.2 0.0
0.8 F 5 0.6 $ Gi 0.4 5 0.2
O.OL’. : -2.0 -1.0
,/ ’
1
’
’
0.0 1.0 R [Al
’ 2.0
’
3.0
Fig. 3. The potential energy per H, molecule (thin line) along the minimum energy path, the effective ID PES of a molecule in the vibrational ground state (thick line) and of a molecule in the first vibrational state (dotted line). The vibrational energy is determined at the levels, $Aw(R) and $ho(R), in a harmonic fit to the potential energy perpendicular to the reaction path, R. The latter is defined as the sum of (AR1 in A, where R is the six-dimensional vector of the Cartesian Soordinates of the two hydrogen atoms. R is zero for Z-2.3 A. wa is the intra-molecular vibrational frequency of the gas phase H, molecule. Each subfigure refers to a different PES: (a) the LDA based PES of Fig. Za, (b) the GGA based PES of Fig. Zb, Cc>the GGA based PES of Fig. 2~.
energy along the reaction path and provide adiabatic, quantum mechanical estimates of the effective 1D potential energy surfaces (thick and dotted lines) [16]. The effective 1D PESs illustrate
Science 304 (19941 131-144
the “bootstrapping” effect [15]; as the vibrational frequency of the molecule is softened as a result of the molecule-surface interaction, a molecule will encounter an effectively smaller barrier for dissociation. The remaining four degrees of freedom, (the two molecular rotations and the two translations parallel to the surface), we now have to model approximately. The first approximation is to consider these degrees of freedom individually, which implies that they only couple weakly to the first two, and that they couple only weakly with one another. The second approximation is to neglect any velocities in these dimensions. This means that we model only molecules impinging normal to the surface having no rotational energy (J = m = 0). Next we have to choose how to model the dynamical behavior in these dimensions. Candidates for simple models are the extreme cases, i.e. either the adiabatic or sudden limit. And for each of these in either a quantum mechanical or a classical picture. As the molecule approaches the surface the potential in the remaining four degrees of freedom increases from zero frequency to some finite value. If, on the timescale of the vibrational period, the frequency increase is small compared with the frequency itself, then we are in the adiabatic limit. In a quantum mechanical picture, this adiabatic limit implies that the molecule has time to “build up” the zero-point energy, ihw;, in each of the four degrees of freedom. The energy has to come from the kinetic or the initial intra-molecular vibrational energy of the scattering molecule, as illustrated in Fig. 3 in the Z-b plane. A similar energy “build up” takes place in a classical picture. Here the vibrational energy, E, divided by the frequency, o, for each of the additional degrees of freedom is an adiabatic i~~~~ia~t 1371. Thus when the frequency increases the energy “bound” in the corresponding mode increases proportionally to the frequency - like in the quantum mechanical picture. The other extreme is the sudden limit. If only a fraction of a vibrational period can be completed, while the frequency change takes place, the molecule does not have time to react on the change in the potential. The molecule will there-
K Gundersen et al. /Surface
fore interact with unchanged coordinates in the third (and higher) degrees of freedom. Classically, this implies that the molecule encounters a surface with a distribution of barrier heights. The effect of the extra degrees of freedom must therefore be included by integrating out the coordinates while weighting the sticking probability for each set of coordinates with the likelihood for the molecule of having these coordinates. For the purpose of the present sticking model, we adapt a classical description within the sudden approximation. That this is adequate can be judged from the following argument: The typical kinetic energies of impinging molecules are about 0.4 eV (for those that stick) and the potential strength (vibrational frequency) in each of the four degrees of freedom rises towards the transition state to values in the range 40-180 meV as given in Table 1 of the previous section. This gives of the order a few tenths of a period of vibration for a classical trajectory in the time the molecule spends in the region of appreciable molecule-surface interaction before it reaches the transition state. Now for the details of the hole model. We start out by considering the sticking coefficient, S,,(E), in the two dimensions of the Hz-surface distance, Z, and the H-H bond length, b. The S,,(E) is chosen to be the sticking coefficient resulting from a 2D quantum dynamical simulation. This we do as a split-time-propagation of a wavepacket [38] in the 2D potential energy surfaces of molecule-surface interaction described in the previous section. The wavepacket is initiated at some large Zvalue. First, the quantum mechanical solution of one vibrational state, v, is found in the 1D intersection of the PES at this large Z-value. This is done as a propagation in imaginary time. If a vibrationally excited state is desired a proper out-projection of lower lying states is done during this initial propagation. Next, the solution is given a Gaussian shaped kinetic energy, E, of width 5-15% of E in the Z-direction. The following 2D propagation is then performed until the products have reached the asymptotic regions of space. In the exit channel there is an absorbing potential making sure that we are not disturbed by back
137
Science 304 (1994) 131-144
reflection due to the periodic boundary conditions that are imposed for simplicity. The sticking coefficient is given by the fraction of the wavepacket that moves down the exit channel. Tests have been carried out to ensure that the sticking coefficient is well converged with respect to the finite time step of the wavepacket propagation and with respect to the finite phase space grid-density on which the potential and the wavepacket are represented. If the additional four ionic degrees of freedom, 13,4, X and Y, where all steered such that the lowest possible barrier was encountered by all impinging molecules then the 2D sticking coefficient, S,,(E), would also be the 6D sticking coefficient. However as argued above, for the kinetic energies of interest we are in the sudden limit, where the instantaneous ~9,4, X and Y-coordinates of the gas (or molecular beam) phase remain essentially unchanged during the reaction event. Therefore to get from S,,(E) to the sticking coefficient, S,,(E), of all six ionic degrees of freedom we now subtract from the energy argument of S,,(E) whatever is needed to overcome the barrier increase, AE,(8, 4, X, Y), due to the values of 13, 4, X and Y. The S&E AE,(8, 4, X, Y)) is then averaged over all configuration space to yield S,,(E). The barrier increase, AE,(B, 4, X, Y), we estimate from a decomposition of the considered four ionic degrees of freedom into normal mode coordinates, ui(8, 4, X, Y>, at the transition state. Each of these normal modes, ui, is associated with a frequency, oi, that describes to second order the behavior of the potential. The wi were determined in the previous section. The resulting barrier increase reads
i=3 (1)
where mu is the proton mass (assuming u properly normalized). We notice that the expansion above will be best for energies not too far from the lowest activation energy. The present description is therefore best suited for describing the onset of sticking. Now only the average over
K. Gundersen et al. /Surface Science 304 (1994) 131-144
138
configuration space remains. This is done by calculating the volume of the “hole” in configuration space, on,(E), available for sticking at kinetic energy E and properly normalizing by the volume of the entire configuration space, 0,. The volume of the “hole” becomes a,(E)
= 2j$,,(E
the final 6D sticking curves and compare with experiment, it is very instructive to consider separately the output from the 2D wavepacket simulations. Also some support will be given to the hole model in Section 4.2, concentrating on a comparison between 3D wavepacket simulations and 3D hole models.
- A&,(0, $7 X, Y)) 4.1. Sticking curves for 20 PES
X
sin 0 d0 d+ dX dY,
(2) where the integration is over the 4D configuration space and sin 0 is the Jacobian with 8 being the polar angle. The factor of two in front of the integral is due to the inversion symmetry of the molecule which means that there are actually two minimum energy reaction paths in all configuration space whereas AE, in Eq. (1) is expanded as if only one minimum energy reaction path was present. The total volume of possible configurations is: %=j
sin 0 dtl d+ dX dY= 47rA,,
(3)
4D
where A, is the area of the (1 x 1) surface cell. The sticking probability in this 6D hole model is finally &n(E)
=%(E)/%*
(4)
Here the hole model has been described for all six dimensions. It can, however, easily be modified to yield a 3D sticking curve, where besides b and 2 only 1 degree of freedom is taken into account. If - as will be the case below - the X-coordinate is considered, the sticking coefficient is given by S,,(E)
= ;
jS,,(
E - $2muw;X”)
dX,
(5)
0
where R, is now the length of the (1 X 1) surface cell in the X-direction.
4. Results and discussion In this section we present and discuss the sticking curves that emerge from the sticking model of the last section, when applied to the PESs of Section 2. However, before we present
Much understanding of the H,/Al(llO) system can be gained already at the level of the first part of the model, the 2D quantum dynamical simulation of the reaction event. In Fig. 4 are shown the 2D sticking curves for the PES of Fig. 2. For each PES is shown the sticking curve for an impinging molecule, which is in its vibrationally ground state, v = 0, and in its vibrationally first excited state, v = 1. The sticking curves for the LDA based PES (lowest barrier) is shown in Fig. 4a. The sticking of vibrational ground state molecules is seen in general to increase with increasing energy. However a pronounced suppression of the sticking takes place in the kinetic energy range from 0.1 to 0.3 eV. This suppression can be understood in terms of the impinging molecules being trapped in the metastable v = 1 states in front of the barrier (see Fig. 3a). The trapping of the molecules in such metastable states strongly affects the sticking behavior and may cause the observed decreased sticking probability. The effect is discussed in detail by Halstead and Holloway [39]. That the effect is observed in the right energy range can be judged from the v = 1 curve of Fig. 3a. The bottom of this dip corresponds to a kinetic energy, 0.09 eV, of a ground state H, and the curvature of the dip along the reaction path coordinate corresponds to a vibrational frequency of 0.25 eV. A metastable level is therefore expected for kinetic energies around 0.22 eV. The results of the 2D dynamical simulations on the two GGA PESs are shown in Figs. 4b and 4c. Again the sticking probability increases with the kinetic energy for small kinetic energies. For larger kinetic energies the curves bend over and the sticking decreases with the kinetic energy. This is the same effect of interaction with
K. Gundersen et al. /Surface
1.0
i$
I
’
’
’
I
’
’
I
0.8
(d
% 0.6 ii 2 0.4 Ii .g 0.2 c/J 0.0
1.0
I
’
g 0.8 tu % 0.6 2 p 0.4 E .g 0.2 cn 0.0
139
H, in the Z-b plane is considered at different kinetic energies. The probability distributions are taken from the wavepacket simulations on the GGA based potential of Fig. 2b. At high energies the reflected beam shows excited H-H vibrations in the interaction region. Unlike the LDA based results, the GGA potentials with a high barrier for dissociation, give rise to differences between the vibrational ground state (v = 0) and the first excited state (v = 1) sticking at small energies; the v = 1 sticking curves has an 0.15 eV lower onset. We will return to this considerable coupling between the vibrational state of the wavepacket and the sticking coefficient when comparing with experiment below. 4.2. Verification of hole model - 30 study
1.0 $ 0.8 Cu “0 0.6 lc p 0.4 P .g 0.2 cn 0.0
Science 304 (1994) 131-144
0.0
0.2
0.4
0.6
0.8
Kinetic Energy [eV] Fig. 4. 2D sticking curves for the ground state (solid line) and the vibrationally first excited state (dotted line). Also shown is the 3D wavepacket simulated sticking curve for the ground state (dashed line). Each subfigure refers to a different PES: (a) the LDA based PES of Fig. 2a, (b) the GGA based PES of Fig. 2b, (c) the GGA based PES of Fig. 2c.
metastable states in front of the barrier as for the LDA PES, but now shifted to higher energies due to the more repulsive character of the two GGA PESs. The trapping in the metastable states implies that impinging molecules may become vibrationally excited due to the interaction with the surface. An indication that the interaction is indeed strong enough for this to happen is given in Fig. 5. Here the probability distribution for the
We have investigated the accuracy of the hole model by comparing 3D quantum calculations with the hole model in its 3D version, treating two degrees of freedom quantum dynamically and one degree of freedom in the classical sudden approximation through Eq. (5). The 3D potential has been determined within the LDA description. The 3D GGA potentials have been obtained by adding the difference between the 3D and 2D LDA potentials to the 2D GGA potentials. This means that the 3D potentials used are, to some extent, no more than model potentials. Since they are only used to ascertain the general applicability of the hole model this introduces no further approximations. The third dimension is chosen to be the translation along the close-packed rows, a very smooth mode which remains a normal coordinate all the way out from the surface. The vibrational frequency in this mode is 40 meV. In order to address the range of applicability of the hole model we also modify the mode to yield a vibrational frequency of 180 meV. Hereby the whole range of potential strengths of the four modes to be included is covered. The 3D quantum calculations are performed in a similar manner to the 2D calculations. The initial wavepacket is chosen to be the previously mentioned 2D wavepacket multiplied by a planewave in the third dimension.
140
(a) At barrier, E=0.2 eV
(b) Reflected wave, E=0.2 eV
(c) At barrier, E=0.4 eV
(d) Reflected wave, E=0.4 eV
(e) At barrier, E=O.6 eV
(f) Reflected wave, E=0.6 eV
Fig. 5. Probability density in the Z-b plane at different incoming energies. The PES is that of Fig. 2b. (a), (b) The scattering wavepacket for E = 0.2 eV. It is seen that the wavepacket is adiabaticaIly reflected from the potential with only very little tunneling. Cc), (d) The scattering wavepacket for E = 0.4 eV. There is a large barrier penetration and no vibrational excitation of the beam in the reaction zone. (e), (f) The scattering wavepacket for E = 0.6 eV. The wavepacket is clearly excited in the interaction region of the potential. Only a small barrier penetration is seen. Notice the excited tail of the outgoing wavepacket indicating that it has been trapped in an excited metastable state.
K. Gundersen et al. /Surface
In Fig. 6 are given the quantum dynamical 2D and 3D sticking curves. The sticking curves are seen to be very sensitive to the added third dimension. This confirms that all of the six molecular degrees of freedom must be included in a proper description of the sticking coefficient. Also included in Fig. 6 are the 3D hole model estimates of the sticking curve. It can be seen that the hole model captures the large decrease of the sticking probability due to the third degree of freedom, and that it is almost quantitative at the lowest energies (up to about 0.4 eV) where it is supposed to work best. The deviations are largest
141
Science 304 (1994) I31 -144
for the high frequency mode, the reason being that here the steering of the molecule towards the lowest energy barrier is strongest. As a consequence, the hole model systematically underestimates the sticking probability for this mode. The slight underestimate for the low frequency mode is, on the other hand, most likely due primarily to the second-order expansion of Eq. (1). 4.3. Results of the 60 hole model In Fig. 7 we show the 6D sticking curves obtained from our 2D quantum curves, using the
r
0.60
I
I
(b)
’
I
I
2D
1 0
0.50
i t .= h =
I
0.40
I-
% s z
I 0.30
i-
.-F
i
% z
0.20 0.005 0.10
Kin. Energy [eV] Fig. 6. The hole model strengths of the potential
applied to 3D sticking in the third dimension.
Kin. Energy [eV]
curves plotted together with the quantum (a) PES of Fig. 2b, (b) PES of Fig. 2c.
3D sticking
curve
for two different
K Gundersen et al. /Surface
Kinetic Energy [eV] Fig. 7. 6D sticking curves obtained from Eqs. (l)-(4) using the v = 0 ZD sticking curves of Fig. 4, (a) corresponds to the LDA PES while (b) and (c) correspond to the two GGA PESs (Figs. 2b and 2c respectively). The difference between (a) and (b) is caused by the inclusion of the non-local exchange-correlation effects, (b) and (cf differ due to the different Ha-H, separations with Cc) being the most converged result.
hole model, Eqs. (l)-(4). Comparing to the 2D sticking curves of Fig. 4 it is seen that the inclusion of all six degrees of freedom is essentia1 for a proper description of the sticking process. The large difference between the LDA and GGA results is striking. The GGA includes, in an approximate way, some non-local exchange-correlation effects not included in the LDA. One would therefore expect the GGA to be more accurate, but there is no reason to believe that it gives a
2
0.0008
r----
i-’
/Level
d
-7
of atomic H
-~E(T~YT~)I
in the beam ”; g
0.0006 0.0004
1 L
2 %
0.0002
L
g g
0.0000
i 0.0 Kinetic Energy [eV]
Fig. 8. The pure and seeded beam calculated from the 2D quantum curves of Fig. 4c using the hole model, Eq. (61, including the v-dependence of the sticking up to v = 2. The experimental result of Berger et al. [3] for the pure beam is also included. Also shown is the level of atomically adsorbed hydrogen found in the experiment. This has been subtracted from the raw experimental data to give the molecular sticking probability and thus gives an indication of the uncertainty.
Science 304 (19941 131-144
quantitative description. The difference between the LDA and GGA curves thus gives an indication of the accuracy by which “ab initio” sticking probabilities can be determined today. In order to compare the calculated sticking curves with the results of the beam experiments of Berger et al. [3] we must include the effect of vibrational excitations in the beam. In a beam experiment, both the translational energy and the degree of vibrationa excitation of the molecules in the beam, are given by the temperature of the nozzle. If the kinetic energy of the beam is changed by seeding techniques, the nozzle temperature and the beam energy can be considered independent variables, and the sticking probability as a function of normal kinetic energy E and nozzle temperature T,, is given by S( E, T,) = [ S,(E) + exp
X[I+exp[-$j
j-$hj”‘(“) +...I +...]-I,
(6)
where wO is the vibrational frequency for H, and S,(E) is the (6D) sticking probability - from Eq. (4) - for a molecule in vibrational state Y at the kinetic energy E. These have been summed with their Boltzman weights at the nozzle temperature to give the total sticking probabili~. For an unseeded beam, the relation between the translational energy and the nozzle temperature is E(T,)
= ;kBTn
(7)
and the sticking probability of the pure beam is therefore given by S(E(T,), T,). The model calculation, using what we believe to be our best results based on the Iargest HZ-H, separation possible and the GGA, is compared to the experimental result in Fig. 8. Considering the uncertainties in the tota energy calculation and in the description of the dynamics and the uncertainties in the experiments (see below) the agreement between theory and experiments is surprisingly good. In Fig. 8 we have also shown the sticking curve for a seeded beam with a fixed nozzIe tempera-
K. Gundersen et al. /Surface
ture of T,, = 1850 K calculated according to Eq. (6). Here we have included up to the second excited vibrational state. Higher excited states do not contribute, on the scale of interest, because their Boltzman factors are less than 10p4. It is clear that there is a difference between the pure and the seeded beam sticking curves but it is small on the scale of the experimental accuracy. The large contribution from atomically adsorbed hydrogen drowns the effect from the excited states and this is thus not seen in the experiment. It is therefore not surprising that Berger et al. suggest that there is no vibrational coupling for the H,/Al(llO) system [3]. In Ref. [13] we used a simplified version of the hole model to get a sticking probability from the LDA potential. This calculation differs from the present one in three respects. First, it was based on the LDA potential. Second, it had only a rough WKB estimate of the quantum dynamics in the Z-b plane, which among other things did not include the coupling between the H-H vibration and the sticking. And finally, it did not include the thermal averaging of Eq. (6). The present result is therefore substantially different from the earlier version.
5. Conclusion In the present paper we have presented a set of ab initio calculations of the H, potential outside an Al(110) surface. Further we have modeled the fully six-dimensional H, dynamics on this potential energy surface. We have shown that the treatment of exchange correlation effects greatly influences the dynamics. We have shown that a full quantum treatment of the dynamics in the Z-b plane combined with the hole model gives a good description of the dynamics, at least for the onset of sticking. Finally, we have compared the hole model calculation of the sticking probability directly to the beam experiments of Berger et al. [3] and shown that a reasonable agreement between theory and experiment can be obtained considering both the theoretical approximations and the experimental uncertainties.
Science 304 (1994) 131-144
143
6. Acknowledgements We are grateful to M.C. Payne for giving us the CASTEP (Cambridge Serial Total Energy Package) code. We are indebted to Stephen Holloway for pointing out an initial misinterpretation of the 2D sticking curves. We thank the Center for Applications of Parallel Computers (CAP) for providing computer resources. The research was done in part under the Center for Surface Reactions of the Danish Research Councils and in part under the Center for Atomic Scale Materials Physics (CAMP).
7. References 111 B.E. Hayden and C.L.A. Lamont, Phys. Rev. Lett. 63 (1989) 2823. 121 H.F. Berger, M. Leisch, A. Winkler and K.D. Rendulic, Chem. Phys. Lett. 175 (1990) 425. [3] H.F. Berger and K.D. Rendulic, Surf. Sci. 253 (1991) 325. [4] H.A. Michelsen and D.J. Auerbach, J. Chem. Phys. 94 (1991) 7502. 151 C.T. Rettner, D.J. Auerbach and H.A. Michelsen, Phys. Rev. Lett. 68 (1992) 1164. [6] H. Hjelmberg, J.K. Norskov and B.I. Lundqvist, Phys. Ser. 20 (1979) 192; P. Johansson, Surf. Sci. 104 (1981) 510. 171 J.K. Nsrskov, A. Houmoller, P. Johansson and B.I. Lundqvist, Phys. Rev. Lett. 46 (1981) 257. 181 J. Harris and S. Andersson, Phys. Rev. Lett. 55 (1985) 1583. [91 J.E. Miiller, Phys. Rev. Lett. 59 (19871 2943. [lOI J.E. Miiller, Surf. Sci. 272 (1992) 45. [ill J.K. Norskov, J. Chem. Phys. 90 (1989) 7461. 1121 C. Engdahl, B.I. Lundqvist, I-J. Nielsen and J.K. Norskov, Phys. Rev. B 45 (1992) 11362. [131 B. Hammer, K.W. Jacobsen and J.K. Norskov, Phys. Rev. Lett. 69 (1992) 1971. [14] P.J. Feibelman, Phys. Rev. Lett. 67 (1991) 461. 1151 J. Harris, Surf. Sci. 221 (1989) 335. I161 M.R. Hand and S. Holloway, J. Chem. Phys. 91 (1989) 7209. [171 S. Kuchenhoff, W. Brenig and Y. Chiba, Surf. Sci. 245 (1991) 389. [181 M.C. Payne, M.P. Teter, D.C. Allan, T.A. Arias and J.D. Joannopoulos, Rev. Mod. Phys. 64 (1992) 1045. I191 B. Hammer, K.W. Jacobsen and J.K. Norskov, Phys. Rev. Lett. 70 (1993) 3971. [201 M. Karikorpi, S. Holloway, N. Henriksen and J.K. Norskov, Surf. Sci. 179 (1987) L41. [211 J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson,
144
K. Gundersen et al. /Smface
M.R. Pederson, D.J. Singh and C. Fiolhais, Phys. Rev. B 46 (1992) 6671. [22] G. Ortiz and P. Ballone, Phys. Rev. B 43 (1991) 6376. [23] X.J. Kong, CT. Chart, K.M. Ho and Y.Y. Ye, Phys. Rev. B 42 (1990) 9357. [24] G. Ortiz, Phys. Rev. B 45 (1992) 11328. 1251 M. Kiirling and J. ~glund, Phys. Rev. B 45 (1992) 13293. [26] A. Garcia, C. Elsiisser, J. Zhu, S.G. Louie and M.L. Cohen, Phys. Rev. B 46 (1992) 9829. [271 B. Hammer, PhD Thesis, Physics Department, Technical University of Denmark (1993). [28] H. Kondoh, M. Hara, K. Domen and H. Nozoye, Surf. Sci. 287/288 (1993) 74. [29] D.M. Ceperly and B.J. Alder, Phys. Rev. Lett. 45 (1980) 566; J.P. Perdew and A. Zunger, Phys. Rev. B 23 (1981) 5048; J.P. Perdew and Y. Wang, Phys. Rev. B 45 (1992) 13244. [30] G.P. Kerker, J. Phys. C 13 (1980) L189.
Science 304 (1994) 131-144
[31] L. Kleinman and D.M. Bylander, Phys. Rev. Lett. 48 (1982) 1425. 1321H.J. Monkhorst and J.D. Pack, Phys. Rev. B 13 (1976) 5188. 1331M.P. Teter, M.C. Payne and D.C. Allan, Phys. Rev. B 40 (1989) 12255. [34] C.-L. Fu and K.-M. Ho, Phys. Rev. B 28 (1983) 5480. [35] R.J. Needs, R.M. Martin and O.H. Nielsen, Phys. Rev. B 33 (1986) 3778. [36] U. Nielsen, D. Halstead, S. Holloway and J.K. Nsrskov, J. Chem. Phys. 93 (1990) 2879. [37] L.D. Landau and E.M. Lifshitz, Quantum Mechanics (Pergamon, Oxford, UK 1977). [SS] J.A. Fleck Jr., J.R. Morris and M.D. Feit, Appl. Phys. 10 (1976) 129. [39] D. Halstead and S. Holloway, I. Chem. Phys. 93 (1990) 2859.