The onset of disorder in Al(110) surfaces below the melting point

Surface Science 220 (1989) L693-L700 North-Holland, Amsterdam

L693

SURFACE SCIENCE LETTERS THE ONSET OF DISORDER IN AI(II0) SURFACES BELOW THE MELTING POINT P. STOLTZE, J.K. NI~RSKOV Laboratory of Applied Physics, Technical University of Denmark, DK-2800 Lyngby, Denmark

and Uzi L A N D M A N School of Physics, Georgia Institute of Technology, Atlanta, GA 30332, USA Received 26 May 1989; accepted for publication 28 June 1989

Molecular dynamics simulations of the disordering of AI(ll0) surfaces below the melting point are presented. Effective medium theory has been used to calculate the inter-atomic many-body interactions of the metallic system. The origin of the anisotropy of the disordering with respect to the direction along the surface deduced from LEED experiments is suggested to be mainly a consequence of the g2 dependence of the Debye-Waller factor even for non-harmonic interactions. Further evidence is presented showing that the disordering of the first layers coincides with the onset of surface vacancy formation.

As the temperature of a solid surface is increased the melting has been shown in a number of cases to be preceded by a disordering, which is usually ascribed to as pre-melting. Ion scattering studies of Pb surfaces [1] have shown a rapid decrease in the short range order starting around 150 K below the melting point. First a single and then a few layers disorder and as the temperature is increased toward the melting point, the thickness of the disordered or pre-melted layer diverges. LEED experiments on Pb and A1(110) surfaces [2,3] have indicated the disappearance of long range order well below the bulk melting point. Furthermore, X-ray scattering studies on Pb(110) and (111) indicate disordering even at half the melting temperature [4]. Low energy He scattering experiments on Pb(ll0) also indicate liquid-like behavior at the surface below the bulk melting point [5]. Similar observations have been made on non-metal surfaces of, for instance, the rare gases [6]. The premelting behavior is strongly dependent on the crystal orientation for metals. On Pb the open (110) surface shows a very pronounced effect whereas the close packed (111) surface shows virtually no effect [1]. 0039-6028/89/$03.50 © Elsevier Science Publishers B.V. (North-Holland Physics Publishing Division)

L694

P. Stoltze et al. / Disordering of Al( l l O) below the melting point

Most recently Prince, Breuer and Bonzel [3] have suggested that for a given surface the disordering is anisotropic. They showed that the decrease in the measured LEED intensity with temperature on a P b ( l l 0 ) surface is faster in the (10) direction parallel to the close packed rows than in the perpendicular (01) direction after the harmonic Debye-Waller part of the intensity decrease has been divided out. Simulations of surface pre-melting have mostly been performed using simple pair potentials most suited for rare gas solids [7-9]. We have performed molecular dynamics simulations for A1(110) and (111) surfaces in which the effective medium theory is used to calculate the interactions [10]. Thus, unlike previous simulations of this phenomena, the potentials used in this work are many-body in nature, derived from density functional theory within the context of effective medium theory [11]. The simulations exhibit signatures of the pre-melting phenomena observed experimentally, i.e., the structure constant for the outermost layers falls off more rapidly with temperature than the deeper lying ones indicating the gradual loss of order starting in the first layer and moving inward with increasing temperature. The crystal face specificity of the phenomenon was also observed in the simulations [10]. Only a limited number of layers were treated in the simulations describing the disordering of the first few surface layers well below the melting point. The asymptotic behavior as the temperature gets very close to the melting point is outside the reach of this study. Theoretical descriptions of surface pre-melting range from Landau theory [12] to thermodynamical models [1]. These approaches have the common characteristic that they are continuum theories. Van der Veen and coworkers [1] are able to give a very detailed account of the growth of the liquid like layer on top of the surface close to the melting point and of the crystal face specificity of this behavior based on measured or calculated surface tensions. What is still lacking is an understanding of the initiation of the disordering in the first few surface layers. The description of the atomic structure of the defects created in this process is outside the scope of the continuum theories. In the present paper we present further results based on molecular dynamics simulations of A1(110). The aim is threefold: (i) to introduce a slightly improved simulation where more neighbors have been included in the force calculations; (ii) to show that while analysis of our data in the manner suggested by Prince et al. [3] results in directional anisotropy of the pre-melting, other ways of analyzing the same data indicate the absence of such an anisotropy; and (iii) to present further evidence that the initial disordering is associated with the formation of surface vacancies. The simulations are performed as described in ref. [10] with one main exception. We have improved the force calculation by including more than just the first neighbors in the density sums as described in detail in appendix D of ref. [11]. In our current simulations we have included all atoms within a sphere

