Theoretical studies of adatom diffusion on metal surfaces

Physics Letters A 180 ( 1993 ) 444-452 North-Holland

PHYSICS LETTERS A

Theoretical studies of adatom diffusion on metal surfaces Keh-Dong Shiang Institute of Physics, Academia Sinica, Nankang, Taipei 11529, Taiwan, R O("

Received 30 June 1993; accepted for publication 25 July 1993 Communicated by V.M. Agranovich

We propose in this paper a theoretical model to investigate surface self-diffusionof singleadatoms on two different low-index planes, closelypacked (001) and densely packed ( 111 ), of face-centered-cubicrhodium, nickel and copper metal crystals. Two realistic model potentials are applied to describe the interatomic interaction of the adatom-substrate systems. The first model is a Morse-type potential, which involves several empirical fittings of bulk properties of solid. The second, newlypopular, potential was introduced by Sutton and Chert, which incorporates many-body effects. With these potentials, conventional molecular dynamics (MD) is employed to obtain trajectories of the atoms. The average squared displacements are computed for a range of initial kinetic energies, and the surface diffusion constants can be obtained by means of the Einstein relation. The estimated random walk exponential prefactors and activation energies exhibit an Arrhenius behavior, and are compared with previous results.

1. Introduction The migration of individual atoms self-adsorbed on their own metal surfaces is a f u n d a m e n t a l aspect of many surface p h e n o m e n a , such as adsorption (physisorption and chemisorption ), desorption, epitaxial growth of thin films, crystal growth, corrosion, surface reconstruction, and surface enhanced chemical reactions. In particular, we restrict our focus on surface diffusion [ 1-4 ] which involves many important mechanisms which are worth studying. Till now experimental and theoretical attention has been paid to the spectroscopy and dynamics of adsorbed atoms. A major experimental tool for investigating surface diffusion is the field ion microscope ( F I M ) [ 5 - 8 ] , since it is capable of visualizing individual metal atoms on metal substrates. In general, FIM experiments determine the diffusion constant of a single adatom at several temperatures and fit the results to a general Arrhenius form. Experimental measurements of the self-diffusion have been primarily restricted to a n u m b e r of appropriate metals due to the high field used for imaging: p l a t i n u m [9,10], rhodium [11], nickel [12], tungsten [ 1 3 - 1 8 ] , and iridium [ 1 9 - 2 7 ] . These different systems were also examined by Ehrlich [28]. In addition, a n u m b e r of 444

theoretical efforts were also performed by Tully et al. [29], Liu et al. [30], McDowell and Doll [ 3 1 - 3 4 ] , Park and Bowman [35], Feibelman [36] and Kellogg [37]. The choice of a reliable potential for a given metallic system is a difficult task due to the lack of knowledge of the interaction forces in real materials. In general, pair potentials [ 38 ] have been popular to describe metallic b o n d i n g because of their computational simplicity. Recently, a significant advancement in the empirical description of the interatomic potential has been to introduce a many-body (cohesive) term in addition to a pairwise potential [3942]. The first term models the effect of the local electronic density. The second term is a two-body (repulsive) potential, which can be fitted as a Morse function or an exponential (repulsive) form. In this work two model potentials, a two-body Morse potential and a many-body Sutton-Chert ( S - C ) potential, are applied to describe the interatomic forces in our systems. In general, the diffusion of an adsorbate on a surface can be characterized by two types of motion [ 1 ]. When the activation energy for diffusion, ta, is less than the thermal energy of the substrate, kBT, the adsorbates migrate freely across the surface. This type

0375-9601/93/$ 06.00 © 1993 Elsevier Science Publishers B.V. All rights reserved.

Volume 180, number 6

PHYSICS LETTERSA

of motion, which is called mobile diffusion, is characteristic of physisorbed species. When Va is significantly greater than kBT, the adsorbates are mainly localized at different surface lattice sites and move from site to site by hopping over the activation barrier. Diffusion via a hopping mechanism is characteristic of chemisorbed species. Due to the periodic arrangement of the substrate atoms, adsorbed atoms experience a periodic potential along the surface. Based on the hard-sphere model of the surface, an adatom is expected to jump along the direction of the least corrugation height. From the energy point of view, since hopping along the directions with the least corrugation height requires passing through atop sites or bridge sites, the binding energy is expected to be large and the barrier to diffusion correspondingly low. In the conventional picture of surface diffusion on metals, atom migration takes place by a series of displacements over the minimum in the potential barrier between adjacent binding sites. As a matter of fact, most experimental data are consistent with this assumption. Nevertheless, there are a few exceptions to this rule. An exception is that on the fcc(001 ) surface of some metals, Pt, Ir, Ni, A1, Cu and W, the surface diffusion might proceed by any other mechanism. The most striking feature is that the surface diffusion may occur by an exchange mechanism. An adatom might possibly take the place of a surface atom which then becomes the new adatom. In this study, however, we are not interested in this special mechanism. For many applications of surface diffusion static methods are quite acceptable. But in studies of hightemperature materials or materials in which there is an exceptionally high atomic mobility, it may become necessary to include thermal motions explicitly. This can be achieved by using the molecular dynamics technique, in which kinetic energy is included in the simulation. This method has been proved to be one of the most important tools in theoretical condensed matter physics. Thus we examine self-diffusion on the (001) and (111 ) faces of rhodium, nickel and copper surfaces from a microscopic perspective using molecular dynamics (MD) techniques. In this study, computer simulations are performed on a system of 648 classical atoms confined within a cubic box of fixed volume and subjected to periodic boundary conditions. Every atom in the cube

20 September 1993

has an image at all combinations of all translations of the parent cube along its edges. The time evolution of the ensemble is then followed by an iterative numerical solution of the equations of motion of the system. This approach involves the introduction of a many-body interaction potential and an investigation of the diffusion constants and activation energies. Comparison of these predicted results with the previous experimental and computed values will also be made. This paper is organized as follows: Section 2 contains the descriptions of the interaction potential models and of how molecular dynamics technique is applied to compute the diffusion parameters. The resuits are then summarized and discussed in section 3. A short conclusion is given in section 4.

2. Method of calculation

2. I. The interaction potential The dynamic processes occurring at a metal surface are governed by the adsorbate-substrate potential which has been a subject of great interest for a long time. Most MD studies of surfaces have been performed using pairwise potential forms such as the well-known Lennard-Jones type potentials. In this study, we select two potential models which are a semi-empirical two-body Morse-type and many-body Sutton-Chen potentials (i.e. Finnis-Sinclair potential ).

2. I. 1. Morse potentials The Morse function energy of interaction between atoms i and j separated by the distance r o is given by [38]

Uij(rij)

= D{exp [ - 2 a ( r o - r o ) ]

- 2 exp[ - a ( rij - r o ) ]) .

( 1)

The parameters, D, a, and ro, are determined by semiempirical fitting of the bulk lattice constant, cohesive energy, and compressibility.

2.1.2. Sutton-Chen potentials Before introducing the Sutton-Chen potential, a brief description of a related potential model, the 445

Volume 180, number 6

PHYSICS LETTERS A

embedded atom method (EAM), must first be given in this section. This popular technique of including many-body effects for metals has been developed recently by Daw and Baskes [39,40]. It is based on density functional theory, which asserts that the energy o f a solid can be written as a unique function of the electron density distribution. This is assumed to be the local density at each atomic site. Thus the electron density of the solid is approximated by the superposition of the atomic electron densities of the constituent atoms. In this approach the total energy of an arbitrary arrangement of atoms is given by

C;o.= E c',,

(2)

i

which is written as a sum of two terms,

U, = ½ ~ Oo(ro) +F,(p~) .

(3)

where p,= Z.g(r,j),

(4)

where/~o, is the total internal energy of an assembly of atoms, U, is the internal energy associated with atoms i, p, is the total electron density at atom i due to the rest of the atoms in the system, F,(p~) is the energy required to embed the atom i into the local electronic charge density p,, Oo(r,j) is the conventional two-body central potential between atom i and a t o m j separated by the distance ro, and~(r,j) is the contribution to the electron density at atom i due to atom j. The EAM functions are determined empirically by fitting to a number of measured properties of the solid such as the equilibrium lattice constants, elastic constants, bulk modulus, heat o f sublimation, and vacancy formation energy. Almost at the same time, Finnis and Sinclair [ 41 ] developed a model which is mathematically equivalent to the EAM, but has been deduced using the second-moment approximation to the density of states in the tight-binding theory. Based on their idea, Sutton and Chen [42] introduced a model potential with an a t o m - a t o m repulsive term and a term proportional to the square root of the local density to represent the cohesive energy o f the electrons. They chose F(p) to be pl/2, SO as to mimic the result of tight-binding theory in which f ( r ) could be interpreted as a sum of squares of overlap integrals. The 446

20 September 1993

basic equations of the Sutton and Chen approach are U,=~(½ ~ V(,'(j)-cp~/2),

(5)

]¢1

where the pairwise repulsive potential, V(r a), is modelled by ~'(r,j)

= (a/r,,)"

(6)

,

and the local density, p, is described as

p,= ~ (a/,',j)".

(7)

Here tSj is the separation between the atoms i and j, c is a positive dimensionless parameter, ~ is a parameter with the dimension of energy, a is the lattice constant, and m and n are positive integers with n > m. The choice of n determines the steepness of the repulsive potential and that of m determines the range. In practice, a better form of the Sutton-Chen potential may be gained by rearranging the above equations into the form

C',o. = ½~Z Z [(a/,',)" -c(p;-'/2 + p . . / 2 ) (a/r,;)"] .

(8)

which is similar to the Mie or Lennard-Jones potential.

2. 2. Dynamic calculations The evolution of the atoms in time and space is determined by the numerical solution of the classical equations of motion which can be integrated using the "'velocity" form of the Verlet algorithm. This standard method was clearly introduced by Allen and Tildesley [43]. Assuming atom i is at position R, at time t and has a velocity V,, then the position and velocity at time t+dt is given by

R ~ ( t + d t ) = R , ( t ) + I ~ ( t ) d t + ~1 F~(t) m dr2, V,(t+ ½dt) = I•(t) + 1 Fi(t) dr, 2

m

(9)

(10)

and l/](t+dt) = I~(t+ k d r ) + 21 k ] ( t + d t ) dr,

