PHYSICAL REVIEW B VOLUME 50, NUMBER 16 15 OCTOBER 1994-II

Surface anharmonicities and disordering on Ni(1QQ) and Ni(11Q)

Yvon Beaudet and Laurent J. Lewis*Departement de Physique et Groupe de Recherche en Physique et Technologie des Couches Minces, Universite de Montreal,

Case Postule 6188, Succursole Centre Vi-lle, Montreal, Quebec, Canada HSC 9J7

Mats PerssontDeportment of Applied Physics, Chalmers University of Technology and University of Goteborg, S-$12 96, Goteborg, Sweden

(Received 15 March 1994)

Recent observations on various low-index metallic surfaces of an extraordinary thermal atten-uation of diffracted peak intensities suggest that anharmonic efFects are important at high tem-peratures, though such observations have also been interpreted in terms of disordering processes(premelting and roughening). In order to address this problem, we have carried out an extensivemolecular-dynamics investigation of the evolution in temperature of the (100) and (110) surfacesof Ni, with the interactions between atoms described using embedded-atom method potentials. Weobserve both surfaces to sufFer an anomalously large thermal expansion, with a concomitant rapidincrease of the mean-square amplitudes of vibration; the latter are found in both cases to be largerin the plane of the surface than out of the plane. An analysis of our simulated low-energy electron-diKraction intensities supports the experimental interpretation of the measured ones in terms ofenhanced mean-square displacements. Indeed, anharmonicities are found to appear at temperatureswell below the onset of disordering, which occurs at 1200 and 1000 K for the (100) and (110) sur-faces, respectively (via in-plane difFusion and adatom-vacancy formation). The onset of disorderingproceeds difFerently on the two surfaces: while it involves essentially only the outermost layer on the(100) surface, both the first and second layer participate in the process on the (110) surface. Wealso calculate the phonon spectra at high-symmetry points of the surface Brillouin zone, which wefind to be in good agreement with experimental values at room temperature. We predict, further,that the temperature-dependent frequency shifts are not significantly softened in comparison withbulk phonons, which suggests that surface phonons are afFected by anharmonicity in a nontrivialway. Several of the observed properties can be understood in terms of the larger displacements thatsurface atoms can a8'ord, because of their reduced coordination compared to bulk atoms.

I. INTRODUCTION

The thermal behavior of crystalline surfaces, and, inparticular, metallic, has raised a lot of attention in re-cent years due to a large extent to re6nements in bothexperimental and theoretical techniques. Low-energyelectron difFraction (LEED), i s elastic He scattering,medium-energy ion scattering, ' and electron energy lossspectroscopy are all techniques that are well suited tothe study of the evolution in temperature of the structure(both static and dynamic) of metallic surfaces. These ex-periments indicate that such quantities as surface ther-mal expansion, mean-square displacernents, and phononline shapes often exhibit behaviors which are "anoma-lous" when compared to those in the bulk. It is our pur-pose to assess, using numerical simulation methods, theimportance of anharmonic effects and disordering pro-cesses on the evolution in temperature of metallic sur-faces.

One direct consequence of anharmonic terms in thepotential seen by atoms is thermal expansion which, inview' of the above, is expected to be much larger in thedirection perpendicular to a surface than in the bulk. In-

deed, recently, LEED experiments have been reportedrevealing that the (100) surface of Ni sufFers a large out-ward relaxation at high temperature, the thermal expan-sion reaching a value considerably larger —by a factor ofabout 20—than in the bulk. We know of only one otherreported case of an anomalously high thermal expansioncoeKcient, namely, Pb(110),s though indirect evidenceof increased anharmonicity is provided by the thermalattenuation of diffracted-beam intensities, which can beunderstood in terms of Debye-tA'aller factors with anoma-lously large mean-square amplitudes of vibration. '

Theoretical understanding of anharmonicity is a long-standing problem, especially in cases where the crys-talline symmetry is broken. In the best of cases, onehas to resort to some sort of quasiharmonic approxima-tion, where the formulation of the problem is harmonicbut quantities such as mean-square displacements are al-lowed to deviate from their harmonic behavior; other-wise, techniques based on many-body theory must beused, where anharmonic terms are at best treated pertu-batively. The numerical method of molecular dynam-ics (MD), however, is not embarassed by these diffi-culties, as it provides an exact solution to the classi-

0163-1829/94/50(16)/12084(20)/$06. 00 12 084 1994 The American Physical Society

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 085

cal many-body problem, i.e., fely includ ng al aWm-monic contributions. The main ingredient of such sim-ulations is the potential of interaction between atoxns.The advent in recent years of realistic, xnaterial-speci6cpotentials for metallic systems, such as the semiempiri-cal embedded-atom method (EAM) (Refs. 11—13) or theefFective-medium theory (EMT) (Ref. 14) (as opposed togeneric, like Lennard-Jones), have opened the way to adetailed understanding of numerous atomic-scale prop-erties (and their temperature dependence), in particular,at surfaces where the broken symmetry is incorporatedat the outset into the many-body formulation of the po-tential.

The MD-EAM (or EMT) tandem, therefore, offers analmost ideal solution to the problem of investigating thedynamics of metal surfaces (modulo finite-size limitationsand quantum efFects at low temperatures); these, in fact,have been the object of several such studies recently, deal-ing with premelting and melting on Ni(110) (Ref. 15) and

(100),is Cu(110) (Ref. 17) and Al(110),is anharmonic-ity, vibrations, and disorder on Cu(ill), is (100),is'2s

(110),is 2 2 Ag(100),2o and (110), to name a few.On the basis of MD studies such as the above and

others, it is becoming more and more evident that theobserved anoxnalies at surfaces are mostly a consequenceof enhanced anharmonicity rather than disordering—either in-plane (premelting) or out-of-plane (roughen-ing), resulting from the loss of symmetry compared tothe bulk, which allows the atoms to probe a larger (andpresumably more anharmonic) portion of the potential.Most calculations have focused on the open (110) sur-faces of face-centered-cubic metals; it is an interestingproblem, therefore, to determine how ubiquitous the an-harmonic efFects are, and what the nature is of the dis-ordering processes that take place at high temperatures,by contrasting the (110) surface of a specific fcc metalwith a more close-packed surface like (100).

We present here an account of our MD-EAM investi-gation of anharmonic efFects and disordering at the (100)and (110) surfaces of Ni, and propose a detailed compar-ison with available experimental data, in particular, theLEED measurexnents of Cao and Conrad. ' Parts of ourwork on the Ni(100) surface have been presented else-where and are not reproduced in detail here. ' Our re-sults are generally in good agreement with the availableexperimental data, and unambiguously reveal the roleplayed by anharmonicity on these surfaces: we 6nd in-deed that the extraordinary attenuation of the LEED in-tensities takes place well before any form of disorder setsin, and observe a concomitant "anomalous" growth of theexpansion coeKcients. Our calculations further show themean-square amplitudes of vibration to exhibit stronganharmonic behavior, in addition to being much largerthan their bulk counterparts; this behavior is a conse-quence of the reduced coordination of surface atoxns, re-sulting in enhanced motion of the atoms away from equi-librium@, which allows them to probe anharmonic regionsof the potential unaccessible to bulk atoms. The ampli-tudes of vibration, further, are found for both systems tobe larger in the plane of the surface than out of the plane,a somewhat counter-intuitive result. The anharxnonic ef-

fects, however, do not lead to pronounced softening ofthe surface phonons (compared to bulk ones) at elevatedtemperatures, which suggests that collective excitationsat the surface are afFected by anharmonicity in a nontriv-ial way.

We present also a detailed study of the high-temperature surface disordering regime; we observe dis-ordering to proceed via in-plane difFusion and adatom-vacancy formation, though in manifestly difFerent wayson the two surfaces: while the onset of disordering in-volves essentially only the outermost layer on the (100)surface, both the 6rst and second layer participate in theprocess on the (110) surface. Thus, on (100), we observein-plane disordering to take place before out-of-plane dis-ordering, by a mechanism involving the concerted xno-tion of several atoxns and the transient occupation of anadatom site; on the (110) surface, disorder sets in simul-taneously in plane and out of plane.

We begin in Sec. II by giving full details of our model.In Sec. III, we present and discuss our results, start-ing with a study (Sec. IIIA) of the evolution in tem-perature of the slabs that we use to describe the sur-faces; we characterize the latter in terms of a layerdensity function, adatom-vacancy formation, van Hoveself-correlation function, and particle trajectories. Next(Sec. IIIB), we proceed with a discussion of surface ther-mal expansion coefBcients which we find indeed to beanoxnalously large. We present in Sec. III C the results ofour calculations of the mean-square amplitudes of vibra-tion. In Sec. III D, we calculate the phonon f'requency dis-tributions at selected high-symmetry points of the two-dixnensional surface Brillouin zones, as well as their varia-tions with temperature; we find good agreexnent with theexperimental values measured at room temperature, andour calculations reveal that the softening of the phononswith texnperature is not anomalously large in compari-son to the corresponding phenomenon in the bulk. InSec. IIIE1, we propose an explicit comparison of ourdata with the LEED measurements of Cao and Conrad,which supports their interpretation of enhanced surfaceanharmonicity in terms of increased vibrational axnpli-tudes. Layer structure factors are discussed in detail inSec. IIIE 2, both as a function of temperature and layernumber, i.e., distance away from the surface. A sumxnaryand conclusions are given in Sec. IV.

II. MODEL AND COMPUTATIONAL DETAILS

The semiempirical EAM represents a good compro-mise between efFective potentials and 6rst-principlesapproaches, i.e., between accuracy of the model andfeasability. In particular, EAM includes —albeit in anapproximate way —the many-body nature of interac-tions in metals, essential to a proper description of theirenergetics. Indeed, as we will see in the course of our dis-cussion, several of our observations can be traced backto this simple fact. Two-body potentials, in contrast,despite their widespread use for metallic systexns, fail toprovide a proper description of quantities such as vacancyforxnation energies, elastic constants and, most impor-

12 086 YVON BEAUDET, LAURENT J. LE%IS, AND MATS PERSSON 50

(c)

lEa1

XJ E

b,

()~ II ~

a,

4L1

Xka

b)

0 ~ ~FIG. l. (a) Geometry of the (100) surface; the basis vec-

tors for the two-diinensional direct lattice are ai ——(a/~2)xand az = (a/v 2)y, with a the lattice parameter (the unitvectors are defined in Ref. 26); full and open circles refer tothe first and second layer, respectively. (b) Geometry of the(110) surface; here, ai = (a/v 2)x and a2 ——ay. (c) and (d)Corresponding Srst Brillouin zones in reciprocal space; indi-cated are the zone boundary points M and 4 for the (100)surface, and X and Y for the (110) surface.

tant for our purpose, surface energies, found to be muchhigher than the experimentally-observed ones. In fact,the use of EAM for surfaces has been extremely success-ful in describing energetics, relaxation, reconstruction,segregation, phonon dispersion, thermal behavior, pre-melting, and melting.