P. Stoltze et al. / Disordering of Al( l l O) below the melting point

L695

of radius R = 8.54 Bohr radii. Increasing the radius of the sphere has no detectable effect on the results. While limiting the contributions in the density sums to the first nearest neighbor shell, as was the case in our previous simulations [10], does not affect significantly the properties of the crystalline material; properties of the molten state, and in particular the melting temperature, calculated within this scheme are sensitive to the position and steepness, of the spline used to cut off the densities. Therefore we recommend the current procedure (i.e. inclusions of contributions from all atoms within a sphere of radius R = 8.54 Bohr radii), which alleviates the above mentioned problems. The simulations have been performed in a computational box consisting of 512 atoms with periodic boundary conditions in the two directions parallel to the surface. The movements of the atoms in the 12 upper layers are determined via numerical integrations, using the Verlet algorithm with a time step of 6 fs, of the classical equations of motion. The 4 layers under the dynamical substrate are kept fixed at their perfect lattice positions. The size of the computational box is determined at each temperature from an independent calculation of the bulk lattice constant. This arrangement is well suited for studying the disordering of the outmost layers, but becomes less well suited close to the bulk melting temperature. After thermalization over 50 ps, we typically sampled averages over 5000 timesteps. The use of a static substrate eliminates the problems associated with ensembles with 2 free surfaces. However, the lower dynamic layers will feel the presence of the static substrate as an ordering field. Close to the bulk melting point the number of melted layers becomes comparable to the total number of dynamic layers and the presence of the static substrate becomes a problem. By comparing to simulations using a large number of dynamical layers and to simulations well above the melting point, we estimate that this ordering field 1.0

0.8 0.6

0.4 0.2

0.0

200

400 600 BOO Temperature (K)

1000

Fig. 1. The absolutc squarc of the calculated " L E E D structure factor" eq. ( I ) for an A l ( l l 0 ) surface s h o w n for various reciprocal lattice vectors g = (2, 0), (0, 2), and (1, 1).

L696

P. Stohze et aL / Disordering of Al( l l O) below the melting point

does not significantly influence any of the results presented here. Doubling the dimensions of the computational box parallel to the surface or propagating the simulation further in time did not produce any change in the observed behavior. The results of the present simulations do not differ qualitatively from those published in ref. [10]. At the quantitative level the main difference is that the melting temperature has increased by approximately 200 K. This is illustrated in fig. 1, which shows the absolute square of the structure factor ( ~ exp(i¢- g) e x p ( -

zJX))

Steeo(g) =

(1)

Here the sum is over all the atoms in the computational box, z is the distance perpendicular to the surface, and the average is over configurations collected after thermalization of the system. When the sampling depth k is chosen to be of the order 1.5 ,~ ]SLEEO(g)] 2 has strong similarities to the intensity measured in a LEED experiment (although multiple scattering effects of major importance in the LEED process are not included in eq. (1); we assume that this does not affect the temperature dependence significantly). Comparison with the similar curve in ref. [10] shows the qualitative behavior to be the same, whereas the current result is shifted up in temperature by approximately 200 K relative to the old one. The comparison to experiment is equally good for both simulations, one being 100 K below the experimental results of 1.0

0.8 -

=

• , •

0.6

0.4

0.2

o.o

' 2~o ' 46o ' 660 ' ado ' lobo Temperature

(K)

Pig. 2. The absolute square of the " L E E D structure factors" f r o m fig. I corrected for the D e b y e - W a l l e r part of the t e m p e r a t u r e dependence. The D e b y e - W a l l e r factor is fitted to the low t e m p e r a t u r e (less t h a n 400 K) p a r t of the curves in fig. 1. The s t r u c t u r e factors clearly show an a n i s o t r o p y as suggested b y Prince et al. [3].

P. Stoltze et a L / Disordering of Al(l lO) below the rneltingpoint

L697