(11)

Volume 180, number 6

PHYSICS LETTERS A

where Fi is the total force exerted on the atom i of mass m by all the other atoms within the sphere with a cutoff radius. The computational procedures are explained as follows. The new positions at time t + dt are calculated using eq. (9), and the velocities at midstep are computed using eq. (10). The forces and accelerations at time t + d t are then computed, and the velocity move completed by applying eq. ( 11 ). When the Sutton-Chen potential is used in the dynamic calculations, two steps must be included in the computations. In the first step the local densities for each site are accumulated, and in the second step the energies and forces between each pair of atoms are computed. The molecular dynamics unit cell used for fcc (001) and fcc ( 111 ) are given in figs. 1 and 2, respectively. Each surface is arranged by stacking between eight or nine defect-free square or rectangular atomic layers, and each layer contains 81 or 72 atoms. Numerical simulations are performed using a conven-

20 September 1993

tional spherical cutoff, m i n i m u m image technique. The cutoff radius is chosen as 14 bohr in the present work. The initial velocity components obey a Maxwellian distribution in the form [44,45] P(Vix) = ( m / 2 7 t k T ) l / 2 e x p ( - m V 2 i x / 2 k a T ) ,

(12)

where P(Vix) is the probability density for velocity component vi~ of atom i in direction x, and similarly for the other two directions. Periodic boundary conditions are applied in the x and y directions parallel to the surface, but not in the z direction, where the motion is free. An adsorbed atom of the same species as the bulk is then added and allowed to diffuse over the surface. Since the vibration period in simple metals is about 1 ps, the time step for the simulation of this study is chosen to be 0.01 ps. We carry out a direct measurement of a diffusion coefficient by tracing the motion of individual adatoms over a time period of 200 ps, which includes three different types of M D simulations. Firstly, the quenching procedure is started from the unrelaxed positions o f atoms on