Our model for the (100) surface consisted of a slab of12 layers, each containing 64 atoms arranged in an fccpattern, oriented perpendicular to the z direction. Thereal-space geometry of the surface is shown in Fig. 1(a),where we also indicate the lattice vectors used to generateeach layer. 2s Likewise, in Fig. 1(b), we show the structureand the lattice vectors used in our model for the (110)surface. Again, in this case, a slab made up of 12 layerswas used, but now with 7x10 unit cells, i.e., 70 atoms,in each layer. In order to assess the importance of finite-size eH'ects, we have carried out preliminary calculationson smaller systems, and found no significant di6'erencesin the average properties, if only in statistical accuracy.Periodic-boundary conditions were used in the z and ydirections, but the surfaces were entirely &ee in the z di-rection, except that the bottom two layers of each slabwere fixed in their perfect fcc structure in order to mimicthe effect of the presence of the bulk (infinite) crystal;the lattice parameter of this "simulated bulk" was deter-mined by carrying out simulations on a "true bulk, " asdetailed below. We have verified that the presence of thefixed layers have little eQ'ect on the dynamics of the sys-tems near the surface, and in particular contributes neg-ligibly to the total linear momentum. 2 Atoms which lie

near the static layers, however, will perceive them as an"ordering Geld" perturbing their dynamics; this is espe-cially true at elevated temperatures, near melting, wherethese atoms will behave in an unrealistic manner. Notealso that since the static substrate has no kinetic energy,it acts as a thermal insulator, and the mobile atoms willnot loose energy to it.

The simulations of the surfaces were carried out in thecanonical (N, V, T) ensemble using a Nose thermostat tocontrol the temperature. The time step was in all cases(i.e., for both the bulk and the surfaces) set equal to 3.0fs which, we find, leads to an excellent conservation ofenergy. Typical surface runs were 90—125 ps long, fol-lowing a period of equilibration of 135 ps. Several tem-peratures in the range 100—1700 K were investigated (Nimelts, experimentally, at 1726 K, and at about 1740 K inthe universal-equation EAM, Ref. 29). It is important tonote that, since our simulations are purely classical, theresults can strictly be compared with experiment only fortemperatures of the order of, or above, the Debye tem-perature 8Li (390 K for Ni, Ref. 30). Nevertheless, wewill use the results of the low-temperature runs to extractthe harmonic limit of various properties.

The lattice parameter for the simulated bulk was de-termined by performing a series of runs under constant-pressure conditions. More specifically, we used the An-dersen (N, P, T) algorithm —with the external pressureset to zero —which constrains variations of the volumeto be uniform. The system, in this case, was chosento consist of 4x4x4 unit cells, i.e., 256 atoms, whichwe found is suKcient to provide reliable and convergentresults. Each run (16 in total) consisted of a long equili-bration of 243 ps, followed by a "production" run of 60ps, during which statistics were accumulated. %e veri-fied that longer production runs did not acct the resultsin any significant manner.

In addition to yielding information on the tempera-ture dependence of the lattice parameter, the simulationof the bulk system provides a reference against which cer-tain surface properties can be compared, as we will seebelow, but also a test of the validity of the model, sincethe properties of the bulk material are known in moredetail than the corresponding ones at the surfaces. Inparticular, since the original motivation of this work wasto assess the experimental observation of a large ther-mal expansion at the (100) surface of Ni, it appearedimportant to us to verify that, to the least, the thermalexpansion coefFicient of the bulk crystal, o.~, was well de-scribed by the model potential in the temperature rangeconsidered here.

We have tested two versions of the EAM: the early"cubic-spline" version of Daw and Baskes, and the morerecent "universal-equation" formulation of Foiles, Baskes,and Daw. [The potential was truncated at a distance of1.5 a(T), which includes atoms up to the fourth nearest-neighbor shell. ] The calculated lattice constants a{T) forboth cases are displayed in Fig. 2(a), together with theexperimental values "recommended" in Ref. 31, whichare an average over a number of diferent measurements(i.e. , the error bar on the experimental nn is ratherlarge). The universal-equation-EAM values follow closely

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 087

I)

I

3.60—

3.55—

3.50 I x I

(b)

~ ~~ 4.4~, y, g Q 4

i

I i I I i I

400 800 1200 i6007 (K)

FIG. 2. (a) Calculated lattice parameter a as a functionof temperature for the two versions of the EAM discussedin the text: universal equation (solid circles) and cubic spline(solid squares). Also shown are the experimental values (opencircles, Ref. 31). (b) Corresponding thermal expansion coef-Bcients, o.xx = (1/a)da/dT.

III. RESULTS AND DISCUSSION

A. Layer structure and atomic migration

X. Layer numbev densitiee

A useful quantity for characterizing the structure ofthe slabs is n(z), the density along the z direction (i.e.,

the experimental values (modulo a slight overestimationarising for a large part &om fitting the zero-temperatureEAM to the room-temperature value of the lattice pa-rameter) and provide a slightly poorer agreement thanthe cubic-spline-EAM values. However, of utmost im-portance here, the 8lope of the temperature dependenceof the lattice parameter, i.e., essentially the thermal ex-pansion coefBcient, is much better reproduced with theuniversal-equation EAM than with the cubic spline. Thisis seen more explicitly in Fig. 2(b), where we plot now thethermal expansion axe(T) = (1/a)da/dT. Our model, ofcourse, neglects other contributions to the thermal ex-pansion such as those arising &om electrons and vacan-cies, which are present in the measured quantity. In thetemperature range of interest in this work the electroniccontribution is 10% and the vacancy contribution isonly significant close to the melting point. 2 We there-fore used, in all subsequent calculations, the universal-equation formulation of the EAM. It should be notedagain that, at temperatures below the Debye tempera-ture, quantum effects set in and result in deviations be-tween calculated and measured o.~, as can be seen in Fig.2(b) at temperatures below 400 K.

perpendicular to the surface), defined such that n(z) dz isthe number of atoms in a slice of thickness dz at z. Thus,in an ordered —or layered —system, the integral of n(z)over a layer yields the number of atoms in that layer.For a system in equilibrium, n(z) can be obtained byan average over many difFerent configurations (i.e. , overtime). The partitioning of space in terxns of layers willbe particularly useful later in calculating various layer-related properties (such as structure factors).

Sets of plots of n(z) at different temperatures for ourtwo surfaces are given in Figs. 3 and 4. At low temper-atures, n(z) consists of a series of sharp, well-resolved,peaks with the function going to zero between them. Thepartitioning of the z axis in layers, therefore, is unam-biguous, as is the identification, for each atom, of thelayer to which it belongs, since atoms do not migrate &omone layer to another. As temperature increases, the posi-tions of the peaks, still well resolved (though broader andsmaller because of enhanced vibrational motion), moveto larger values. In fact, the evolution in temperature ofthe peaks will be used in the next section to evaluate thesurface thermal expansion coeKcient. At elevated tem-peratures, vibrational amplitudes become so large thatdisorder sets in, with atomic migration taking place be-tween layers, as evidenced by the fact that the minimabetween the peaks close to the surface rise to nonzerovalues; this is reminiscent of the structure of a liquid asobserved in the radial distribution function. The system,as we will see in more detail below, crosses over to a stateof "premelting" upon increasing the temperature. (Re-call that bulk Ni xnelts at 1740 K in this model. ) Thesurface layers, therefore, sufFer more disorder than layersaway Rom the surface, an effect related to premelting,and enhanced by the ordering presence of the static sub-strate.

A most striking result is that disorder takes place muchearlier, and much more efFectively, on the (110) surface incoxnparison to the (100) surface, which is evident whencomparing the n(z) for the two slabs in the region ofthe surface; we will return to this point later. It shouldbe noted that, nevertheless, the concept of layers is welldefined even at the highest temperatures, in particularfor the (100) surface.

2. Adlayer fo~ation

The topmost layer is, consistent with the above ob-servations, the one most a8'ected by disordering efI'ects.Since layers beneath the first are fully occupied (at lowenough temperatures), this leads to the formation of anadlayer on top of the first layer, the adatoms leaving be-hind an equivalent number of vacancies. This is evidentin Figs. 3 and 4 at high temperatures, where an adatompeak develops on top of the surface layer. Naturally,as temperature increases further, the migration processwill reach layers deeper in the substrate. This is clearlydemonstrated in Fig. 5, where we plot, for the two sur-faces, the relative number of vacancies as a function oftemperature, nI v (T), for each layer E up to the fourth;

12 088 YVON BEAUDET, LAURENT J. LE%'IS, AND MATS PERSSON 50

Ni(100)

900

I)

I)

It

1

0.08(a) Ni(100)

. i IL i iliii iiilii Ill I

10001

p iliiiiiilililliiiiliIi1aOO-

g~ 0.06'U

0.04—

0.02

0.002.3—4

Ni(110)' li I

'

Iil QOO-

I, & , I I, Il I Il Ll LJ li ~~ppp-

1

„,ii, i Ii ii Ii Il Il It' iapp-

1

i u Il il 0 il Il il II

g- i

Ii'

I I ~i '14oo-

-„", I(ILlIlkk J1500

1

, iIi, li illII J

, „,I,,

i ~,1

, , 448

z (k)

12 16

FIG. 4. Same as Fig. 3, but for the (110) slab.

, Ll li/lliilliILiIiilig I

1

, Iil iilililIIIIIIA

iii li!i iil Ll Ll Ll IA

16OO-

, I I i I i I Ll I I il il IkI

'

1700

, » Ii Ii LI il It0 8 16 2

z (A)

FIG. 3. Distribution n(z) of atoms through the (100) slabalong the direction perpendicular to the surface at varioustemperatures, as indicated.

0.40 —(

0.30—0.20

0.10

0.00 ' e4

800 1000 1200 1400 1600 1800T (K)

FIG. 5. Relative number of vacancies ng~ in layers 1—4, asindicated, vs temperature, for (a) the (100) and (b) the (110)surface; layer numbered 1 is the topmost layer of the perfectsurface. Also shown (open circles) is the relative number ofatoms in the allayer, n~.

also plotted is the relative number of adatoms in the ad-layer n~(T) equal, by definition, to the total number ofvacancies left behind. Clearly, again, disorder sets in ata much lower temperature on the (110) surface than onthe (100) surface (1000 vs 1300 K, say, but see below),and much more electively on the former than on the lat-ter (notice the difference in scales in Fig. 5). Further,the topmost layer is the sole provider of adatoms up totemperatures of about 1600 and 1250 K on the (100) and(110) surfaces, respectively, above which deeper layersstart to contribute to the adlayer.

Experimentally, in their LEED study of the Ni(110)surface, 2 Cao and Conrad evaluate the number of vacan-cies to be of the order of 4% at 1000 K, and between 5%and 15% at 1200 K. This is in remarkable accord withour numbers of 4% and 13%. Further, Cao and Con-rad also conclude &om their study that roughening takesplace at approximately 1300 K, i.e., approximately thesame temperature at which out-of-plane disorder reachesthe second layer in our calculation. Strictly speaking,the roughening transition is de6ned as that temperatureat which the height-height correlation function diverges(logarithmically); it is not possible to identify preciselysuch a transition in a small-scale MD simulation such asthe present one. However, Hakkinen et a/. have shown

[for the Cu(110) surface], using a combination of lattice-gas Monte Carlo and MD, that the roughening mecha-nism is related to the formation of thermally-generatedadatom-vacancy clusters.

For the (100) surface, now, Cao and Conrad see no out-of-plane disorder till about 1470 K, consistent with ourobservations: the adlayer becomes substantial at about1500 K in our model, as can be seen in Fig. 5(a), itsoccupancy being less than about 0.2 % below that tem-

12 08950

surface-dynamical roor e study of surfaces

b""" 'h d'" """"M

t ereforehtC tesian irections n

)