Blankenhagen et al. [2] and one being 100 K above. We note that the potential calculation used here does not involve fitting of the parameters and we consider an agreement between a calculated and experimental melting temperature to within 10% to be surprisingly good. In fig. I we have included the results for several reciprocal lattice vectors g. The g---(2, 0) [13] structure factor is clearly seen to decrease considerably faster than the others. If we fit the low temperature part of the curves to a Debye-Waller form e x p ( - 2 W ( g ) T ) and divide this part out of ] SLEED(g) I 2 as suggested by Prince et al. [3] we get the result shown in fig. 2. In accordance with the experimental results of Prince et al. [3] the Debye-Waller corrected g = (2, 0) intensity falls off before and faster than the others indicating an anisotropy in the disordering process. One can, however, make another analysis of the results in fig. 1. The Debye-Waller constant W(g) is proportional to the absolute square g2 of the momentum transfer g in the harmonic approximation. A more detailed analysis shows [14] that this still holds to a good approximation even for an anharmonic crystal, except that in this case W(g) becomes temperature dependent. It is therefore tempting to plot I S " [ 2 given by

(2)

In I S " ( T ) I = In I SLEEp(g) 12/g 2,

as a function of temperature in order to eliminate the g dependence of the structure factors. This is done in fig. 3. The curves for all the different momentum transfers now basically fall on top of each other. On the basis of fig. 3 it cannot be concluded that there is an anisotropy in the disordering in the present simulation for AI(ll0). Rather the explanation of the anisotropy observed in figs. 1 and 2 appears to be that the g dependence of the structure

1.0

0.B

~__ 0.6 0.4

0.2

o

26o

460

660

Temperature

(K)

860

,o'oo

Fig. 3. The absolute square of the "structure factor" S"(T) from eq. (2) for the same three reciprocal lattice vectors as in figs. 1 and 2. Now the anisotropy has basically disappeared.

P. Stoltze et al. / Disordering of Al( l l O) below the melting point

L698

2

1.0

10 -1

0 . 5 :> 0 C~

10 -~ o/

10 -s

2()0

4 0, 0 : ,

~6-0' -0, - ~

, BOO

1 0 0 00.0

Temperature (K) Fig. 4. The absolute square of the g = (0, 2) layer structure factor eq. (3) for the first two layers as a function of temperature. Also shown (dashed) is the average occupation of the layer on top of the surface. The decrease in structure factor for the first layer is seen to be closely correlated to the number of atoms that have moved out of the crystal and onto the first layer.

factors amplifies the fall off of the g = (2, 0) i n t e n s i t y d u e to a significantly larger length of the m o m e n t u m transfer vector in this case. W e now turn to the question of the origin of the p r e - m e l t i n g behavior. In fig. 4 we show the a b s o l u t e square of the l a y e r structure c o n s t a n t S~ for the two o u t e r m o s t layers of the crystal,

s,(g)

=

,

(3)

where z~ defines the m i d p o i n t b e t w e e n layer l a n d l + 1 as d e t e r m i n e d b y a b u l k calculation a n d O(z) is the step-function. T h e g r a d u a l d i s o r d e r i n g is clearly seen. A n a l y z i n g the d a t a further we find that the d e c r e a s e in the structure factor for the o u t e r m o s t layer is c o r r e l a t e d with the o n s e t of the m i g r a t i o n of a t o m s f r o m the o u t e r m o s t layer to p o s i t i o n s on the surface. In fig. 4 we also show the n u m b e r of a t o m s o b s e r v e d on t o p of the surface. This is the n u m b e r of vacancies c r e a t e d in the first a n d (to a lesser e x t e n d ) lower layers. T h e decrease of the o r d e r in the first layer is seen to be d i r e c t l y r e l a t e d to the onset of defect f o r m a t i o n . This o b s e r v a t i o n is further s u p p o r t e d b y the s n a p s h o t s shown in fig. 5. T h e figure shows the i n s t a n t a n e o u s p o s i t i o n s of the a t o m s of the two o u t m o s t layers at two t e m p e r a t u r e s , one j u s t b e l o w the t e m p e r a t u r e w h e r e the first layer d i s o r d e r s a n d one j u s t above. T h e layer n u m b e r has b e e n i n d i c a t e d on the i n d i v i d u a l atoms. W h e r e a s the layer n u m b e r b e c o m e s of l i m i t e d signifi-

P. Stoltze et al. / Disordering of Al(l lO) below the melting point

L699

1 2 12 1