[001]

[111]

(

N

c3 "'~

X

Fig. 1. The top view and side view of the fcc(001 ) surface structure.

Fig. 2. The top view and side view of the fec( 111 ) surface structure. B is hcp site, C is fcc site. 447

Volume 180, number 6

PHYSICS LETTERS A

the surface. Then the surface system is allowed to relax to the m i n i m u m energy configuration. This is accomplished by setting the velocity of each atom to zero whenever the scalar product of the velocity and force becomes negative. This procedure is able to quench rapidly the system into the relaxed zerotemperature structure. After the relaxed structures of the surface are determined, a 10 ps constanttemperature (canonical) equilibrium run is made to ensure that the system is equilibrated at the desired temperature. During this procedure, atomic velocities are adjusted to correspond to the desired temperature every time step for the first 1000 steps. Afterwards, the system is left undisturbed in a 190 ps constant-energy (microcanonical) simulation.

2. 3. Determination of diffusion constants We now turn to the measurement of the diffusion constant. This is clone by tracking the motion o f individual adatoms on surfaces. We place a single adatom at a random position above the surface and follow the trajectories of the atoms for 200 ps. The selfdiffusion constants, D, can be extracted by examining the long-time behavior of the mean square displacement which for a diffusive process obeys the Einstein relation D= (Ar2(t)) 1 N 2dt - 2dt ~ [ r ' ( t ) - r ' ( 0 ) [ 2

(13)

where t is the time, r,(t) is the position of the ith adatom at time t, Ar2(/) is the time-dependent squared displacement of the adsorbate, N is the number of adatoms whose trajectories we traced, and d is the dimensionality of the diffusion space; d is two in the present case. The diffusion constants D are calculated assuming several different initial kinetic energies. In this approach, the temperature can be assigned by the ensemble-averaged value of each final kinetic energy.

3. Computational results and discussion It has to be noted that when an atom approaches a surface, it must disturb the local environment of the lattice. Interactions with the substrate cause the neighboring substrate atoms to relax from their nor448

20 September 1993

real lattice positions. This will lead to a m i n i m u m energy configuration. This is so-called "surface relaxation". A more complicated situation, "surface reconstruction", is the rearrangement of atoms at a crystal surface to a structure with a symmetry different from that of the underlying crystal planes. The temperature-dependent reconstruction also is an important aspect we may consider. However, we are now not interested in investigating these phenomena, but they may be worth studying in the future. Based on the result of surface relaxation, our simulation is prepared and carried out with this stable starting configuration, which has been presented in section 2.2. Roughly speaking, the position of the binding sites can be understood as the places where the adatom can most easily form bonds with a number of substrate atoms. The (001) surface of a face-centeredcubic lattice has square symmetry and each atom has coordination number eight with four nearestneighbors in the surface layer. The adsorption of atoms is generally at the fourfold hollow sites with a large binding energy. However, adsorption on an atop site is highly unlikely due to its unfavorable binding. An adsorbed atom resting in a fourfold equilibrium site (indicated as A) on the (001) plane is shown in fig. 1. As is evident from fig. 2, the (111 ) surface allows two-dimensional diffusion through a hexagonal geometry of two relatively shallow but distinct adsorption sites marked as B and C, corresponding to hcp and fcc stacking of the next layer, respectively. The popular Arrhenius expression for the adatom moving on the surface can be described as

D=Do e x p ( - t~/kBT) ,

(14)

where Do is the exponential prefactor, and }a is the activation energy. If we measure the logarithm of D versus the reciprocal of the surface temperature T, then the slope gives the activation energy and the intercept corresponds to the pre-exponential factor. These results are illustrated in figs. 3-8 as In D versus 1 / T and fall approximately along straight lines consistent with thermally activated processes. As Arrhenius behavior is expected from these lines, it is very straightforward to carry out the calculations of exponential prefactor Do and activation energy l'~. The diffusion characteristics, including the exponential prefactors and activation energies, are summa-

Volume 180, number 6

PHYSICS LETTERS A

-10

-7.6

20 September 1993

L. . . . . . . . . . . . . .

'

'

L

-11

tm

-

Rh on Rh(001)

-

~-

, I

....

Ni on N i ( l l l )

~ m8

-8 !

M_orsePoSentlial

-12

-

-13

-

-8.2

q

-14

-

-8.4

I

J T

4 -15

-8.6

-16

-8.8

100

200

150

250

300

i

E

350

~

150

1/T

200

250

300

1

I

350

400

450

500

liT

Fig. 3. Arrhenius plot of diffusion results for Rh(001 ). All the values are in atomic units.

Fig. 6. Arrhenius plot of diffusion results for N i ( 111 ). All the values are in atomic units.

-7.8 . . . . . . . . .

-7

I 7

"

Q ~ K _

Rh on Rh(ll 1)

Cu on Cu(001)

-8

-8.2

t21 cm

~,CX~,,

~

-8.4~ -8.6

." Morsepotential "'"C)~ S-C Potential

~~

- 7.5 D

~"~

~

-8

Morsepotential -~ntial

! 2

~

~-

-8.5

-8.8 ~-

d

~

J

-9

-9 100

150

200

250

300

350

400

250

300

350

liT

Fig. 4. Arrhenius plot of diffusion results for Rh( 111 ). All the values are in atomic units. -9 -10

~

,

,

.

.

.

.

450

500

Fig. 7. Arrhenius plot of diffusion results for Cu(001 ). All the values are in atomic units.

.

-8

[L k

-11

400 I/T

_

- 8.1

~

C u on .C_.u( 1 l I )

- 8.2

~

."

_.

Morsepotential

~ 3

q2

-8.4 -13

-8.5 -14 -15

~ L

~

---~

Morsepotential S-Cpotential

~

d

-8.6

"~±

a

-8.7

-16

:-

-8.8 150

200

250

300

350

400

450

500

l/T

Fig. 5. Arrhenius plot of diffusion results for Ni(001 ). All the values are in atomic units.

250

300

350

400

450

500

1/T Fig. 8. Arrhenius plot of diffusion results for Cu ( 111 ). All the values are in atomic units. 449

Volume 180, number 6

PHYSICS LETTERS A

an adjacent m i n i m u m , d is the dimensionality of the diffusion space, and n is the n u m b e r of j u m p directions available to the adatom. For an fcc(001 ) structure, d is two and n is four. For an fcc( 111 ) structure, d is two and n is three. By comparing eq. ( 14 ) with eq. (15), the well-known r a n d o m walk diffusion parameters thus are

rized in table I. As clearly displayed in this table, the calculations give the parameters which are in the correct quantitative order. However, the estimated activation energies for diffusion on the loosely packed (001) planes are higher than for diffusion on the close packed ( 111 ) planes. A more likely explanation is that close packing of ( 1 l 1 ) face enhances the mobility of surface adatoms. It is apparent that the adatom spends less time trapped in these shallow potential wells. On the other hand, the theoretically predicted exponential prefactors of r h o d i u m and copper are significantly higher than the experimentally d e t e r m i n e d values. Typical values of the exponential prefactor are in the neighborhood of 10-3 cm2/s for the (001) surfaces and of 10 -4 c m 2 / s for the ( l l l ) surfaces. G i v e n the energy of two configurations, ~ and V~, the diffusion constant in the harmonic approximation is written as [30,48] nva z

D= ~ff- exp[(~]- l~)/kaT]

~'~ = H - t~

( ]6 )

and rl v a z

Do-

2d

(17)

"

The attempt frequency of the simple harmonic oscillator was computed by v= ~

1

(c/M)

(18)

~/2,

where c is the force constant evaluated from the second derivative of the three-dimensional potential surface with respect to z at the m i n i m u m site of x and y. In principle, the value of u can be estimated via eq. (17) as the exponential prefactor is given. For instance, if the Do value of rhodium is assumed to be 10 -3 cm2/s, it will yield an attempt frequency value of the simple harmonic oscillator of around l 01 ~ s t which is in excellent agreement with the comm o n value of molecular vibration.

(15)

,

20 September 1993

where I~ is the potential energy with an adatom located at the lattice site, and t~ is the potential energy with the adatom located at the saddle point along the diffusion path. In this equation, u is the so-called "attempt frequency" which in the simplest approximation is just the harmonic frequency of the adatom in the potential m i n i m u m , a is the distance to

Table 1 Arrhenius parameters for atomic self-diffusion Surface

Face

Diffusion prefactor Do (cm2/s)

Aclivation energy 1~ (eV) present results

rhodium nickel

(001 (111 (001 (Ill

exp.

theor.

Morse

1.00)< 10-3a) 2.00)<10-4a)

4.06)<10-3b) 7.10X 10-4b) 1.60)<10 3c) 5.00)< 10-4 c) 1.20)< 10-3 c) 3.00X 10-4c)