esl s, an we definne,

perature and th'7 erefore ncontrast

not statlstlcall y szgnlficant. Ipp nce o en-

d l, bot e seen i

'FA

""""""'"(100) d (110

at the adlao e

er foa so i entified

orxnation w

'va ion en-

ia atom-divacancri uted to the

x ( ~ ppo e oa atom

N

G.(r, t) = —) (0 r —r —, , a1G —&

—t;(t) —,(0)])), a.1

where r, (t) is the ath Cartesian

e average i.e.e, an denote s

, t)d = 1. Th)

't etd

hbrian o ps corre

rnperatures

por 1onof thb ofth h'e igh symmetry of thee 100) surface we

e noted that,

.",-h-. a. . t = '-.r, S) = G.(-r, t)

3. pan Ho@ee self-cor relations

mportant insights on the dynamics of the slab

of the

ime t given that th

e wo contributions G'f'h f

""Pic s ructure

rce ot e

a oms awad b fattzce position . t

qui x riumy om their e l b

bshed a ount of hw ed e rst pub-

, r, t in the analysis of

Ni(100)Ii '

I' '

I' ll '

I

1000-- 1400-0.4

0,2

0.00.4 1100-

0.2—

0.00.4

II I I

I

I I ~

I

1200-- 1600-

SURFACE ANH S AND DISpRDER ING PN

0}— 0.2

0.00.4

0.2

I'

I]

'I

1300--I I I

I'

I

I

II\

I'l

1700-

I I

8 101/(ksT) (1/eV)

12

O.o—4 —2 0 2

r (A)

0 2r (A)

FIG. 6. ArrArrhenius plot ofo o o atoms in t

linear fits ases, as indicated. The

di d ' th

FIG. 7. van H. van Hove self-correlatcaro f on

tedt

ed, and dashed liures, as indicated. th

spectively., an z directions, re-

12 090 YVON BEAUDET, I.AURENT J. LE%'IS, AND MATS PERSSON 50

only, where r = x, y. In practice, however, statisticalHuctuations arising from the 6nite size of our model pre-vent these identities &om holding precisely.

By definition, in a crystalline system, G, (r, t), plot-ted against r, is characterized by a single peak aboutthe origin, whose width is a measure of the mean ampli-tude of vibration of the atoms around their equilibriumpositions. This is exactly what we observe in our two sys-tems at "low" temperatures (e.g. , Fig. 7 at 1000 K; wewill return to a detailed discussion of these amplitudes inSec. IIIC). At elevated temperatures, however, in con-formity with our observations in n(z), features appear in

G, (r, t) that are not strictly vibrational in nature, andin fact correspond to diffusional processes. The onset ofsuch di6'usion events takes place, as can be seen fromFigs. 7 and 8, at 1200 K on the (100) surface (thoughit becomes significant only at about 1500 K), and 1000K on the (110) surface. The diffusional motion is quitedifferent on the two surfaces, and we now analyze it indetail.

First, on the (100) surface, for temperatures in therange of 1200—1400 K, we find small but mell-defined

peaks in the x and y directions at positions correspond-ing to the nearest-neighbor distance within the surfacelayer; no secondary peaks is observed however in the zdirection. This result is a little surprising since it sug-gests that diÃusional motion within the layer takes placein absence of vacancies. The occurence of x or y diffu-sion at these temperatures is, however, extremely rare,and a correspondingly minute adatom population has infact been observed in nA —cf. Fig. 5. (Note that ngbeing a "zero-time" property, it can be calculated usingthe full extent of the equilibrium portion of the run; fortime-dependent properties, in order to preserve the sta-tistical quality, correlations are limited to a &action ofthe length of the run; this calculational difference leadsto small discrepancies such as the above. ) In fact, as we

will see in the next section, a detailed analysis of the tra-jectories reveals that this behavior involves the correlatedmotion of surface atoms via the transient occupation ofa hollow adatom site. Our results, nevertheless, are inline with Cao and Conrad's observation of the onset ofin-plane disorder at about 1300 K.

Decimation of the surface layer, by the formation of anadlayer or by diffusion into the second layer, starts (onthis time scale) at approximately 1500 K, as signaled bythe appearance in G, (z, t) of peaks in the positive andnegative z directions, centered around the surface inter-layer distance. At this temperature, further, the intensityof the central peak for lateral motion has dropped downconsiderably in comparison to the corresponding peakfor the perpendicular motion, and overlaps appreciablywith the secondary peaks, which now have a high rela-tive intensity. This behavior gets more pronounced attemperatures above 1500 K, where now the central peakfor perpendicular motion also exhibits considerable lossof intensity in favor of the satellite peaks.

The behavior in temperature of G,(r, t) on the (110)surface differs considerably &om the corresponding oneon the more close packed (100) surface, thus suggest-ing a different scenario for disordering. First, on this

Ni(110)I

1000--

0.2—

00 I

I'

I Ii

0.4—

0.2

0.0s.' 0.4—

l

120

Il

II

Il

I

160

0,2

0.0

0.2

0.0—2 0 2

r (A)

-4 —2 0 2

r ()(.)

FIG. g. game as Fig. 7, but for the (110) slab.