Fig. 5. Snapshots of the configuration of atoms in the first layers at two temperatures above and below the disordering temperature of the first layer. The atoms are labeled by their layer, 1 and 2 for the first and second layer, and 0 for the layer on top of the first one. cance for s t r o n g l y d i s o r d e r e d systems where there are m a n y a t o m s in the region a r o u n d the p l a n e s d i v i d i n g the layers, it gives a r o u g h i m p r e s s i o n of the height of the a t o m even in these case. U p to the t e m p e r a t u r e w h e r e [$1[ 2 begins to d r o p the (110) surface structure with the close p a c k e d rows in the (10) d i r e c t i o n is intact. T h e p i c t u r e is q u a l i t a t i v e l y different a b o v e the diso r d e r i n g t e m p e r a t u r e for the first layer. N o w the l a y e r a b o v e the surface is significantly p o p u l a t e d a n d the u n d e r l y i n g layer reacts s t r o n g l y to the a t o m s a b o v e the surface b y d i s t o r t i n g locally a r o u n d them. It is this d i s t o r t i o n that i n d u c e s the decrease of the structure factor. T h e loss o f a t o m s in the o u t e r m o s t

L700

P. Stoltze et al. / Disordering of Al(l l O) below the melting point

l a y e r d o e s n o t b y itself give rise to a d e c r e a s e in i n t e n s i t y ( g i v e n t h e n o r m a l i z a t i o n in eq. (3)) if the r e m a i n i n g a t o m s w o u l d stay at t h e l a t t i c e positions. A s p o i n t e d o u t in o u r p r e v i o u s p a p e r [10] a d i s o r d e r i n g p r o c e s s t h a t i n v o l v e s the f o r m a t i o n o f s u r f a c e d e f e c t s like v a c a n c i e s r e a d i l y e x p l a i n s w h y t h e d i s o r d e r i n g o f the c l o s e p a c k e d s u r f a c e s b e l o w t h e m e l t i n g p o i n t is m u c h w e a k e r t h a n o n t h e o p e n surfaces. T h e e n e r g y i n v o l v e d in t a k i n g a n a t o m o u t o f the s u r f a c e l a y e r a n d p l a c i n g it a b o v e it is m u c h h i g h e r o n the c l o s e p a c k e d t h a n o n t h e o p e n s u r f a c e s (1.3 eV o n A I ( l l l ) as c o m p a r e d to 0.3 e V o n A1(110)). T h i s w o r k was s u p p o r t e d b y t h e D a n i s h R e s e a r c h C o u n c i l u n d e r F T U c o n t r a c t n u m b e r 5.17.1.1.11, E E C c o n t r a c t n u m b e r S T 2 J - 0 2 8 6 - C ( T T ) (P.S., J . K . N . ) a n d the U S D e p a r t m e n t o f E n e r g y g r a n t n u m b e r F G 0 5 - 8 6 E R 4 5 2 3 4 (U.L.).

References [1] J.W.M. Frenken, P.M.J. Maree and J.F. van der Veen, Phys. Rev. Letters 54 (1985) 134; B. Pluis, A.W. Denier van der Gon, J.W. Frenken and J.F. van der Veen, Phys. Rev. Letters 59 (1987) 2678; for a review see: J. van der Veen, B. Pluis and A.W. Denier van der Gon, in: Chemistry and Physics of Solid Surfaces, Vol. VII (Springer, Berlin, 1988). [2] P. von Blankenhagen, W. Schommers and V. Voegele, J. Vacuum Sci. Technol. A 5 (1987) 649. [3] K. Prince, U. Breuer and H.P. Bonzel, Phys. Rev. Letters 60 (1988) 1146. [4] P.H. Fuoss, L.J. Norton and S. Brennan, Phys. Rev. Letters 60 (1988) 2046. [5] J.W.M. Frenken, J.P. Toennies and Ch. W~511,Phys. Rev. Letters 60 (1988) 1727. [6] D.-M. Zhu and J.G. Dash, Phys. Rev. Letters 57 (1986) 2959; 60 (1988) 432. [7] J.Q. Broughton and L.V. Woodcock, J. Chem. Phys. 79 (1983) 5119. [8] S. Valkealahti and R.M. Nieminen, Phys. Scripta 36 (1987) 646. [9] V. Rosato, G. Ciccotti and V. Pontikis, Phys. Rev. B 38 (1986) 1860. [10] P. Stoltze, J.K. N~rskov and U. Landman, Phys. Rev. Letters 61 (1988) 440. [11] K.W. Jacobsen, J.K. N~rskov and M.J. Puska, Phys. Rev. B 35 (1987) 7423. [12] R. Lipowsky and W. Speth, Phys. Rev. B 28 (1983) 3983; R. Lipowsky, U. Breuer, K.C. Prince and H.P. Bonzel, Phys. Rev. Letters 62 (1989) 913. [13] The recipocal lattice vectors are given in units for 2~r/a where a is the lattice constant. The (10) reflection is bulk forbidden. [14] A.A. Maradudin and P.A. Flinn, Phys. Rev. B 129 (1963) 2529.