2.13)<10-3 1.36)<10-3 0.876a) 6.40X10-4 6.09)<10-4 0.156al 1.06)<10-3 2.87)<10-3 0.629d) 5.04× 10-4 8.47)<10-4 0.329a) 7.80× 10 -3 8.62× l0 -3 0.280e) 5.24)<10-4 6.40×10 -4 1.000f) 0.093 f)

1.00)< 10-s e)

copper (00! (111 a) From ref. [11]. 450

present results

b) From ref. [31].

c) From ref. [30].

S-C

d) From ref. [12].

exp.

theor.

Morse

S-C

0.902b) 0.269b) 0.630c~ 0.056c) 0.380~) 0.026c)

0.799 0.089 0.298 0.016 0.249 0.064

0.832 0.106 0.530 0.086 0.243 0.059

e) From ref. [46].

r) From ref. [47].

Volume 180, number 6

PHYSICS LETTERSA

4. Conclusion

As already mentioned, surface diffusion plays an i m p o r t a n t role in u n d e r s t a n d i n g surface-related phenomena. The goal in this study is to know how an adsorbed atom migrates on a rough or smooth crystal plane at the atomic scale. A classical calculation has been performed for the dynamic motion of an adatom on a metal surface using a two-body Morse potential or a many-body S u t t o n - C h e n potential. The results of the Einstein equation and Arrhenius expression show interesting dynamical features: some of the initial kinetic energies are spent in surmounting surface energy barriers a n d their relaxed energies become translational energies which can be assigned to be the corresponding temperatures. As expected, the temperature dependence of the diffusion constants do exhibit Arrhenius behavior for the adatom. The estimated exponential prefactors agree reasonably with the experimental and theoretical values. More significantly, typical values of the random walk exponential prefactors of fcc(001) and fcc(1 1 l ) surfaces are generally predicted in the neighborhood of 10 -3 and l 0 -4 cm2/s, respectively. Generally speaking, the overall mobility of adatoms on the (111 ) surface is much higher than the mobility of adatoms on the rougher (001) surface due to activation energies increasing with the atomic roughness of the surface. This finding will allow us to conclude more clearly that the activation energies are generally larger for j u m p s between sites on atomically rough surfaces than for smooth surfaces. Via the application of proper model potentials, the M D simulations of an adatom on the perfect Rh, Ni and Cu surfaces might provide a realistic picture of the adatom motion. We feel that the level of agreem e n t indicated in table 1 for six different systems suggests that comparison of calculated molecular dynamics with previously measured Arrhenius parameters is meaningful. However, a possible application of the angular dependence of three-body potential or a detailed e x a m i n a t i o n of the energetics of surface diffusion may be needed to verify the present results.