surface, the onset of disordering (i.e. , nonvibrational mo-tion) takes place at the same T in all three Cartesiandirections, as is visible already in the T = 1000 K panelof Fig. 8. Second, adatoms are found to preferentiallydiffuse in a direction parallel to the rows of atoms on thetop layer (directed along z), manifest by the existence inG, (x, t) of structure at distances further away from thecentral peak than in the case of G, (I(, t) (e.g. , at 1000K); further, the smoother shape of G, (x, t) suggests thatdiffusion along x is almost continuous, while it proceedsby jumps along y; this will become clearer in the nextsection, where we analyze the trajectories de visu. Themagnitudes of the secondary peaks are small and subjectto large statistical Huctuations at this temperature, butincrease rapidly with increasing temperature. At 1300 K,the central peaks for lateral motion have lost much moreintensity than the corresponding one for perpendicularmotion, showing that self motion is much more rapid inplane than out of plane. At 1500 K, peaks start to ap-pear in the z direction at distances corresponding to twointerlayer spacings, signaling the formation of a secondadlayer, and the motion of atoms deep into the substrate.We would interpret this result (rather loosely) as the on-set of the roughening transition. At still higher T, therapid self-motion of the atoms results in the secondary-peak structure fading out almost completely, particularlyin the z direction; the central peaks, further, have lostconsiderable intensity, and long-distance tails develop inall three functions. Thus, even on a time scale as short as45 ps, the surface is in a liquidlike state, i.e., premolten.

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 091

Zhoj ectori es

Ni(110), T=900 K

+4444&4 f4-~~t4$~I~

0.5 -4 p 1 4 0 4 p 4 g +kg%% $44/+eyeyaQqe-"W444»y$4 4

, , +~ekeyg4~I

/

I(

I[

I/

I

1 0 5&%%&4% &4-

I I, I, I, I

—1.0 -0.5 0.0 0.5 1.0x (nm)

1.0 1 ~ 5z (nm)

FIG. 9. Trajectories of surface atoms for Ni(110) at 900 K,projected onto the (x, y), (z, z), and (y, z) planes.

A detailed description of the dynamical processes tak-ing place at the surface can, in principle, be obtainedfrom a careful examination of the trajectories of individ-ual atoms. It is difficult, in practice, to extract quantita-tive information from such a study, but it can neverthe-less provide useful insights on the behavior of the system.We present, in Figs. 9—12, in order to understand bettersome of the features already noted in G, (r, t), selectedtrajectory plots for atoms that belonged, initially (i.e., atthe beginning of equilibration), to the top layer of eachof our two slabs; the trajectories are projected onto theplanes (z, y) and (z, z) for (100), and onto (z, y), (z, z),and (y, z) for (110). In all cases, only the equilibriumportion of the runs are shown.

The behavior of the atoms in the vibrational-motionregime is illustrated by the trajectory plots of Fig. 9,for the (110) surface at 900 K. These "slices" show thatthe surface atoms undergo bounded motion around theirequilibrium positions, except for occasional excursions asubstantial distance away from the lattice point. Theslightly elongated shape of the spots in the (z, y) planereveals that the amplitude of the motion is larger in the ythan in the z direction, as expected, while, interestingly,it is about the same in the x and z directions; the latterresult will be discussed in more detail in Sec. III C. Theanharmonic character of the vibrational motion is clearlyvisible in the (y, z) slice: the trajectory spots exhibit ahighly asymmetric shape along the z direction and aresignificantly elongated along y; such spots, in a harmonicworld, would have inversion symmetry.

An example of a trajectory giving rise to the unex-pected secondary-peak structure observed in G, (r, t) forthe (100) surface in the temperature range 1200—1400 K,namely, in-plane diffusion in absence of out-of-plane mo-tion (cf. Fig. 7), is presented in Fig. 10. The (z, y) slice

1.0

0.5

0.0

—0.5

Ni(100), T=1300 K

~64WIp )v4f$i& sIL.&&

+~+It 4~%+-~'~ 4

%4@4~~4% 4 4%'g-

-~~~ 4~~=~»

2.0

1.5

laemvdSM

-~ '&p &~ ~ ~~ 'IW 1IP

I i I i I i I

-0.5 0.0 0.5 1.0x (nm)

FIG. 10. Trajectories of surface atoms for Ni(100) at 1300K, projected onto the (z, y) and (z, z) planes [(y, z) is equiv-alent to (z, z) by symmetryj.

shows the concerted "circular" motion of four atoms viathe temporary occupation of a hollow adatom site (centerleft). Similar events have been observed at other tem-peratures, for instance at 1200 K where one particularevent involved six atoms. Other types of concerted mo-tion involve atoms in a stringlike configuration, as seenagain in the (z, y) slice (lower right). Such events con-tribute significantly to G, (z, t) and G, (y, t), but negli-gibly to G, (z, t) since several atoms hop to new latticesites within the layer during the process, while the hol-low adatom site is occupied only temporarily by a singleatom.

The process of disordering and the nature of the cor-responding adatom motion on the (110) surface is illus-trated in Fig. 11. First, the (z, y) slice shows a numberof atoms in hollow adatom sites as well as subsurfacesites. The spots for subsurface atoms are evidently morelocalized (and more symmetric) than the correspondingspots for adatoms. The latter are elongated and streakedin the direction of the rows, indicating that difFusion ismuch easier along than across the rows; in fact, we findno event of diffusion across a row, except by an exchangemechanism involving an atom belonging to a row. Thisexplains why the peak structure in G, (r, t) is betterdefined along y than along x, i.e., the system becomes"liquidlike" in the latter before it does in the former.

At still higher temperatures, the surfaces are led intoa highly disordered state, even though the underlyingstructure of the lattice is visible at all temperatures ex-amined (i.e. , up to 1700 K). This is exemplified in Fig.12 for the (100) surface at 1600 K. Thus, even thoughdiffusion events are taking place at an extremely highrate, the layer structure can still clearly be identified;this is true even of the adlayer, where motion is fastest.The surface, therefore, can be said to be in a state ofdynamic disorder but static order.

12 092 YVON BEAUDET, LAURENT J. LEWIS, AND MATS PERSSON SO

Ni(110), T=1100 K

1.0

0.5

0.0

—0.5„ssi ~ij~l~1.5

1.0

II

I/

I[

I[

I

I i I i I i l i i

—1.0 -0.5 0.0 0.5 1.0

1.0 1.5z (nm)

x (nm)

FIG. 11. Same as Fig. S, but at 1100 K.

B. Thermal expansion

Ni(100), T=1600s, Me i

1.01 i i--'

. @&i I i

0 5E

'ti 1 ~

,-p~=

f

0 0 s~

Thermal expansion is a direct manifestation of anhar-monicity. For a surface, one expects the thermal expan-sion coefFicient to be different in the direction perpen-dicular to the surface than in the direction parallel toit. Indeed, the broken inversion symmetry at the surfacecauses an asymmetry in the potential seen by surfaceatoms, as was elegantly demonstrated by Ditlevsen et al.for the case of Cu(110) (Fig. 2 in Ref. 22): the potentialis quite generally expected to be softer in the region nearthe surface, thereby causing atoms to be displaced &omequilibrium further away than they are deep in the bulk(this can actually be seen &om Figs. 3 and 4, the peaksgetting broader near the surface) and allowing them toprobe more anharmonic regions of the potential. Thismotion will, of course, depend strongly on temperature,

1.30 Ni{110}

1.P.8

1.26

1 +24

10

8 —(b

and in general also on the particular surface under con-sideration [e.g. , (100) vs (110)]. Parallel to the surface,now, the lattice parameter is essentially that of the bulk.In our model, in fact, the bulk lattice parameter, whichvaries with temperature (cf. Fig. 2), is used as input tothe slab calculations, since it determines the structure ofthe 6xed layers.

For layered systems, in the direction perpendicularto the surface, thermal expansion is directly observablein the dependence on temperature of the interlayer dis-tance, E. Results for the (100) surface were given in Ref.24; they are shown in Fig. 13(a) for the (110) surface.More speciffcally, we show both Eg2 (the distance be-tween the top layer and the one just below), and E2s(between second and third layers) as well as, for com-parison, properly-scaled values for the bulk (cf. Fig. 2),

The corresponding thermal expansion coefBcients,a,~(T) = (1/E;~. )dE;~/dT, are displayed in Fig. 13(b).

While the thermal expansion coefFicient for the bulkshows only a weak temperature dependence, as predictedby quasiharmonic theory, even at values approachingmelting, the situation is quite different for the surfaces.For the (100) surface we reported, in Ref. 24, a thermalexpansion coefficient ni2 about 30% larger than n~ atT = OD, twice as large at 900 K, and 10 times that of thebulk at 1700 K. For the second interlayer spacing, on theother hand, we find o.23 to follow closely the bulk valuesat low temperatures, but to depart signi6cantly Rom a~at temperatures above 1200 K; at 1600 K, in fact, o.2q

becomes as large as a~2 at 1600 K; we will return to thispoint in our discussion of the layer structure factors.

For the (110) surface, now, Fig. 13, a relaxation in-

wards of the top layer is observed at all temperatures [ofthe order of 2.5 %, compared to, respectively, —5% and

—0.5

II

2.0I

'l 0

0 400 800T (K)

I i I

1200 1600

I i I i I

-0.5 0.0 0.5 1.0x (nrn)

FIG. 12. Same as Fig. 10, but at 1600 K.

FIG. 13. (a) Calculated interlayer distances Eq2 (squares)and E2s (triangles) vs temperature for the Ni(110) surface;also shown for comparison is the average interlayer distancefor bulk Ni (circles). (b) Corresponding thermal expansioncoefBcients. Note that Zq2 suHers an inward relaxation at alltemperatures except, perhaps, the highest one.

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 093

C. Mean-square amplitudes of vibration

In a strictly harmonic system, the mean-square ampli-tudes of vibration (MSAV), (u ), with a the spatial direc-tion of the motion, increase linearly with temperature, 8

20I

'I

'I

—9% from high-energy ion scattering (HEIS) and LEEDmeasurementsM], except perhaps at the highest one. Itis clear, however, upon examining Fig. 4, that the veryconcept of layers, for this surface, becomes doubtful atsuch high temperatures, even though some structure isstill visible in n(z); this results in large uncertainties onthe estimated layer positions, manifest in the scatter ofthe points in Fig. 13(a). Nevertheless, the enhancementwith temperature of the interlayer spacing Zz2, comparedto the bulk, is evident. We note here that, contrary tothe (100) surface, the anomalous thermal expansion is lo-calized to the topmost layer: 823 follows closely the bulkvalues at all temperatures. We have fitted fourth-orderpolynomials to the data up to 1400 K in order to cal-culate the thermal expansion coefFicient shown in Fig.13(b): clearly, o.i2 rises much faster than n~ or cr23.

We compare explicitely, in Fig. 14, the calculated ther-mal expansion coefficients for both surfaces, and for thebulk. It is seen that expansion is larger at the (110) sur-face, and the onset of "anomalous behavior" takes placeat lower temperature, than at the (100) surface, in linewith our observations of earlier disordering on (110).

The enhanced inward relaxation of the (110) surface,compared to (100), can be understood simply, again, interms of coordination: because of its reduced coordina-tion, the (110) surface suffers a "deficit" of electron den-

sity (relative to the ideal bulk case) and compensates bymoving the surface layer closer to the layer beneath it,as observed in our calculations. Such an argument, even

though it proves useful, is, however, only qualitative, as"glue" models (such as EAM) neglect the details of thedistribution in space of the atoms, which differ from onesurface to another. 3

so that the quantity P—:(u2)/T is a constant. Anydeviation from linearity, therefore, is a measure of anhar-monicity. The MSAV's are shown in Fig. 15 for the (110)surface [cf. Ref. 24 for the corresponding results on (100)]as well as for the bulk (averaged over the three directions,since they are equivalent). The calculated bulk MSAV'sare close to the measured values, viz. , 5.54 x 10 s A2 at400 K compared to 6.4 x 10 s A2 at T = 370 K frominelastic neutron scattering. It is important to note atthis point that the bulk (u )'s sufFer deviations from lin-earity that increase with temperature, so that even thebulk system is clearly anharmonic.

For surfaces, now, just as for the other properties stud-ied so far, the loss of coordination results in the restoringforces being less than in the bulk, and consequently themean-square displacements being larger. For the (100)surface, as reported in Ref. 24, the MSAV's are foundto be essentially isotropic for temperatures below 900 K,while above 900 K we observe (u~~) approximately 20%larger than (u2&), i.e., (u2&)/(u2~~) —0.8. Similar obser-vations have been made in the MD study of Cu(100)mentioned above. On the experimental side, in-planeamplitudes larger than out-of-plane ones by about 30%have been observed in ion-back scattering measurementsof the Cu(100) surface, but at relatively lower T com-pared to the bulk Debye temperature. To our knowl-edge, MSAV's have not been measured on Ni surfaces,and it would be of considerable interest that experimen-tal verification of our predictions be performed.

The situation is somewhat difFerent on the (110) sur-face, where the anisotropy develops already at the lowesttemperatures. Just as with the (100) surface, in-planeamplitudes are larger than out-of-plane amplitudes, atleast at high temperatures. In contrast to the (100) sur-face, however, in-plane MSAV's are no longer equivalentalong the x and y directions, because of the asymmetryof the surface lattice vectors ai and a2 (cf. Fig. 1). Thus,vibrations along the rows ((u2)) are substantially smallerthan across the rows ((u„)). This was already visible inthe trajectory plots of Sec. IIIA4.

Our results for the perpendicular vibrations are in ex-cellent agreement with those of Cao and Conrad who ob-serve, in the low-temperature regime (T ( 900 K), P~ =

15

10

f--"-I444+++++++i I i I ( I

0.06

0.04

V0.02

I

Ni(110) II

Ir/ "ej

..0~ 'i--g ..Q

a -0g. .y e.-8

e=",--0 'F

400 800 1200 1600T (K)

FIG. 14. Comparison between calculated thermal expan-sion coefBcients for the (100) surface (squares), the (110) sur-face (triangles), and the bulk (circles).

0 400 800 1200 1600

FIG. 15. Mean-square amplitudes of vibration (u ), (u„),and (u, ) (611ed triangles, squares, and circles, respectively)vs temperature for the (110) surface, compared to the bulk(open circles —averaged over the three Cartesian directions).

12 094 YVON BEAUDET, LAURENT J. LEWIS, AND MATS PERSSON SO

2.7 x 10 A K, compared to our 2.8 x 10 A2 K(averaged over temperatures below 900 K). Whether thisagreement is fortuitous or not is unclear: our results, infact, disagree strongly with Gao and Conrad's for the in-plane MSAV's. While we find P~ = 3.0 x 10 s A2 Kand P„=4.0 x 10 A K, Cao and Conrad report val-

ues of 17 and 15x10 A K, repectively, i.e. , approx-imately six times larger than the out-of-plane MSAV's,and five times larger than our own values. Yang andRahman, in contrast, find, for the Cu(110) surface,(uz)/(u„) = 0.6, quite comparable to our value of 0.7.We cannot explain the discrepancy between calculationsand experiment.

It remains that, on both (100) and (110) surfaces, in-plane MSAV's dominate out-of-plane MSAV's. This is aneffect of the many-body nature of the interatomic forcesand cannot generally be understood in terms of a simplepairwise interaction.

to K and to the surface); thus, v,. is the projection of v,onto the nth polarization vector.

The Fourier transform of Ce (K, t) yields the densityof phonon states (at fixed wave vector):

CP(K, ~) = (2m)' f CP(K, t)e 'dt'. (3.4)

+e (K ~) = 2~1"e (K ~)l (3.5)

where ve (K, ur) is the Fourier transform of ve (K, t):

Yang et a/. have remarked that the above expressionmay lead to negative-&equency peaks in the phonon spec-tra, because of the finite length of the MD runs, i.e.,truncation errors; this is a problem, in particular, at lowtemperatures where correlation times are so long that thecorresponding functions do not rapidly approach zero.They propose, instead, to use the following approxima-tion:

D. Surface phonons+OO

ve (K, ur) = (2z) 've (K, t)e'"'dt. (3.6)

In addition to thermal expansion and "anomalous" am-plitudes of vibration, anharmonicity infIuences markedlythe dynamics of the lattice. Anharmonic terms in the in-teratomic potentials result in the phonons being coupledto one another, to a degree which increases with tem-perature, in turn causing energy to be exchanged andphase coherency to be lost. Phonons, therefore, acquirea finite width in &equency space, and are also shiftedin energy, as temperature increases. A study of this be-havior in temperature is therefore of direct relevance tounderstanding anharmonicity.

Method of calculation

In a molecular-dynamics simulation, informationabout the lattice dynamics is implicitely included in thecorrelations in the motion of atoms, and can thereforebe obtained from quantities such as the velocity (or cur-rent) autocorrelation function. Specializing to the case ofa two-dimensional system, e.g. , here, the layer designatedby the index I and containing Ng atoms (determined dy-namically), the velocity autocorrelation function is de-fined, for momentum transfer K = (K,K„,K, = 0)(i.e. , in the surface), as the ensemble average:

Ce (K, t) = {ve (K, t)ve (—K, 0)), (3 2)

where the reciprocal-space, time-dependent, velocity isgiven by

Nc

ve (K, t) = ) v, (t)e ' (3.3)

with R;(t) is the position of atom i at time t In the.above, o. is the polarization of the phonon mode withrespect to K, which we choose to be either 'ongitudinal(parallel to K), shear horizontal (perpendicular to K, inthe plane of the surface), or shear vertical (perpendicular

It is easy to demonstrate that (3.4) and (3.5) are strictlyidentical in the limit of infinite run times; following Yanget al , we h.ave used (3.5).in the calculations presentedin the next sections.

The spectra were in all cases calculated using the fulllength of the runs. In order to simulate —to some ex-tent —the experimental situation, and to eliminate er-rors arising &om the use of a finite time integration, theraw spectra were convoluted with a Gaussian functionhaving a full width at half maximum of 0.2 THz. Thus,the widths of the peaks result in part from anharmoniceffects, and in part &om the convolution procedure.

2. Lou-tempemtum spectv a

We present, in this section, the calculated phonon spec-tra for our two surfaces in the low-temperature limit(100 and 300 K). At these temperatures, the thermaldisplacements of the atoms are so small that the inter-atomic forces are dominated by the harmonic terms; at300 K, further, i.e., close to the Debye temperature forNi, quantum effects are small and the comparison of thecalculated spectra with experiment will be more accu-rate. All numbers quoted in this section will thereforerefer to the 300 K spectra; the temperature dependenceof the phonons is discussed in the next section. Since we

are interested in the behavior of the surface phonons lo-calized to the outermost layers, we concentrate our effortson the two-dimensional —Brillouin-zone boundary pointsX and M for the (100) surface, and X and Y for the(110) surface, as indicated in Figs. 1(b) and 1(d).

The low-temperature spectra for the M and X pointsof the (100) surface Brillouin zone are displayed in the

upper two panels of Figs. 16 and 17, respectively, while

Fig. 18 shows the 100 K spectra for the second layer, i.e.,just beneath the surface.

At the M point of the top layer, Fig. 16 reveals quite

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 095

clearly the shear vertical Sq(z) mode at 4.88 THz, locatedfar from the projected bulk band edge (at 6 THz),and the degenerate horizontal surface resonance Rq [oftenlabeled Lq(z + y) and Lz(x —y)j at 6.7 THz. In thesecond layer, Fig. 18(a), the vertical Sz(z2) mode is seenat 6.53 THz, as well as the Lj and L~ modes at 6.7 THz,polarized in the same direction as the corresponding toplayer modes.

At the X point, now, Fig. 17, we see the shear hori-zontal SL(y) mode at 2.62 THz, the shear vertical S4(z)mode at 3.85 THz, and the longitudinal Ss(z) mode at7.60 THz. The corresponding spectrum for the secondlayer is shown in 18(b), and we observe changes in thepolarizations of the modes as compared to the surface:the S4 mode (at 3.85 THz) is now longitudinally polar-ized, while the Ss mode (at 7.60 THz) becomes verti-cal. Additionnally we see, in the second layer, the shearhorizontal Ss(y) mode at 4.2 THz, which is symmetryforbidden in the 6rst layer; by the same argument, theSq mode is absent in the second layer. We also clearlysee the vertical resonance R4 at 4.65 THz (longitudinallypolarized in the top layer, but hardly visible).

A comparison of our calculated &equencies with theresults of electron energy-loss spectroscopy measure-ments ' is provided in Table I. Our results are quite

0.012

0.008

0.004

0.000

3 0.008

0.004

O, OOQ

0,008~ I~4

0.004A

0.000

0,008

0.004

0.000

, ~ I

~ ~

~ ~

~ ~

II

II~ ~

~ ~ II~ ~ I~ ~

~ & I~ ~

I~ ~

I I

I I

I.'I

I~ I

II

'e

I', ll

~ I, .r =.I

I~ ' I

II

~ ~

~ ~

~ ~

~ ~

~ ~

~ ~

~ ~

~~

~ ~

~ ~

~ ~

~ ~

~ ~

~ ~

100

II

I

300

I[

I

600

II

II

I

~ ~

~~

~ ~

~ ~

~ ~

~ ~

~ ~

~ ~

~~

~~

II

I

1200

~ I

Ni(100) —X

0.012

0.008

0.004

0.000

Ni(100) —M

I I" I'

I

100I

II

I

II

I

I

I'

I'

I

II3OO

0.012 II

II

II

II

I

& (a) Ni(100) —M

0 2 4 6 8 10u (THz)

FIG. 17. Same as Fig. 16, but at point I of the first Bril-louin zone of the (100) surface.

3 0.008—

~ 0004—N4P

0.000N

0.008~&Na 0.004—QA

0.000

0.008

I'

I'

I

I

.II

I

II

e

\r' L i JI

II

II

1200

0.008

30.004—

U

0.000

0

~ ~t 0.008A

IIII

III

I II

I II

I

I

I [ I Pl

~ ~

~ ~

~ ~

~ ~

~ ~

~ ~

I~ ~

~ ~

(b) Ni(100) —X

~ ~

II

~ 'I

II

[I

II

0.004

0 000 I~ = . L I w j

0 2 4 6 8 10u (THz)

0.004

:II'.I I

4I4 I

~ ~

I I~ I I

~ ~ ILI LI

FIG. 16. Phonon spectrum at point M of the first Bril-louin zone of the (100) surface at selected temperatures, asindicated. The full, dotted, and dashed lines correspond to,respectively, longitudinal, shear horizontal, and shear verti-cal component of the spectra (see text). The units for thedensities of states are arbitrary.

6 8u (THz)

10 12

FIG. 18. Phonon spectrum for the aecond layer of the (100)slab at points (a) I and (b) X of the surface Brillouin zoneat 100 K.

12 096 YVON BEAUDET, LAURENT J. LEWIS, AND MATS PERSSON 50

ModeS4(z)Se(x)

MD Expt. Mode MD Expt.3.85 3.96 + 0.06 Si(z) 4.88 4.65 6 0.157.60 7.56 + 0.15 Sq(zz) 6.53 6.6

Ni(110)

Mode

S~(z)Sz(y)S~(x)MSp(x)

XMD3.944.807.434.18

Expt.4.115.076.84.26

YMode MD

Sg(y) 2.18S,(*) 3.07Ss(z) 3.59S,(y) 5.75

Expt. '2.4

3.483.7

References 41 and 42.References 43 and 44.

generally in excellent agreement with the experimentalvalues, mostly within the error bars from both numericaland experimental data, thus reassuring us of the valid-ity of our MD-EAM-slab model —at least in the (low-temperature) harmonic limit.

Moving to the (110) surface, the results are displayedin Figs. 19 and 20 for the X and Y' points of the Brillouin

TABLE I. Comparison between calculated (MD) and mea-sured (Expt. ) surface phonon frequencies, as indicated. Allfrequencies are at 300 K, and are measured in THz. The po-larization of each mode is given in parenthesis (see text); allmodes are for the top layer, except z2, which is for the secondlayer.

Ni(100)X

zone, respectively. For the X point, we find the verticalSq(z) mode at 3.94 THz, the longitudinal MSo(x) mode4.18 THz, the horizontal Sz(y) mode at 4.80 THz, andthe longitudinal S7(x) mode at 7.43 THz. For the Ypoint, we have the longitudinal Sq(y) mode at 2.18 THz,the horizontal Sz(x) mode at 3.07 THz, the vertical Ss(z)mode at 3.59 THz, and the longitudinal Ss(y) mode at5.75 THz. Our values are compared with experimentin Table I and again confirm the validity of our model.

tA'e proceed now to a discussion of the dependence intemperature of the phonon frequencies and linewidths,i.e. , anharmonic efFects.

8. Temper ature dependence

In order to illustrate the presence of enhanced an-harmonicity as a function of temperature, we present,in Figs. 16, 17, 19, and 20, the spectra deduced fromour runs at temperatures of 600 and 1200 K. The pre-cise evolution of the &equencies and widths for severalof the modes discussed above are given in Figs. 21 and22. There are evidently uncertainties in the values calcu-lated arising, for instance, &om the difFiculty in properlyde6ning the average &equency for a line that may con-tain smaller features, as well as the width when morethan one line overlap. The trends, nevertheless, makeno doubt: frequencies decrease and widths increase withtemperature, and we estimate the error in our results tobe of the order of the scatter in the points relative to asmooth behavior. Note again that, because of the arti-

0.008

0.004

0.000 II

3 0.008

0.004

0,000

0.008~ W

0.0040Gh

lI

I~ 'I

I

~1

II

II~ ~

I

I

II

300

600

Ni(110) — X0.012

100IIII

0.012

0.008—

0.004

0.000

0.008—hC

0.004NQ

0.000

0.008—

0.004—Q

Ni(110) —Y

I'

I'

I

100I

I

~ I

'I

I—~;~~L &r iI

' 'I

~ ~

~ I

I I

h~ ~ II~ I

~ ~

; ] I

( I

I', \

II

II

II

II

600~ ~ II

0.000

0.0081200

0.000 II

II

II

II

I

O.OOB

0.004 0.004

0.0000 2 4 6 B 10

~ (THz)

FIG. 19. Same as Fig. 16, but at point X of the first Bril-louin zone of the (110) surface.

0.0000 2 4 6 8 10

u (THz)

FIG. 20. Same as Fig. 16, but at point Y of the first Bril-louin zone of the (110) surface.

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 097

Ni(100)I

'I

'I

9-- -e —e e-e-e-y-e-e e q&

I

12000 4 M I

0 400

6cial convolution function used for smoothing the spec-tra, the widths of the phonon lines cannot be compareddirectly to experiment, but their retatiee increase withtemperature is a measure of enhanced anharmonicity.

It is of interest to determine whether there are trendsin the temperature-dependent frequencies and widths of

Ni(110)

800T (K)

FIQ. 21. Temperature dependence of (a) the &equeucy ~and (b) the width b, of the various phouous calculated atthe (100) surface (cf. 'Ihble II): Sz(y) (empty circles), S4(z)(filled circles), Ss(z) (filled triangles), L', (x —y) (empty tri-angles), Sq(z) (empty squares), and Sz(zz2) (filled squares);the widths are not calculated for all phonons.

the different modes with respect to polarization and crys-tal face. In order to investigate this, we present in Ta-ble II the temperature-averaged slopes of both quantitiesfor each mode considered, assuming an approximately-verified linear behavior (as expected from perturbativetreatments, Ref. 45). We note first that there are nopronounced difFerences in the magnitudes of d luego/dTand dA/dT between the two different crystal faces. Fur-ther, on the (100) face, there is no systematic behav-ior among the different polarizations, while on the (110)face the shifts and widths are larger for modes polarizedin the horizontal direction than in the vertical direction.There seems to be no obvious correlation between themagnitudes of the shifts and the widths; however, anenhanced d inn/dT for one mode does not necessarilyimply a larger db, /dT. On the (100) face, all the db, /dThave about the same value, while on the (110) face thewidths vary substantially between diH'erent modes. Inparticular, the S2(y) and the Sv(z) modes at the X pointhave the largest dA/dT.

In order to assess the relative importance of sur-face versus bulk anharmonicity, we compare the surfaced ln~/dT with its bulk counterpart. (As far as we know,there is no experimental information available about theT dependence of b, in the bulk. ) We have not carried outa separate MD calculation of d ln u/dT for bulk phonons;however, an estimate of it can be obtained indirectly froman analysis of the calculated temperature dependence ofthe MSAV's. This analysis is based on the connectionbetween vibrational amplitudes and phonon &equenciesin a quasiharmonic description. As we have remarkedearlier, MSAV's, in the strictly harmonic case, increaselinearly with temperature, i.e., for direction a and layer 8,P g(T) = (u &2( T))/ T~ /3~&, a constant independent ofT: any deviation is a direct manifestation of anharmonic-ity. We may use, then, the slope of P g, d/3 g(T)/dT,as a measure of anharmonicity. Within quasiharmonictheory, which incorporates anharmonic efFects througha dependence on temperature of the phonon frequenciesvia an assumed dependence on the equilibrium positions,we have for a layer in the bulk, as demonstrated in theAppendix:

~ ~ ~ ~ ~-CL ~-

I

1.6 — (b)H

i 1.2

TABLE II. Slopes of the frequencies and linewidths, as-suming a linear temperature dependence, for the variousmodes considered, as indicated. d la ~/dT is in units of 10K ' and dA/dT in units of 10 THz K

Ni(100)

0.8

0.4

0 200 400 600 800 1000T (K)

FIG. 22. Same as Fig. 21, but for the (110) surface:Sz(z) (filled squares), Sz(y) (empty triangles), Sq(y) (filledhexagons), Sz(z) (empty squares), Sz(x) (empty circles),S3(z} (empty hexagons), MSO (x} (filled triangles}, aud Ss (y)(filled circles).

Mode

S~(y)S4(z)s6(2:}Ni(110)

ModeS~(z)Sz(y)S~(*)MSO(x)

Xd tucks/dT

-0.9-0.9-0.8-1.1

dA/dT3.311.11.2.4

Xd ln ~/dT dB/dT

-0.8 7.8-1.4-0.7

Mod.eLi(* —y)

Sg(z)S2(z2}

7.56.9

Modes (ySg(z)S3(z)Ss(y)

Yd lu~/dT dZ/dT

-1.5 2.0-1.7 4.3-0.9 2.0-1.5

Id ln ~/dT d6/dT

-1.3-0.9-1.0

12 098 YVON BEAUDET, LAURENT J. LEWIS, AND MATS PERSSON 50

dPg d ln erg

dT dTAB)

(3.7)with our observations, in Sec. IIIC, that the MSAV's atthe surface are roughly 2—5 times larger than in the bulk.

Here, A labels a phonon mode with &equency up, and

()q~ is an average over the phonon modes in the bulk.The high symmetry of the bulk lattice results in theMSAV's being isotropic, and in effect P q = PB. Atthe surface now, because the symmetry is broken, con-tributions to dP r(T)/dT possibly arise from changes ofpolarizations on going &om one layer to another, or fromthe difFerent Cartesian directions. Such individual con-tributions vanish, however, if one considers, rather, thesum of deviations from the bulk, AP = P&(Ps —P~),where Pt is the average of P r over the three Cartesiandirections. As shown in the Appendix, we have then

db, P — ~ d ln urp

dT dT i.t)d in~,

(3.8)

with ()p s an average over the phonon modes weightedover the polarization vectors in the layer 8. In a quasihar-monic theory, the changes in the equilibrium distances isdriven by the decrease in the free energy with decreas-ing phonon frequencies, so that a thermal expansion isaccompanied by a decrease of the frequencies with tem-perature, as we have just verified. As a result —andnot surprisingly —dP(T)/dT is positive in the bulk, aswas indeed found in Sec. III C. In fact, we also find thatdb, P/dT is positive, which suggests that either the f're-

quency shifts or the mean-square displacements are en-hanced at the surface.

Equation (3.7) yields, for the bulk, at temperaturesbelow 900 K, using the MSAV's of Sec. III C (from whichwe can extract P~ and dP~/dT), (d 1nurp/dT)g~—1.3 x 10 K, a value which is in excellent agreementwith the corresponding experimental number, —1.4 x10-4 K-'.4'

Now, the value of d 1nu/dT for the different surfacemodes are, it turns out, very much comparable to thisaverage bulk value, as can be seen in Table II. Evi-dently, therefore, there is no systematic or pronouncedenhancement of the frequency shifts at surfaces. Thecalculated values are, in fact, consistent with the re-sults of a recent analysis by Benedek and Toenniesof the temperature-dependent shifts of bulk and surfacephonons for Al and Cu, which reveals that the surfaceanharmonicity, as gauged by d ln&up/dT, is at most ofthe order of twice the bulk anharmonicity.

Our calculations demonstrate, therefore, that surfacephonons behave in a way entirely similar to those in thebulk, i.e. , are not particularly manifest of anharmonic ef-fects. Referring to Eq. (3.8), we must conclude, that theobserved anharmonic "anomalies" at the surface, as mea-sured by dP(T)/dT, are almost completely contained inthe prefactor P& in Eq (3.8), esse. ntially the MSAV's, andare therefore not related to an extraordinary or anoma-lous softening of the phonons. This is perfectly consistent

E. Layer structure factors

4. LEED intenaitie8

Because, in a crystal, atomic vibrations break, to someextent, the periodicity of the lattice, difFraction experi-ments provide essentially direct information about theequilibrium dynamics of the atoms, more specificallyabout their vibrational amplitudes. In a kinematic the-ory of difFractive scattering, the elastic intensity for mo-mentum transfer K, I(K), is proportional to the elasticpart of the dynamic structure factor, S'(K)—:S(K, m =0); the latter can be expressed as o

S"(K) = l(~(K)) I' (3.9)

where p(K) = (1/N) P,. i e ' ' ' is the Fourier trans-form of the atomic density of the system (which containsN atoms), and () is an ensemble average. In the har-monic limit, S"(K) reduces to the usual Debye-Wailerexpression,

S (K) = exp[ —((K u) )]. (3.10)

In this limit, therefore, the difFracted intensity dependson temperature directly (and only) through the Debye-Waller factor:

I(K) = Io(K) exp[—((K u) )], (3.11)

where Io(K) is the elastic intensity of the rigid (classi-cal zero-temperature) lattice. It is useful, in the presentcontext, to define a layer structure factor (LSF), S&i(K),related to the Fourier transform of the atomic density inthe first layer, pt(K), as in (3.9).

Cao and Conrad have used LEED to probe the dy-namics of both the (100) and the (110) surfaces of Ni. i 2

Because of the multiple scattering processes of the elec-trons &om the ion cores, a detailed interpretation ofthe LEED spectra is in general diKcult. In an efFort

to achieve this, we have simulated the measured spec-tra within the framework of the kinematic theory pre-sented above. More precisely, noting the extreme surfacesensitivity of LEED, we have assumed that thermal at-tenuation arises solely &om the atoms in the first (top)layer (again defined dynamically), and accordingly con-

sider the Si"(K) structure factor in the discussion below.Our evaluation of the LEED intensities for the (110)

surface is shown in Fig. 23 [cf. Ref. 24 for the correspond-ing results on (100)] for the same moxnentum-transfer val-

ues as in the experiments; here the momentum-transfervectors, which we label (ki, k2, ks), are expressed in termsof the reciprocal-lattice basis vectors (see Ref. 52), i.e. ,

K = kibi + k2b2 + ksbs, thus, we have (002), (001),(02&), (022), (102), and (102). In order to make contactwith experiment more closely, we have renormalized thecalculated Sii(K) by the strictly harmonic contribution

exp[ —K P (K)T], where P (K) = (1/IC ) limT~o((u.K)2)/T is independent of temperature. The renormal-

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 099

I)

I/

If

I

'g+ E'o0.6 W ~s'.'g

0.4 D

p 2 —Ni(110)a%

) I i I i I i

0 400 00 i200 1600

I

C4

Q

8T (K)

FlG. 23. Elastic part of the surface dynamic structure fac-

tor, S"(K},renormalized by the harmonic Debye-Wailer fac-

tor, exp[ —K P (K}j (see text), at the (110) surface, for thesame six wave vectors as in Ref. 1, namely (001) (filled circles),(002) (empty circles), (023) (filled squares), (022) (emptysquares), (103) (empty triangles), and (102) (filled triangles).A purely harmonic theory yields a value of unity, indicated

by the straight dashed lines.

Layer dependence

We consider here the layer structure factor in moredetail in order to gain insight into the behavior of the slab

ized Sii(K), therefore, should equal unity in the strictlyharmonic regime, and deviations signal the presence ofanharmonic contributions.

In a quasiharmonic approxixnation, one considers theharmonic expression (3.10) to be valid, but the harmonic(u2), which varies linearly with T, is replaced by the ac-tual (u2), which includes anharmonic terms. The resultsof our calculations (for both surfaces) in this approxima-tion are displayed in Fig. 24, where we have used theMSAV's of Sec. III C.

As can be seen from these results, both surfaces areessentially harmonic at low temperatures, while signif-icant deviations develop as temperature increases, sim-ulating quite precisely the behavior observed in exper-iment. (Note that the "anisotropy" in the intensitiesarises principally &om the wave vector dependence of theDebye-Wailer factor. ) Most important, we observe thosedeviations, i.e., anharmonic manifestations, to take placewell before disordering —difFusion and adatom-vacancyformation —which, as we have seen in Sec. III A, appearsat temperatures of about 1200 and 1000 K for the (100)and (110)surfaces, respectively. Thus, in agreement withCao and Conrad, we find that the large attenuations ofthe LEED spectra originate primarily &om anharmonicbehavior, and not &om structural disorder.

Experimentally, the quasiharmonic approximation isutilized in order to deterxnine the MSAV's Rom thediffracted intensities. A comparison between the resultsof Figs. 23 (cf. also Ref. 24) and 24 shows, indeed, thatthe anharxnonic character of the structure factor is wellreproduced by the quasiharmonic approximation, exceptmaybe with some underestimation at temperatures closeto the onset of in-layer disordering. The quality of theagreement provides support to our analysis, in Sec III C,of the MSAV's, itself based on the quasi-harmonic ap-proximation.

I'

I'

t100.8

060.4

p p (a) Ni(100)

0.0C4

Q

V)

I'

I'

I

0.8 — '-&q& o.

0.6— gaQ~ KI

0.4

0 o (b) Ni(110)

0.0 I i I i I i I

400 800 1200 i600T (K)

FIG. 24. (a) Same as Fig. 23, but using the quasiharmonicapproximation (see text) for both surfaces. We also give in

(b) the results for the (100) surface, where we have used thesame four wave vectors as in Ref. 2, namely (transformedinto our coordinate system} (003) (filled circles), (004) (filled

triangles), (103) (empty circles), and (114) (empty triangles).

as a whole under increasing temperature. More precisely,we examine the layer (and temperature) dependence ofthe LSF for several inequivalent reciprocal-lattice vectorson the two surfaces, namely (expressed, again, in terms ofour basis vectors, Ref. 52), for the (100) surface, K-+„- ——

(110), Kioo = (023), and Kioo = (001), and, for the

(110) surface, Kiio = (201), Kiio = (021) and Kiio =(001). In order to eliminate the efFect arising Rom thelengths of the lattice vectors when comparing differentdirections, as well as the weak temperature dependenceof the lattice parameter, the logarithms of the I SF's willbe reduced by the factor W(K) = ~KI2/( —)2.

We show, in Figs. 25 and 26, the reduced LSF's forthe top four layers of each slab, as well as for the seventhlayer, which is sufficiently distant from both the surfaceand the static substrate so as to approximate (as we willsee) the behavior of the bulk. In all cases, the attenua-tion is first approximately linear in temperature. Higher-order corrections appear at intermediate temperatures;the vibrational motion of the ordered layers is responsiblefor these contributions. At still higher temperatures, avery rapid drop is observed, related to the appearance ofmigration processes; at this point, therefore, disorderingprocesses cross over &om being vibrational to diffusionalin nature, as demonstrated by the behavior of the vanHove self-correlation function in Sec. IIIA3. Evidently,the effect is more pronounced for layers at or near thesurface {the system is "premelting"), and is essentiallyabsent, in our model, at a depth corresponding to theseventh atomic plane, whose LSF varies more or less lin-early over the whole texnperature range. This indicatesthat the highest temperature examined here, 1700 K, isbelow the melting point for the bulk. (Indeed, bulk melt-

12 100 YVON BEAUDET, LAURENT J. LE%'IS, AND MATS PERSSON 50

Ni(100)I

743

Ni(110)I

7

(a) (0,2.—3)I, II I

I

(a) (

8g —1

av7

-S —(b) (1,1,

I I I I I

1I

I

7432

—1

20v7

4

7432

-3 —(c) (0.0.1)

I i I i I i I

400 800 1200 1600T (K)

FIG. 25. Natural logarithm of the reduced layer structurefactor S'(K) at the (100) surface as a function of tempera-ture, for the @rave vectors indicated. The number attached toeach curve designates the layer number, 1 being the topmostone.

ing occurs at 1740 K in the EAM version used here. ) Itis interesting to note, in Fig. 26(c), at low and interme-diate temperatures, that layers seem to pair up along thez direction, and, in particular, the top two layers behavein a distinct manner already at the lowest temperatures.This behavior is most certainly related to the fact thaton the (110) surface, in contrast to (100), both the topand second layer suffer a loss of coordination.

%e can extract, &om Figs. 25 and 26, the I SF's for thetop layer only of each surface and, in each case, for thethree nonequivalent directions considered; this is done inFig. 27. The function behaves in a way entirely consistentwith the surface MSAV's, as expected &om the quasihar-monic approximation. Thus, for the (100) surface, theLSF is isotropic at low temperature, and the "degener-acy" is lifted at about 900 K, with the amplitudes along xand y clearly larger than along z. For the (110) surface,now, the anisotropy is visible already at low tempera-tures, re6ecting the fact that the interatomic distancealong y is larger than along the other two directions (uvs a/v 2).

It is informative to consider again the I SF's, but nowplotted as a function of layer number (i.e., depth) for se-lected temperatures. We do this in Fig. 28 for our twosurfaces at four different temperatures: low (500 K); in-termediate (1000 K); and high (1500 and 1700 K); thetop layer is numbered "l." At low and intermediate tem-

-3 —(c) (0,0,1)

I i I i I

400 800 1200 1600T (K)

FIG. 26. Same as Fig. 25, but for the (110) surface.

a -1—2

V)

(a) Ni(100)I i I

w~m'% ~

~ -1-

V)

(b) Ni(11

0 400 800 1200 1600T (K)

FIG. 27. Structure factor for the top layer only of (a) the(100) slab, and (b) the (110) slab. The full, dotted, anddashed lines correspond to, respectively, the K„-,K-+-, and

K; Bragg points for the (100) surface, and K-, K„-,andK; for the (110) surface.

peratures, disorder is essentially constant throughout theslab, except for a sudden increase (i.e. , a drop in the LSF)at the erst layer. For the K,-~ wave vector, as noted ear-lier, the attenuation also affects significantly the secondlayer, which evidently couples strongly to the first one,as noted above. For the (100) surface, disorder normal to

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 101

I

0 50010001500

(a) Ni(100)1700—

I I I I I I I I I I I

iI

II

II

II

'I

—.- 5001000

—(b) 1700—

10 8 6 4 2 0

Layer 8

FIG. 28. Layer structure factor as a function of layer num-

ber at various temperatures (as indicated) for (a) the (100)surface, snd (b) the (110) surface. Layer 1 is the topmost

(surface) layer. The full, dotted, snd dashed lines correspond

to, respectively, the K„-, K-+„-, and K; Bragg points for

the (100) surface, and K-, K„-, and K; for the (110)surface.

IV. SUMMARY AND CONCLUSIONS

We have used molecular dynamics and interatomic po-tentials from the embedded-atom method to study thestatic and dynamic properties of the (100) and (110) sur-faces of Ni in a range of temperatures extending all theway to the bulk melting point. In both cases, our calcu-lations reveal the existence of two temperature regimes,which we have characterized by a detailed investigationof layer densities, adatom-vacancy —pair formation, vanHove self-correlations, particle trajectories, thermal ex-pansion, vibrational amplitudes, surface phonons, andlayer structure factors. First, at low and intermedi-ate temperatures, up to about 1200 K and 1000 K forthe (100) and (110) surfaces, respectively, we have the"anharmonic-surface" regime, where the amplitudes ofvibration are such that new, unexplored, and increas-ingly anharmonic regions of the interatomic potentialare being probed. Upon increasing the temperature fur-ther, the displacements become so large that migrationprocesses become possible, opening the way to a "dis-

the surface is always less than in the directions parallelto it.

It is clear &om this study of LSF's, that the tempera-ture at which a layer is struck by a given amount of dis-order, or at which substantial deviations &om the normalharmonic behavior take place, increases upon going awayfrom the surface into the bulk. This is direct evidence forenhanced anharmonicity at the surface in comparison tothe bulk.

ordered" surface regime, where the dynamics is domi-nated by in-plane diffusion and the formation of adatom-vacancy pairs. "Premelting, " in particular, for the (110)surface, gradually reaches deeper in the slab. At thesehigh temperatures, it is no longer possible to define sur-face phonons.

Our calculations not only reveal the adequacy of ourMD-EAM-slab model to describe the problem at hand,but also provides important insights in the dynamics ofthe surfaces at elevated temperatures, and, in particular,explicit manifestations of anharmonicity. First, the sur-face thermal expansion coefBcients are significantly largerthan their bulk counterpart, increasing dramatically withtemperature; in the "disordering" regime, they becomeanomalously large, in good agreement with the recentLEED results of Cao and Conrad. The mean-square am-plitudes of vibration are also found to increase anharmon-ically with temperature, and to be much larger than inthe bulk. We observe, in addition, out-of-plane displace-ments to be smaller than in-plane ones on both surfaces,in particular, at temperatures close to disordering. Theactual values of the MSAV's, we find, are in good accordwith the values reported by Cao and Conrad, except forin-plane ones on the (110) surface, a discrepancy we can-not explain at present. Our analysis of the layer structurefactors, in particular, within the &amework of a quasi-harmonic approximation, supports the experimental in-terpretation of the LEED peak intensities in terms ofenhanced anharmonic mean-square displacements ratherthan disordering. The enhanced mean-square displace-ments at high temperature are found to be a consequenceprincipally of the larger magnitude that can be affordedby surface atoms, because of the reduced coordination incomparison to the bulk. In contrast, the phonon &equen-cies and widths exhibit no pronounced softening, whichsuggests that collective excitations at the surface are af-fected by anharmonicity in a nontrivial way: we observean essentially linear dependence of surface phonon &e-quencies and widths with temperature. While this tem-perature dependence is not known &om experiment, ourlow-temperature values are in excellent agreement withthe electron energy loss spectroscopy values. It wouldtherefore be of considerable interest to carry out a de-tailed experimental study of the temperature dependenceof phonon line shapes, as well as mean-square displace-ments.

The crossover from the anharmonic to the disorderedregime is manifestly difkrent on the two surfaces. Thedisordering involves essentially only the outermost layeron the (100) surface, while on the (110) surface the sec-ond layer also participates; the latter reBects diHerencesin atomic coordination on the two surfaces. Our analysisreveals that disordering on the (100) surface first occursin the plane, by a mechanism involving the concertedmotion of several atoms with the transient occupation ofan adatom site. The out-of-plane disorder, by adatom-vacancy pair formation, precipitates only at higher tem-peratures. The activation energy for the adatom-vacancyformation is much lower on the (110) surface and, asa result, disordering by in-plane diffusion and adatom-vacancy formation take place both concurrently. The in-

12 102 YVON BEAUDET, LAURENT J. LEWIS, AND MATS PERSSON 50

plane diffusion occurs predominantly in the channels be-tween the rows in the adlayer, while diffusion across therows proceeds via an exchange mechanism with the atomsin the rows. The observed concentration of adatoms onthis surface is in remarkable accord with the numbers ob-tained experimentally by Cao and Conrad. Furthermore,the observed "roughening" temperature obtained in theirexperiments is essentially the same at which we observeout-of-plane disorder to strike the second layer.

In summary, our calculations support the experimen-tal interpretation of LEED intensities, before the onsetof disordering, in terms of enhanced mean-square dis-placements. At higher temperatures, anharmonicity be-comes even more pronounced, but disordering accountsfor a large portion of the observed effects. The enhance-ment of anharmonicity at the surface, however, is notmanifest in all properties. While, for instance, surfacethermal expansion is "anomalous, " no pronounced soft-ening of the surface phonons (compared to bulk ones)at elevated temperatures is observed. The increased vi-brational amplitudes of surface atoms results from theirreduced coordination, allowing them to probe more an-harmonic regions of the potential.

ACKNOWLEDGMENTS

(AI)

where m is the mass of the atom and k~ Boltzmann'sconstant. In a strictly harmonic theory, the amplitudessquared of each individual mode A increases linearly withT; any deviations &om harmonic behavior therefore isrevealed by considering the quantity P e

——(u e(T))/T.The temperature derivative of P e is determined from uqand eg e by a simple difFerentiation of (Al):

dP e

d,Tka - I d~x le~ el dies el

l

—2m „( dT cu~s dT (v~2 )

By introducing the following weighted average over thephonon modes and Cartesian directions:

monicity of the interatomic potential. The thermally-averaged positions are then determined kom minimizingthe corresponding &ee energy. The frequency (dp andthe o.th Cartesian component of the polarization vector,e~ ~, of an atom in layer 8 for a phonon mode A will thenvary with T. In this approximation, the mean-squaredisplacement (u e(T))/T (in the classical regime) is stillgiven by the standard harmonic result,

This work was supported by grants &om the NaturalScience and Engineering Research Council (NSERC) ofCanada and the "Fonds pour la formation de chercheurset 1'aide a la recherche" of the Province of Quebec toL.J.L., as well as by the Swedish Natural Science Re-search Council (NFR) to M.P. Y.B. is thankful toNSERC for financial support. We are indebted to the"Services informatiques de l'Universite de Montreal" forgenerous allocations of computer resources.

APPENDIX

We present in this appendix a derivation, within aquasiharmonic approximation, of the results stated inEqs. (3.7) and (3.8) for the temperature dependence ofthe MSAV's.

In the quasiharmonic approximation. , the restoringforce on an atom is assumed to be linear in the displace-ment &om equilibrium but, in contrast to a strictly har-monic model, force constants vary parametrically withthe lattice parameter so as to account for the anhar-

the result for dP e/dT in (A2) can be directly expressedas

dPe — d Inn', d ln leg el

dT dT dTAcx

We now discuss briefly the contribution of the termpe(dl lnFp~el /dT)p to dpe/dT in (A4). In the bulk,P e = P = P~ is isotropic because of the cubic symmetry,while the translation symmetry causes P leg el = I,implying that the contribution of the above term is zero.This proves the result in Eq. (3.7). The broken cu-bic and translation symmetry at the surface makes P eanisotropic and to vary with E, so that die& el /dT is notnecessarily zero for a particular layer. However, whenconsidering a sum over all layers of Pe, the contributionsdlepoel /dT from different E and n will cancel out sincethe polarization 6elds ~p g for each individual mode A arenormalized over E and n. Such a summation of Pe over E

applied to Eq. (A4), upon subtracting the correspondingbulk terms, leads directly to Eq. (3.8).

' Author to whom correspondence should be sent. Electronicaddress: lewisphyscn. umontreal. ca

t Electronic address: tfympfy. chalmers. seY. Cao and E. Conrad, Phys. Rev. Lett. B4, 447 (1990).Y. Cao and E. Conrad, Phys. Rev. Lett. BS, 2808 (1990).H. Gobel and P. von Blanckenhagen, Phys. Rev. B O'T, 2378(1993).E.H. Conrad, L.R. Allen, D.L. Blanchard, and T. Engel,Surf. Sci. 184, 227 (1987).

G. Armand, D. Gorse, J. Lapujoulade, and R.J. Manson,Europhys. Lett. 3, 1113 (1987).P. Zeppenfeld, K. Kern, R. David, K. Kuhnke, and G.Comsa, Phys. Rev. B 38, 12 329 (1988); G. Armand and P.Zeppenfeld, ibid 40, 5936 (1989); P. Zepp. enfeld, K. Kern,R. David, and G. Comsa, Phys. Rev. Lett. B2, 63 (1989).H.-J. Ernst, E. Hulpke, and J.P. Toennies, Phys. Rev. 84B, 16 081 (1992).J.%'.M. Frenken, F. Huusen, and J.F. van der Veen, Phys.

50 SURFACE ANHARMONICITIES AND DISORDERING ON. . . 12 103

Rev. Lett. 58, 401 (1987).H. Durr, R. Schneider, and Th. Fauster, Phys. Rev. B 43,12 187 (1991).A.P. Baddorf and E.W. Plummer, Phys. Rev. Lett. 66,2770 (1991);J. Electron Spectrosc. Relat. Phenom. 54/55,541 (1990).M.S. Daw and M.I. Baskes, Phys. Rev. B 29, 6443 (1984).S.M. Foiles, M.I. Baskes, and M.S. Daw, Phys. Rev. B 33,7983 (1986).For a comprehensive review of the RAM theory and appli-cations, see M.S. Daw, S.M. Foiles, and M.I. Baskes, Mater.Sci. Rep. 9, 251 (1993); the case of surfaces and interfacesis addressed by S.M. Foiles, in Equilibrium Structure andProperties of Surfaces and Interfaces, edited by A. Gonisand G.M. Stocks (Plenum, New York, 1992), p. 89.J.K. N@rskov, Phys. Rev. B 26, 2875 (1982); K.W. Jacob-son, J.K. N@rskov, and M.J. Puska, ibid. 35, 7423 (1987).E.T. Chen, R.N. Barnett, and U. Landman, Phys. Rev. B41, 439 (1990).E.T. Chen, R.N. Barnett, and U. Landman, Phys. Rev. B40, 924 (1989).R.N. Barnett and U. Landman, Phys. Rev. B 44, 3226(1991).P. Stoltze, J.K. N@rskov, and U. Landman, Phys. Rev. Lett.61, 220 (1988); Surf. Sci. Lett. 220, L693 (1989);P. Stoltze,J. Chem. Phys. 92, 6306 (1990).H. Hakkinen and M. Manninen, Phys. Rev. B 46, 1725(1992).L. Yang, T.S. Rahman, and M.S. Daw, Phys. Rev. B 44,13725 (1991).L. Yang and T.S. Rahman, Phys. Rev. Lett. 67, 2327(1991).P.D. Ditlevsen, P. Stoltze, and J.K. N@rskov, Phys. Rev. B44, 13002 (1991).T.S. Rahman and Z.J. Tian, J. Electron Spectrosc. Relat.Phenom. 64/65, 651 (1993).Y. Beaudet, L. Lewis, and M. Persson, Phys. Rev. B 47,4127 (1993).Y. Beaudet, L. Lewis, and M. Persson, in Material8 Theoryand Modeling, edited by J. Broughton, P. Bristowe, and J.Newsam, MRS Symposia Proceedings No. 291 (MaterialsResearch Society, Pittsburgh, 1993), p. 525.Our unit vectors x, y, and z are related to the unit vectorse, e„, and e of the conventional unit cell for the fcc lattice(the simple cubic lattice with a four-atom basis) as follows:For the (100) surface, x = (e +e„)/v 2, y = (e —e„)/J2,and z = e, ; for the (110) surface, z = (e —e„)/v 2, y = e„and z = (e + e„)/v 2.G. Boisvert (private communication).For a review of MD technology, see, for instance, M.P.Allen and D.J. Tildesley, Computer Simulation of Liquids(Clarendon, Oxford, 1987).S.M. Foiles and J.B.Adams, Phys. Rev. B 40, 5909 (1989).R.J. Birgenau, J. Cordes, G. Dolling, and A.D.B. Woods,Phys. Rev. 136, A1359 (1964).Y.S. Touloukian, R.K. Kirby, R.E. Taylor, and P.D.Desai, Therrnophysical Properties of Matter, ThermalEzpansion Metallic Elements —and Alloys (Plenum, NewYork, 1975), Vol. 12.G. Grimvall, Thermophysical Properties of Materials, Se-

lected Topics in Solid State Physics (North-Holland, Ams-terdam, 1986), Vol. 28.H. Hakkinen, J. Merikoski, M. Manninen, J. Timonen, andK. Kaski, Phys. Rev. Lett. 70, 2451 (1993); J. Merikoski,H. Hakkinen, M. Manninen, J. Timonen, and K. Kaski,Phys. Rev. B 49, 4938 (1994).See, for instance, J.-P. Hansen and I.R. McDonald, Theoryof Simple Liquids (Academic, London, 1986).G.P. Srivastava, The Physics of Phonons (Adam Hilger,Bristol, 1990).W. Schommers, in Structure and Dynamics of Surfaces I,edited by W. Schommers and P. von Blanckenhagen, Top-ics of Current Physics Vol. 41 (Springer, Berlin, 1986), p.199; K.H. Rieder, ibid. , p. 17.I.J. Robertson, V. Heine, and M.C. Payne, Phys. Rev. Lett.70, 1944 (1993).A.A. Maradudin and P.A. Flinn, Phys. Rev. 129, 2529(1963).Q.T. Jiang, P. Fenter, and T. Gustafsson, Phys. Rev. B 44,5773 (1991).J.P. Hansen and M.L. Klein, Phys. Rev. B 13, 878 (1976).M.L. Xu, B.M. Hall, S.Y. Tong, M. Rocca, H. Ibach, S.Lehwald, and J.E. Black, Phys. Rev. Lett. 54, 1171 (1985).M. Rocca, S. Lehwald, H. Ibach, and T.S. Rahman, Surf.Sci. 171, 632 (1986).S. Lehwald, F. Wolf, H. Ibach, M. Hall, and D.L. Mills,Surf. Sci. 192, 131 (1987)M. Balden, S. Lehwald, and H. Ibach, Phys. Rev. B 46,4172 (1992).A.A. Maradudin and A.E. Fein, Phys. Rev. 128, 2589(1962).We Bnd no evidence in the widths for a crossover from T toT dependence at high temperatures, as was observed onCu(110) by Yang and Rahman (Ref. 21) using EAM andby Ditlevsen et al. (Ref. 22) using EMT; a similar effectwas observed on Ni and Cu(110) using the Finnis-Sinclairpotential by D.D. Koleske and S.J. Sibener, Surf. Sci. 298,215 (1993).This value is based on averaged phonon frequencies mea-sured at 576 and 676 K by inelastic neutron scattering, asreported in G.A. De Wit and B.N. Brockhouse, J. Appl.Phys 39, 451 .(1968).G. Benedek and J.P. Toennies, Phys. Rev. B. 46, 13643(1992).M. Gester, D. Kleinhesselink, P. Ruggerone, and J.P. Toen-nies, Phys. Rev. B. 49, 5777 (1994).S.W. Lovesey, Theory of Neutron Scattering from Con-densed Matter (Clarendon Press, Oxford, 1984).' J.B. Pendry, Low Energy Electron Digraction (AcademicPress, London, 1974).In addition to the two-dimensional direct-lattice basis vec-tors defined in Fig. 1, we have as = (a/2~2)(x —y+ ~2z)on the (100) surface, and as ——(a/2~2)( —z + v 2y + z) onthe (110) surface. The three-dimensional reciprocal-latticebasis vectors, therefore, are given by bz = (2x/a)(~2i z), —b2 = (2z'/a)(~2y+ z), and bs = (4z/a)z for the (100) sur-face, and bq ——(2z'~2/a)(x + z), bz = (2z./a)(y —~2z),and bs = (4v 2z'/a)z for the (110) surface.A.A. Maradudin, E.W. Montroll, and G.H. Weiss, SohdState Physics (Academic Press, New York, 1963), SuppL 3.