Acknowledgement

This research has been financially supported by the

20 September 1993

Institute of Physics of the Academia Sinica in Taiwan u n d e r the project Surface Physics.

References

[ 1] S.J. Lombardo and A.T. Bell, Surf. Sci. Rep. 13 ( 1991 ) 1. [2 ] T. Ala-Nissilaand S.C. Ying, Prog. Surf. Sci. 39 (1992) 227. [ 3 ] G. WahnstrSm, Rate equations, rate constants, and surface diffusion, in: Interaction of atoms and molecules with solid surface, eds. V. Bortolani, N.H. March and M.P. Tosi (Plenum, New York, 1990). [4 ] G. Ehrlichand K. Stolt, Annu. Rev. Phys. Chem. 31 (1980) 603. [5]E.W. Miiller and T.T. Tsong, Field ion microscopy, principles and applications (Elsevier, Amsterdam, 1969). [6] T.T. Tsong, Phys. Today (May 1993) p. 24. [7] G. Ehrlich, Phys. Today (June 1981 ) p. 44. [8] T.T. Tsong, Prog. Surf. Sci. l0 (1980) 165. [9] D.W. Bassett and P.R. Webber,J. Phys. C 9 (1976) 2491. [ 10] D.W. Bassett and P.R. Webber, Surf. Sci. 70 ( 1978) 520. [ ! 1] G. Ayrault and G. Ehrlich, J. Chem. Phys. 60 (1974) 281. [ 12] R.T. Tung and W.R. Graham, Surf. Sci. 97 (1980) 73. [ 13] G. Ehrlichand F.G. Hudda, J. Chem. Phys. 44 (1966) 1039. [ 14] P.L. Cowan and T.T. Tsong, Phys. Lett. A 53 ( 1975) 383. [ 15] T.T. Tsong, P. Cowan and G.L. Kellogg,Thin Solid Films 25 (1975) 97. [ 16] W.R. Graham and G. Ehrlich, Thin Solid Films 25 ( 1975) 85. [ 17 ] W.R. Graham and G. Ehrlich, Surf. Sci. 45 (1974) 530. [ 18] D.W. Bassett and M.J. Parsley,J. Phys. D 2 (1969) 13. [ 19] T.T. Tsong and G.L. Kellogg,Phys. Rev. B 12 (1975) 1343. [20] D.A. Reed and G. Ehrlich, Philos. Mag. 32 (1975 ) 1095. [21 ] S.C. Wangand G. Ehrlich,Phys. Rev. Lett. 62 (1989) 2297. [22] S.C. Wangand G. Ehrlich, Surf. Sci. 224 (1989) L997. [23] T.T. Tsong and C.L. Chen, Phys. Rev. B 43 ( 1991 ) 2007. [24 ] C.L. Chertand T.T. Tsong,Phys. Rev. Len. 64 (1990) 3147. [25] C.L. Chen and T.T. Tsong, Surf. Sci. 246 ( 1991 ) 13. [26] C.L. Chen and T.T. Tsong, Phys. Rev. B 41 (1990) 12403. [27] C.L. Chen and T.T. Tsong,J. Vac. Sci. Technol.A 10 (1992) 2178. [28] G. Ehrlich, Surf. Sci. 246 ( 1991 ) 1. [29] J.C. Tully, G.H. Gilmer and M. Shugard, J. Chem. Phys. 71 (1979) 1630. [30] C.L. Liu, J.M. Cohen, J.B. Adams and A.F. Voter, Surf. Sci. 253 (1991) 334. [31 ] H.K. McDowell and J.D. Doll, J. Chem. Phys. 78 (1983) 3219. [32] H.K. McDowell and J.D. Doll, Surf. Sci. 121 (1982) L537. [33] J.D. Doll and H.K. McDowell, J. Chem. Phys. 77 (1982) 479. [34] J.D. Doll and H.K. McDowell, Surf. Sci. 123 (1982) 99. [35] S.C. Park and J.M. Bowman, J. Chem. Phys. 80 (1984) 2191. [36] P.J. Feibelman, Phys. Rev. Lett. 65 (1990) 729. 451

Volume 180, number 6

PHYSICS LETTERS A

[ 37 ] C.L Kellogg and P.J. Feibelman, Phys. Rev. Lett. 64 (1990) 3143. [38] P.G. Flahive and W.R. Graham, Surf. Sci. 91 (1980) 449. [39] M.S. Daw and M.I. Baskes, Phys. Rev. Lett. 50 (1983) 1285. [40] M.S. Daw and M.I. Baskes, Phys. Rev. B 29 (1984) 6443. [41 ] M.W. Finnis and J.E, Sinclair, Philos. Mag. A 50 (1984) 45. [42] A.P. Sutton and J. Chen, Philos. Mag. Lett. 61 (1990) 139. [43] M.P. Allen and D.J. Tildesley, Computer simulation of liquids (Clarendon, Oxford, 1987).

452

20 September 1993

[ 44 ] D. Kahaner, C. Moler and S. Nash, Numerical methods and software (Prentice-Hall, Englewood Cliffs, 1989). [45 ] J.H. Noggle, Physical chemistry, 2nd Ed. (Scott-Foresman, Glenview, 1989) p. 37. [46] H.J. Ernst, F. Fabre and J. Lapujoulade, Phys. Rev. B 46 (1992) 1929. [ 47 ] P. Wynblatt and N.A. Gjostein, Surf. Sci. 12 (1968) 109. [48] A.V. Voter and J.D. Doll, J. Chem. Phys. 80 (1984) 5832.